/***********************************************************
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