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