Cで平方根を計算してみた
昨日始めたばかりのPythonで√５をうまく計算できなかったので、久しぶりにC言語でもやってみた。
ソースコードはこちら
/* main.c */
#include <stdio.h>
double root_digit(double square, double first, double step) {
double r = first;
for (int i = 1; i < 10; i++) {
printf(" r = %1.10f step = %1.10f (r + step)^2 = %1.10f \n", r, step, ((r + step) * (r + step)));
if (square < ((r + step) * (r + step))) {
return r;
}
r = r + step;
}
}
double root(double square) {
double r = 0.0;
double step = 1.0;
double t = 0.0;
for (int i = 1; i < 15; i++) {
t = root_digit(square, r, step);
if (t != NULL) {
r = t;
step = step / 10;
}
}
return (r);
}
void main() {
double root2 = 0.0;
for (double i = 2.0; i < 10.0; i++) {
printf("Root(%d) = %1.10f\n", i, root(i));
}
}極力、元となるPythonのコードと対比するように組んでみました。
実行結果
// Root(2)からRoot(4)までの結果は省略
r = 0.0000000000 step = 1.0000000000 (r + step)^2 = 1.0000000000
r = 1.0000000000 step = 1.0000000000 (r + step)^2 = 4.0000000000
r = 2.0000000000 step = 1.0000000000 (r + step)^2 = 9.0000000000
r = 2.0000000000 step = 0.1000000000 (r + step)^2 = 4.4100000000
r = 2.1000000000 step = 0.1000000000 (r + step)^2 = 4.8400000000
r = 2.2000000000 step = 0.1000000000 (r + step)^2 = 5.2900000000
r = 2.2000000000 step = 0.0100000000 (r + step)^2 = 4.8841000000
r = 2.2100000000 step = 0.0100000000 (r + step)^2 = 4.9284000000
r = 2.2200000000 step = 0.0100000000 (r + step)^2 = 4.9729000000
r = 2.2300000000 step = 0.0100000000 (r + step)^2 = 5.0176000000
r = 2.2300000000 step = 0.0010000000 (r + step)^2 = 4.9773610000
r = 2.2310000000 step = 0.0010000000 (r + step)^2 = 4.9818240000
r = 2.2320000000 step = 0.0010000000 (r + step)^2 = 4.9862890000
r = 2.2330000000 step = 0.0010000000 (r + step)^2 = 4.9907560000
r = 2.2340000000 step = 0.0010000000 (r + step)^2 = 4.9952250000
r = 2.2350000000 step = 0.0010000000 (r + step)^2 = 4.9996960000
r = 2.2360000000 step = 0.0010000000 (r + step)^2 = 5.0041690000
r = 2.2360000000 step = 0.0001000000 (r + step)^2 = 5.0001432100
r = 2.2360000000 step = 0.0000100000 (r + step)^2 = 4.9997407201
r = 2.2360100000 step = 0.0000100000 (r + step)^2 = 4.9997854404
r = 2.2360200000 step = 0.0000100000 (r + step)^2 = 4.9998301609
r = 2.2360300000 step = 0.0000100000 (r + step)^2 = 4.9998748816
r = 2.2360400000 step = 0.0000100000 (r + step)^2 = 4.9999196025
r = 2.2360500000 step = 0.0000100000 (r + step)^2 = 4.9999643236
r = 2.2360600000 step = 0.0000100000 (r + step)^2 = 5.0000090449
r = 2.2360600000 step = 0.0000010000 (r + step)^2 = 4.9999687957
r = 2.2360610000 step = 0.0000010000 (r + step)^2 = 4.9999732678
r = 2.2360620000 step = 0.0000010000 (r + step)^2 = 4.9999777400
r = 2.2360630000 step = 0.0000010000 (r + step)^2 = 4.9999822121
r = 2.2360640000 step = 0.0000010000 (r + step)^2 = 4.9999866842
r = 2.2360650000 step = 0.0000010000 (r + step)^2 = 4.9999911564
r = 2.2360660000 step = 0.0000010000 (r + step)^2 = 4.9999956285
r = 2.2360670000 step = 0.0000010000 (r + step)^2 = 5.0000001006
r = 2.2360670000 step = 0.0000001000 (r + step)^2 = 4.9999960757
r = 2.2360671000 step = 0.0000001000 (r + step)^2 = 4.9999965229
r = 2.2360672000 step = 0.0000001000 (r + step)^2 = 4.9999969701
r = 2.2360673000 step = 0.0000001000 (r + step)^2 = 4.9999974173
r = 2.2360674000 step = 0.0000001000 (r + step)^2 = 4.9999978646
r = 2.2360675000 step = 0.0000001000 (r + step)^2 = 4.9999983118
r = 2.2360676000 step = 0.0000001000 (r + step)^2 = 4.9999987590
r = 2.2360677000 step = 0.0000001000 (r + step)^2 = 4.9999992062
r = 2.2360678000 step = 0.0000001000 (r + step)^2 = 4.9999996534
}
^
main.c:17:0 関数の戻り値が間違っています。戻り値の型を確認してください。 doubleやはり同じあたりでエラーになってる。
どうも、PythonだからとかCだからという理由では無いような気がする。
√6の計算でも同じような結果。
二乗する計算を、pow関数を使ってみても同じ。
うーん、わからん。
GW中にPythonを動かしてみようという趣旨から、横道にそれ過ぎた感があるので、この件はそろそろ打ち止めにしようかな、、、。
追記 - Ayumiさんからのアドバイスで問題解決！
凄いすごい！一気に解決です！
/* main.c */
#include <stdio.h>
#include <math.h>
double root_digit(double square, double first, double step) {
double r = first;
for (int i = 1; i < 20; i++) {
if (square < pow((r + step), 2)) {
return r;
}
r += step;
}
}
double root(double square) {
double r = 0.0;
double step = 1.0;
double t = 0.0;
for (int i = 1; i < 15; i++) {
t = root_digit(square, r, step);
if (t != NULL) {
r = t;
step /= 10;
}
}
return (r);
}
void main() {
double root2 = 0.0;
for (double i = 2.0; i < 10.0; i++) {
printf("Root(%d) = %1.10f\n", i, root(i));
}
}root_digitがLoopする回数を20にしたら、解決！
Root(2) = 1.4142135624
Root(3) = 1.7320508076
Root(4) = 2.0000000000
Root(5) = 2.2360679775
Root(6) = 2.4494897428
Root(7) = 2.6457513111
Root(8) = 2.8284271247
Root(9) = 3.0000000000計算結果もあっているみたい。
計算途中結果をダンプして眺めてみていたのですが、２進数で小数点以下を表現することによる誤差のため、本来0であるべき桁に9が入っていたりして、そのため各桁で最大10回テストしないといけないのに最大9回しかできていなかったので平方根の計算が収束していなかったことが判明。
このため、root_digitがLoopする回数を20とすることでうまく収束したみたい。回数を11としても問題なく計算できた。（ただし、√99まで。√100をするには、回数を11よりも大きくしないといけない。）
多分、この辺りが上述のエラーの原因かと思われます。
ふー、スッキリ！
蛇足編 - Ayumiさんの記事を参考に、再帰呼び出しに変更
上記のC言語のコードを、AyumiさんのPythonのコードを参考にしながら再帰呼び出しに変えたら、すごくスッキリ。root_digitを何回Loopさせるのが妥当かなんて考える必要が無いのがいい。
√1000まで問題なく計算できることを検証済み。
正直、わたしは再帰呼び出しコードを書くのは苦手。
諸先輩のコードを見れば理解できるけど、スラスラと自らは書けない。
再帰呼び出しコードをスラスラ書く人は、絶対に私とは異なる思考回路をお持ちだと思う。
/* main.c */
#include <stdio.h>
#include <math.h>
double root_digit(double square, double r, double step) {
if (square <= pow((r + step), 2)) {
return r;
}
return root_digit(square, r + step, step);
}
double root(double square, double r, double step, int n) {
r = root_digit(square, r, step);
if (n <= 0) {
return r;
}
return root(square, r, step / 10, n - 1);
}
void main() {
double root2 = 0.0;
for (double i = 2.0; i < 10.0; i++) {
printf("Root(%d) = %1.10f\n", i, root(i, 0.0, 1.0, 10));
}
}再帰呼び出しで、いつも心配してしまうのは、最終的に解が定まって、呼び出した側に戻って来れるかということ。なるほど、このロジックなら再帰呼び出しの上限を10回までと予め決めておけるので安心。
ここまで読んでいただき、ありがとうございました。
これまでの収益は全て、それを必要としておられる方々へ、支援機関を通して寄付させていただきました。この活動は今後も継続したいと思っています。引き続きよろしくお願いいたします。