/***********************************************************
cardano.c -- 3次方程式
***********************************************************/
#N88BASIC
CONST PI = 3.14159265358979323846264 /* 円周率 */
#define CHECK
Declare Function acos CDECL Lib"msvcrt"(x As Double) As Double
Declare Function fabs CDECL Lib"msvcrt"(x As Double) As Double
Function cuberoot(x As Double) As Double
Dim s, prev
Dim pos As Long
If (x = 0) Then cuberoot = 0:Exit Function
If (x > 0) Then pos = 1 Else pos = 0: x = -x
If (x > 1) Then s = x Else s = 1
Do
prev = s: s = (x / (s * s) + 2 * s) / 3
Loop while (s < prev)
If (pos) Then cuberoot = prev Else cuberoot = -prev
End Function
Sub cardano(a, b, c, d)
Dim p, q, t, a3, b3, x1, x2, x3
b = b / (3 * a): c = c / a: d = d / a
p = b * b - c / 3
q = (b * (c - 2 * b * b) - d) / 2
a = q * q - p * p * p
debug
If (a = 0) Then
q = cuberoot(q): x1 = 2 * q - b: x2 = -q - b
Print "x = ";x1;", ";x2;" (重解)"
#ifdef CHECK
Print "f(x1) = ", x1 * (x1 * (x1 + 3 * b) + c) + d
Print "f(x2) = ", x2 * (x2 * (x2 + 3 * b) + c) + d
#endIf
Else If (a > 0) Then
If (q > 0) Then
a3 = cuberoot(q + Sqr(a))
Else
a3 = cuberoot(q - Sqr(a))
End If
b3 = p / a3
x1 = a3 + b3 - b: x2 = -0.5 * (a3 + b3) - b
x3 = fabs(a3 - b3) * Sqr(3.0) / 2
Print "x = ";x1;"; ";x2;" +- ";x3;" i"
#ifdef CHECK
Print "f(x1) = "; x1 * (x1 * (x1 + 3 * b) + c) + d
Print "f(x2+-x3i) = ";((x2+3*b)*x2-x3*x3+c)*x2-(2*x2+3*b)*x3*x3+d;" +- ";((3*x2+6*b)*x2-x3*x3+c)*x3;" i"
#endIf
Else
a = Sqr(p): t = acos(q / (p * a)): a *= 2
x1 = a * Cos( t / 3) - b
x2 = a * Cos((t + 2 * PI) / 3) - b
x3 = a * Cos((t + 4 * PI) / 3) - b
Print "x = ";x1, x2, x3
#ifdef CHECK
Print "f(x1) = ", x1 * (x1 * (x1 + 3 * b) + c) + d
Print "f(x2) = ", x2 * (x2 * (x2 + 3 * b) + c) + d
Print "f(x3) = ", x3 * (x3 * (x3 + 3 * b) + c) + d
#endIf
End If
End Sub
Dim a, b, c, d
Input "a b c d = ";a,b,c,d
If (a = 0) Then
Print "3次方程式ではありません."
Else
cardano(a, b, c, d)
End If
