C言語で一次方程式，二次方程式，三次方程式，四次方程式の解の公式を教えて！
こういった悩みにお答えします．
本記事の信頼性
- リアルタイムシステムの研究歴12年．
- 東大教員の時に，英語でOS（Linuxカーネル）の授業．
- 2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校（UNC）コンピュータサイエンス学部で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
- プログラミング歴15年以上，習得している言語: C/C++，Python，Solidity/Vyper，Java，Ruby，Go，Rust，D，HTML/CSS/JS/PHP，MATLAB，Verse（UEFN）, Assembler (x64，ARM)．
- 東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．
- 2020年1月～現在はアメリカのノースカロライナ州チャペルヒルにあるGuarantee Happiness LLCのCTOとしてECサイト開発やWeb/SNSマーケティングの業務．2022年6月～現在はアメリカのノースカロライナ州チャペルヒルにあるJapanese Tar Heel, Inc.のCEO兼CTO．
- 最近は自然言語処理AIとイーサリアムに関する有益な情報発信に従事．
- （AI全般を含む）自然言語処理AIの論文の日本語訳や，AIチャットボット（ChatGPT，Auto-GPT，Gemini（旧Bard）など）の記事を50本以上執筆．アメリカのサンフランシスコ（広義のシリコンバレー）の会社でプロンプトエンジニア・マネージャー・Quality Assurance（QA）の業務委託の経験あり．
- （スマートコントラクトのプログラミングを含む）イーサリアムや仮想通貨全般の記事を200本以上執筆．イギリスのロンドンの会社で仮想通貨の英語の記事を日本語に翻訳する業務委託の経験あり．
こういった私から学べます．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
目次
C言語で一次方程式，二次方程式，三次方程式，四次方程式の解の公式
C言語で一次方程式，二次方程式，三次方程式，四次方程式の解の公式を紹介します．
※五次以上の方程式の解の公式がないこと（代数的に解く方法がないこと）は，アーベル・ルフィニの定理により証明されています．
五次以上の方程式を代数的ではなく超越的に解く方法を知りたいあなたはこちらの動画を観ましょう．
本記事では複素数を理解していることを前提としています．
一次方程式の解の公式
一次方程式の解の公式を紹介します．
\begin{eqnarray*}
ax + b &=& 0\ (a \neq 0) \\
x &=& -\frac{b}{a}
\end{eqnarray*}
一次方程式の解の公式のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <float.h> void linear_equation_formula(double *x, double a, double b) { if (abs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } *x = -b / a; } int main(void) { double a, b, x; printf("ax + b = 0 (a != 0)\n"); printf("Please input a and b: "); scanf("%lf%lf", &a, &b); linear_equation_formula(&x, a, b); printf("x = %lf\n", x); return 0; } |
実行結果は以下になります．
こちらの結果と同じになることを確認しましょう．
|
1 2 3 4 5 6 7 8 9 |
$ gcc linear_equation_formula.c $ a.out ax + b = 0 (a != 0) Please input a and b: 2 -2 x = 1.000000 $ a.out ax + b = 0 (a != 0) Please input a and b: 0 0 Error: a == 0 |
ニ次方程式の解の公式
ニ次方程式の解の公式を紹介します．
\begin{eqnarray*}
ax^2 + bx + c &=& 0\ (a \neq 0) \\
x &=& \frac{-b \pm{\sqrt{b^2 - 4ac}}}{2a}
\end{eqnarray*}
ニ次方程式の解の公式の実装で利用するsqrt関数を知りたいあなたはこちらからどうぞ．
ニ次方程式の判別式Dは以下になります．Dが負の場合は複素数の解を持ちます．
$$D = b^2 - 4ac$$
ニ次方程式の解の公式で実数解のみを持つ場合のコード
ニ次方程式の解の公式で実数解のみを持つ場合のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <float.h> #include <math.h> #define N 2 void quadratic_equation_formula(double x[N], double a, double b, double c) { double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } dis = b * b - 4 * a * c; if (dis < 0.0) { fprintf(stderr, "Error: D < 0.0\n"); exit(2); } x[0] = (-b - sqrt(dis)) / (2 * a); x[1] = (-b + sqrt(dis)) / (2 * a); } int main(void) { double a, b, c, x[N]; int i; printf("ax^2 + bx + c = 0 (a != 0)\n"); printf("Please input a, b, and c: "); scanf("%lf%lf%lf", &a, &b, &c); quadratic_equation_formula(x, a, b, c); for (i = 0; i < N; i++) { printf("x%d = %lf\n", i + 1, x[i]); } return 0; } |
実行結果は以下になります．
こちらと同じ結果になることを確認しましょう．（実数解のみを持つ場合は同じ結果になります．）
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
$ gcc quadratic_equation_formula.c -lm $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 2 3 -5 x1 = -2.500000 x2 = 1.000000 $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 1 2 1 x1 = -1.000000 x2 = -1.000000 $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 1 1 1 Error: D < 0.0 |
二次方程式の解の公式で複素数の解を持つ場合のコード
二次方程式の解の公式で複素数の解を持つ場合のコードは以下になります．
C99規格で採用されたcomplex型（複素数型）を利用しています．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <float.h> #include <math.h> #include <complex.h> #define N 2 void quadratic_equation_formula(double complex x[N], double a, double b, double c) { double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } dis = b * b - 4 * a * c; x[0] = (-b - csqrt(dis)) / (2 * a); x[1] = (-b + csqrt(dis)) / (2 * a); } int main(void) { double a, b, c; double complex x[N]; int i; printf("ax^2 + bx + c = 0 (a != 0)\n"); printf("Please input a, b, and c: "); scanf("%lf%lf%lf", &a, &b, &c); quadratic_equation_formula(x, a, b, c); for (i = 0; i < N; i++) { printf("x%d = %lf", i + 1, creal(x[i])); if (cabsf(cimag(x[i])) >= DBL_EPSILON) { printf("%+lfi", cimag(x[i])); } printf("\n"); } return 0; } |
実行結果は以下になります．
複素数の解を持っていることがわかります．
こちらと同じ結果になることを確認しましょう．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
$ gcc quadratic_equation_formula2.c -lm $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 2 3 -5 x1 = -2.500000 x2 = 1.000000 $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 1 2 1 x1 = -1.000000 x2 = -1.000000 $ a.out ax^2 + bx + c = 0 (a != 0) Please input a, b, and c: 1 1 1 x1 = -0.500000-0.866025i x2 = -0.500000+0.866025i |
三次方程式の解の公式（カルダノの公式）
三次方程式の解の公式（カルダノの公式）を紹介します．
\begin{eqnarray*}
ax^3 + bx^2 + cx + d &=& 0\ (a \neq 0) \\
\end{eqnarray*}
ここで，\(A_1\)，\(A_2\)，\(A_3\)，p，qを以下とします．
\begin{eqnarray*}
A_1 &=& \frac{b}{3a} \\
A_2 &=& \frac{c}{a} \\
A_3 &=& \frac{d}{a} \\
p &=& A_1^2 - \frac{A_2}{3} \\
q &=& \frac{A_1(A_2 - 2A_1^2) - A_3}{2}
\end{eqnarray*}
判別式Dは，p，qを利用すると以下になります．
\begin{eqnarray*}
D = p * p * p - q * q
\end{eqnarray*}
D > 0の場合，xの実数解\(x_1\)，\(x_2\)，\(x_3\)を\(A_1\)，p，qを利用して表すと以下になります．
\begin{eqnarray*}
x_1 &=& 2\sqrt{p} \cos(\frac{\arccos \frac{q}{p\sqrt{p}}}{3}) - A_1 \\
x_2 &=& 2\sqrt{p} \cos(\frac{\arccos \frac{q}{p\sqrt{p}} + 2\pi}{3}) - A_1 \\
x_3 &=& 2\sqrt{p} \cos(\frac{\arccos \frac{q}{p\sqrt{p}} + 4\pi}{3}) - A_1 \\
\end{eqnarray*}
D = 0の場合，xの実数解\(x_1\)，\(x_2\)，\(x_3\)を\(A_1\)，qを利用して表すと以下になります．（\(x_2\)と\(x_3\)は重解です．）
\begin{eqnarray*}
x_1 &=& 2\sqrt[3]{q} - A_1 \\
x_2 &=& -\sqrt[3]{q} - A_1 \\
x_3 &=& -\sqrt[3]{q} - A_1
\end{eqnarray*}
D < 0の場合，Qを以下のように定義します．
\begin{eqnarray*}
Q &=& \left\{
\begin{array}{ll}
\sqrt[3]{q + \sqrt{-D}} & (q \geq 0)\\
\sqrt[3]{q - \sqrt{-D}} & (q < 0)
\end{array}
\right.
\end{eqnarray*}
そして，xの実数解\(x_1\)，xの2つの複素数解\(x_2\)と\(x_3\)を\(A_1\)，p，Qを利用して表すと以下になります．
\begin{eqnarray*}
x_1 &=& Q + \frac{p}{Q} - A_1 \\
x_2 &=& -\frac{Q + \frac{p}{Q}}{2} - A_1 + \frac{|Q - \frac{p}{Q}| \sqrt{3})}{2}i \\
x_3 &=& -\frac{Q + \frac{p}{Q}}{2} - A_1 - \frac{|Q - \frac{p}{Q}| \sqrt{3})}{2}i \\
\end{eqnarray*}
三次方程式の解の公式で実数解のみを持つ場合のコード
三次方程式の解の公式で実数解のみを持つ場合のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <float.h> #include <math.h> #define PI 3.1415926535897932384626433832795 #define N 3 double get_cubic_root(double x) { if (x > 0.0) { return pow(x, 1.0 / 3.0); } else { return -pow(-x, 1.0 / 3.0); } } void cubic_equation_formula(double x[N], double a, double b, double c, double d) { double p, q, r, t; double a1, a2, a3; double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / (3 * a); a2 = c / a; a3 = d / a; p = a1 * a1 - a2 / 3; q = (a1 * (a2 - 2 * a1 * a1) - a3) / 2; dis = p * p * p - q * q; if (dis < 0.0) { fprintf(stderr, "Error: D < 0.0\n"); exit(2); } else if (fabs(dis) <= DBL_EPSILON) { r = get_cubic_root(q); x[0] = 2 * r - a1; x[1] = x[2] = -r - a1; } else { r = sqrt(p); t = acos(q / (p * r)); r *= 2; x[0] = r * cos(t / 3) - a1; x[1] = r * cos((t + 2 * PI) / 3) - a1; x[2] = r * cos((t + 4 * PI) / 3) - a1; } } int main(void) { double a, b, c, d; double x[N]; int i; printf("ax^3 + bx^2 + cx + d = 0 (a != 0)\n"); printf("Please input a, b, c, and d: "); scanf("%lf%lf%lf%lf", &a, &b, &c, &d); cubic_equation_formula(x, a, b, c, d); for (i = 0; i < N; i++) { printf("x%d = %lf\n", i + 1, x[i]); } return 0; } |
実行結果は以下になります．
こちらの結果と同じになることを確認しましょう．（実数解のみを持つ場合は同じ結果になります．）
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
$ gcc cubic_equation_formula.c -lm $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 -2 -11 12 x1 = 4.000000 x2 = -3.000000 x3 = 1.000000 $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 -7 15 -9 x1 = 3.000000 x2 = 1.000000 x3 = 3.000000 $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 2 3 4 Error: D < 0.0 |
三次方程式の解の公式で複素数の解を持つ場合のコード
三次方程式の解の公式で複素数の解を持つ場合のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <float.h> #include <math.h> #include <complex.h> #define PI 3.1415926535897932384626433832795 #define N 3 double get_cubic_root(double x) { if (x > 0.0) { return pow(x, 1.0 / 3.0); } else { return -pow(-x, 1.0 / 3.0); } } void cubic_equation_formula(double complex x[N], double a, double b, double c, double d) { double p, q, r, t; double a1, a2, a3; double pp, qq; double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / (3 * a); a2 = c / a; a3 = d / a; p = a1 * a1 - a2 / 3; q = (a1 * (a2 - 2 * a1 * a1) - a3) / 2; dis = p * p * p - q * q; if (dis < 0.0) { if (q >= 0.0) { qq = get_cubic_root(q + sqrt(-dis)); } else { qq = get_cubic_root(q - sqrt(-dis)); } pp = p / qq; r = -(qq + pp) / 2 - a1; t = (fabs(qq - pp) * sqrt(3.0)) / 2; x[0] = qq + pp - a1; x[1] = CMPLX(r, t); x[2] = CMPLX(r, -t); } else if (fabs(dis) <= DBL_EPSILON) { r = get_cubic_root(q); x[0] = 2 * r - a1; x[1] = x[2] = -r - a1; } else { r = sqrt(p); t = acos(q / (p * r)); r *= 2; x[0] = r * cos(t / 3) - a1; x[1] = r * cos((t + 2 * PI) / 3) - a1; x[2] = r * cos((t + 4 * PI) / 3) - a1; } } int main(void) { double a, b, c, d; double complex x[N]; int i; printf("ax^3 + bx^2 + cx + d = 0 (a != 0)\n"); printf("Please input a, b, c, and d: "); scanf("%lf%lf%lf%lf", &a, &b, &c, &d); cubic_equation_formula(x, a, b, c, d); for (i = 0; i < N; i++) { printf("x%d = %lf", i + 1, creal(x[i])); if (cabsf(cimag(x[i])) >= DBL_EPSILON) { printf("%+lfi", cimag(x[i])); } printf("\n"); } return 0; } |
実行結果は以下になります．
こちらの結果と同じになることを確認しましょう．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
$ gcc cubic_equation_formula2.c -lm $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 -2 -11 12 x1 = 4.000000 x2 = -3.000000 x3 = 1.000000 $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 -7 15 -9 x1 = 3.000000 x2 = 1.000000 x3 = 3.000000 $ a.out ax^3 + bx^2 + cx + d = 0 (a != 0) Please input a, b, c, and d: 1 2 3 4 x1 = -1.650629 x2 = -0.174685+1.546869i x3 = -0.174685-1.546869i |
四次方程式の解の公式（フェラーリの公式）
四次方程式の解の公式（フェラーリの公式）を紹介します．
\begin{eqnarray*}
ax^4 + bx^3 + cx^2 + dx + e &=& 0\ (a \neq 0) \\
\end{eqnarray*}
ここで，\(A_1\)，\(A_2\)，\(A_3\)，\(A_4\)，p，q，rを以下とします．
\begin{eqnarray*}
A_1 &=& \frac{b}{a} \\
A_2 &=& \frac{c}{a} \\
A_3 &=& \frac{d}{a} \\
A_4 &=& \frac{e}{a} \\
p &=& A_2 - 6\left(\frac{A_1}{4}\right)^2 \\
q &=& A_3 - 2A_2\left(\frac{A_1}{4}\right) + 8\left(\frac{A_1}{4}\right)^3 \\
r &=& A_4 - A_3 \left(\frac{A_1}{4}\right) + A_2 \left(\frac{A_1}{4}\right)^2 -3 \left(\frac{A_1}{4}\right)^3
\end{eqnarray*}
また，与えられた四次方程式に関するuの三次分解方程式は以下になります．
\begin{eqnarray*}
u (p + u)^2 - 4ru = q^2
\end{eqnarray*}
\(m = \sqrt{u - p}\)とする場合，xの解は以下になります．
\begin{eqnarray*}
x_1 &=& \frac{1}{2}\left(-m + \sqrt{- u - p + \frac{2q}{m}} \right) - \frac{A_1}{4} \\
x_2 &=& \frac{1}{2}\left(-m - \sqrt{- u - p + \frac{2q}{m}} \right) - \frac{A_1}{4} \\
x_3 &=& \frac{1}{2}\left(m + \sqrt{- u - p - \frac{2q}{m}} \right) - \frac{A_1}{4} \\
x_4 &=& \frac{1}{2}\left(m - \sqrt{- u - p - \frac{2q}{m}} \right) - \frac{A_1}{4} \\
\end{eqnarray*}
四次方程式の判別式はとても複雑なので，学びたいあなたは解の様子を読みましょう．
四次方程式の解の公式で実数解のみを持つ場合のコード
四次方程式の解の公式で実数解のみを持つ場合のコードは以下になります．
三次方程式の解の公式のコードのget_cubic_root関数やcubic_equation_formula関数を利用しています．
また，are_all_roots_real関数では，四次方程式の判別式で実数解を持つ場合のみを判定しています．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <stdbool.h> #include <float.h> #include <math.h> #define PI 3.1415926535897932384626433832795 #define N 4 double get_cubic_root(double x) { if (x > 0.0) { return pow(x, 1.0 / 3.0); } else { return -pow(-x, 1.0 / 3.0); } } void cubic_equation_formula(double x[N], double a, double b, double c, double d) { double p, q, r, t; double a1, a2, a3; double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / (3 * a); a2 = c / a; a3 = d / a; p = a1 * a1 - a2 / 3; q = (a1 * (a2 - 2 * a1 * a1) - a3) / 2; dis = p * p * p - q * q; if (dis < 0.0) { fprintf(stderr, "Error: D < 0.0\n"); exit(2); } else if (fabs(dis) <= DBL_EPSILON) { r = get_cubic_root(q); x[0] = 2 * r - a1; x[1] = x[2] = -r - a1; } else { r = sqrt(p); t = acos(q / (p * r)); r *= 2; x[0] = r * cos(t / 3) - a1; x[1] = r * cos((t + 2 * PI) / 3) - a1; x[2] = r * cos((t + 4 * PI) / 3) - a1; } } bool are_all_roots_real(double a, double b, double c, double d, double e) { double dis, pp, rr, delta0, dd; dis = 256 * pow(a * e, 3) - 192 * pow(a * e, 2) * b * d - 128 * pow(a * c * e, 2) + 144 * pow(a * d, 2) * c * e - 27 * pow(a * d * d, 2) + 144 * pow(b * e, 2) * a * c - 6 * pow(b * d, 2) * a * e - 80 * pow(c, 2) * a * b * d * e + 18 * a * b * c * pow(d, 3) + 16 * pow(c, 4) * a * e - 4 * a * pow(c, 3) * pow(d, 2) - 27 * pow(b * b * e, 2) + 18 * pow(b, 3) * c * d * e - 4 * pow(b * d, 3) - 4 * pow(b, 2) * pow(c, 3) * e + pow(b * c * d, 2); pp = 8 * a * c - 3 * pow(b, 2); rr = pow(b, 3) + 8 * pow(a, 2) * d - 4 * a * b * c; delta0 = pow(c, 2) - 3 * b * d + 12 * a * e; dd = 64 * pow(a, 3) * e - 16 * pow(a * c, 2) + 16 * a * pow(b, 2) * c - 16 * pow(a, 2) * b * d - 3 * pow(b, 4); if (dis < 0.0) { /* 2 real roots and 2 (1-pair) complex roots. */ return false; } else if (fabs(dis) <= DBL_EPSILON) { if (pp < 0.0 && dd < 0.0 && fabs(delta0) > DBL_EPSILON) { /* 1 real double root and 2 real roots. */ return true; } else if (dd > 0.0 || (pp < 0.0 && (fabs(dd) > DBL_EPSILON || fabs(rr) > DBL_EPSILON))) { /* 1 real double root and 2 (1-pair) complex roots. */ return false; } else if (fabs(delta0) <= DBL_EPSILON && fabs(dd) > DBL_EPSILON) { /* 1 real triple root and 1 real root. */ return true; } else if (fabs(dd) <= DBL_EPSILON) { if (pp < 0.0) { /* 2 real double roots. */ return true; } else if (pp < 0.0 && fabs(rr) <= DBL_EPSILON) { /* 2 (1-pair) complex double roots. */ return false; } else if (fabs(delta0) <= DBL_EPSILON) { /* 1 real quadruple root (= -b/4a). */ return true; } else { /* not reached but just in case. */ fprintf(stderr, "Error: not reached.\n"); exit(1); } } else { /* not reached but just in case. */ fprintf(stderr, "Error: not reached.\n"); exit(2); } } else { if (pp < 0.0 && dd < 0.0) { /* 4 real roots. */ return true; } else if (pp < 0.0 || dd > 0.0) { /* 4 (2-pairs) complex roots. */ return false; } else { /* not reached but just in case. */ fprintf(stderr, "Error: not reached.\n"); exit(3); } } /* not reached but just in case. */ fprintf(stderr, "Error: not reached.\n"); exit(4); } void quartic_equation_formula(double x[N], double a, double b, double c, double d, double e) { double a1, a2, a3, a4, p, q, r; double xx[N - 1]; double u; double m; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / a; a2 = c / a; a3 = d / a; a4 = e / a; p = a2 - 6 * pow(a1 / 4, 2); q = a3 - (2 * a2 * a1) / 4 + 8 * pow(a1 / 4, 3); r = a4 - (a3 * a1) / 4 + a2 * pow(a1 / 4, 2) - 3 * pow(a1 / 4, 4); if (are_all_roots_real(a, b, c, d, e)) { cubic_equation_formula(xx, 1, -p, -4 * r, 4 * p * r - q * q); } else { fprintf(stderr, "Error: at least one root is complex.\n"); exit(1); } u = xx[0]; m = sqrt(u - p); x[0] = (-m + sqrt(-u - p + (2 * q) / m)) / 2 - a1 / 4; x[1] = (-m - sqrt(-u - p + (2 * q) / m)) / 2 - a1 / 4; x[2] = (m + sqrt(-u - p - (2 * q) / m)) / 2 - a1 / 4; x[3] = (m - sqrt(-u - p - (2 * q) / m)) / 2 - a1 / 4; } int main(void) { double a, b, c, d, e; double x[N]; int i; printf("ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0)\n"); printf("Please input a, b, c, d, and e: "); scanf("%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e); quartic_equation_formula(x, a, b, c, d, e); for (i = 0; i < N; i++) { printf("x%d = %lf\n", i + 1, x[i]); } return 0; } |
実行結果は以下になります．
こちらの結果と同じになることを確認しましょう．（実数解のみを持つ場合は同じ結果になります．）
|
1 2 3 4 5 6 7 8 9 10 11 12 |
$ gcc quartic_equation_formula.c -lm $ a.out ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0) Please input a, b, c, d, and e: 1 -7 5 31 -30 x1 = 1.000000 x2 = -2.000000 x3 = 5.000000 x4 = 3.000000 $ a.out ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0) Please input a, b, c, d, and e: 1 2 3 4 5 Error: at least one root is complex. |
四次方程式の解の公式で複素数の解を持つ場合のコード
四次方程式の解の公式で複素数の解を持つ場合のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <stdbool.h> #include <float.h> #include <math.h> #include <complex.h> #define PI 3.1415926535897932384626433832795 #define N 4 double get_cubic_root(double x) { if (x > 0.0) { return pow(x, 1.0 / 3.0); } else { return -pow(-x, 1.0 / 3.0); } } void cubic_equation_formula(double complex x[N], double a, double b, double c, double d) { double p, q, r, t; double a1, a2, a3; double pp, qq; double dis; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / (3 * a); a2 = c / a; a3 = d / a; p = a1 * a1 - a2 / 3; q = (a1 * (a2 - 2 * a1 * a1) - a3) / 2; dis = p * p * p - q * q; if (dis < 0.0) { if (q >= 0.0) { qq = get_cubic_root(q + sqrt(-dis)); } else { qq = get_cubic_root(q - sqrt(-dis)); } pp = p / qq; r = -(qq + pp) / 2 - a1; t = (fabs(qq - pp) * sqrt(3.0)) / 2; x[0] = qq + pp - a1; x[1] = CMPLX(r, t); x[2] = CMPLX(r, -t); } else if (fabs(dis) <= DBL_EPSILON) { r = get_cubic_root(q); x[0] = 2 * r - a1; x[1] = x[2] = -r - a1; } else { r = sqrt(p); t = acos(q / (p * r)); r *= 2; x[0] = r * cos(t / 3) - a1; x[1] = r * cos((t + 2 * PI) / 3) - a1; x[2] = r * cos((t + 4 * PI) / 3) - a1; } } void quartic_equation_formula(double complex x[N], double a, double b, double c, double d, double e) { double a1, a2, a3, a4, p, q, r; double complex xx[N - 1]; double u; double complex m; if (fabs(a) <= DBL_EPSILON) { fprintf(stderr, "Error: a == 0\n"); exit(1); } a1 = b / a; a2 = c / a; a3 = d / a; a4 = e / a; p = a2 - 6 * pow(a1 / 4, 2); q = a3 - (2 * a2 * a1) / 4 + 8 * pow(a1 / 4, 3); r = a4 - (a3 * a1) / 4 + a2 * pow(a1 / 4, 2) - 3 * pow(a1 / 4, 4); cubic_equation_formula(xx, 1, -p, -4 * r, 4 * p * r - q * q); u = creal(xx[0]); m = csqrt(u - p); x[0] = (-m + csqrt(-u - p + (2 * q) / m)) / 2 - a1 / 4; x[1] = (-m - csqrt(-u - p + (2 * q) / m)) / 2 - a1 / 4; x[2] = (m + csqrt(-u - p - (2 * q) / m)) / 2 - a1 / 4; x[3] = (m - csqrt(-u - p - (2 * q) / m)) / 2 - a1 / 4; } int main(void) { double a, b, c, d, e; double complex x[N]; int i; printf("ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0)\n"); printf("Please input a, b, c, d, and e: "); scanf("%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e); quartic_equation_formula(x, a, b, c, d, e); for (i = 0; i < N; i++) { printf("x%d = %lf", i + 1, creal(x[i])); if (cabsf(cimag(x[i])) >= DBL_EPSILON) { printf("%+lfi", cimag(x[i])); } printf("\n"); } return 0; } |
実行結果は以下になります．
こちらの結果と同じになることを確認しましょう．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
$ gcc quartic_equation_formula2.c -lm $ a.out ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0) Please input a, b, c, d, and e: 1 -7 5 31 -30 x1 = 1.000000 x2 = -2.000000 x3 = 5.000000 x4 = 3.000000 $ a.out ax^4 + bx^3 + cx^2 + dx + e = 0 (a != 0) Please input a, b, c, d, and e: 1 2 3 4 5 x1 = -1.287815+0.857897i x2 = -1.287815-0.857897i x3 = 0.287815-1.416093i x4 = 0.287815+1.416093i |
まとめ
C言語で一次方程式，二次方程式，三次方程式，四次方程式の解の公式を紹介しました．
一次方程式と二次方程式は簡単なコードですが，三次方程式と四次方程式はとても複雑で難しいので，理解できるまで何度も読み直しましょう！
連立一次方程式の解法を学びたいあなたはこちらからどうぞ．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
