目的

3次方程式の解を求める。

使用法

cubic(a, b, c, d)

引数

a, b, c, d    a x^3 + b x^2 + c x +d = 0 の係数

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/cubic.R", encoding="euc-jp")

# カルダーノ法により，3次方程式の解を求める
cubic <- function(a, b, c, d)
{
        cubic.root <- function(x) return(sign(x)*abs(x)^(1/3))                       # a x^3 + b x^2 + c x +d = 0 の係数
        res <- NULL
        res$coefficients <- c(a, b, c, d)
        if (a == 0) {
                return("3次の項の係数がゼロです")
        }
        b <- b/(3*a)
        c <- c/a
        d <- d/a
        p <- b^2-c/3
        q <- (b*c-2*b^3-d)/2
        a <- q^2-p^3
        if (a == 0) {
                q <- cubic.root(q)
                res$ans <- c(2*q-b, -q-b)
        }
        else if (a > 0) {
                a3 <- cubic.root(q+sign(q)*sqrt(a))
                b3 <- p/a3
                x <- complex(real=-(a3+b3)/2-b, imaginary=abs(a3-b3)*sqrt(3)/2)
                res$ans <- c(a3+b3-b, x, Conj(x))
        }
        else {
                a <- 2*sqrt(p)
                t <- acos(q/(p*a/2))
                res$ans <- c(a*cos(t/3)-b, a*cos((t+2*pi)/3)-b, a*cos((t+4*pi)/3)-b)
        }
        class(res) <- "cubic"
        return(res)
}
# cubic クラスのオブジェクトを表示する
print.cubic <- function(obj)
{
        put0 <- function(x) return(paste(ifelse(x >= 0, "", "-"), ifelse(abs(x) == 1, "", abs(x)), sep=""))
        put1 <- function(x) return(paste(ifelse(x >= 0, "+", "-"), ifelse(abs(x) == 1, "", abs(x)), sep=""))
        put2 <- function(x) return(paste(ifelse(x >= 0, "+", "-"), abs(x), sep=""))
        printf("%sx^3%sx^2%sx%s=0\n", put0(obj$coefficients[1]), put1(obj$coefficients[2]),
                                      put1(obj$coefficients[3]), put2(obj$coefficients[4]))
        sapply(obj$ans, print)
}


使用例

> cubic(1,1,1,1)
x^3+x^2+x+1=0
[1] -1+0i
[1] 0+1i
[1] 0-1i

> cubic(1,-5,7,-3)
x^3-5x^2+7x-3=0
[1] 3+0i
[1] 1+0i
[1] 1-0i

> cubic(1,-1,2,1)
x^3-x^2+2x+1=0
[1] -0.3926468+0i
[1] 0.696323+1.43595i
[1] 0.696323-1.43595i

> cubic(0,1,1,1)
[1] "3次の項の係数がゼロです"