「3次方程式」の編集履歴（バックアップ）一覧はこちら

「3次方程式」(2010/07/17 (土) 14:05:28) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

#asciiart(blockquote){
/***********************************************************
	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  /* $\Sqr[3]{x}$ */
	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)
			/* $ax^3 + bx^2 + cx + d = 0$, $a \ne 0$ */
	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
	
}

#asciiart(blockquote){
/***********************************************************
	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
	
}