裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

二次方程式の解を求めるプログラムのレビュー

2019年08月02日 | ブログラミング

二次方程式の解を求めるプログラム例はネット上でもよく見られる。
本来(?)お手本を示すべきものであるが,問題を含むものが少なくない(多いということだ)。

結論を先に述べておく。

(0) 全部で 40 サイトを対象とした
(1) float を使っているのは 8 サイト。そのうち C によるものは 6 サイト。C は 12 サイトだが,そのうち float を使っているのは 6 サイト。C はよく使われているが,不適切なサイトも多い
(2) 桁落ちを考慮しているサイトは 8 サイト(20%)。
(3) 虚数型を使用しているのは 3 サイトに過ぎない。
(4) 虚数解の場合を対象としていないのは 17 サイト。実部虚部を表示するのは 19 サイト。

==以下は各サイトのレビュー==

●● https://webkaru.net/clang/quadratic-formula/

 は C 言語によるものだ。

いくつか問題点がある。
(1) 使用しているのが double ではなく float --- float を使うのに何の意味があるのか。float など使う癖がつくと,後々困る。
(2) 実数解を求めるとき,2 つの解を解の公式だけで求めている
    kai1 = ( -b + sqrt(discriminant) ) / (2*a);
    kai2 = ( -b - sqrt(discriminant) ) / (2*a);
これは,次のような方程式を正しく解けない。
2次方程式の定数を入力してください。
a = 1
b = 1.0000000001
c = 0.0000000001
2次方程式の解: x = 0.00, -1.00
答えは printf で書かれているので少数以下を 2 桁に指定しているのでこうなる。また,上述のように float を使っているので有効桁数は 6 桁ほどしかない。そこで,float は全て double にする。係数を端末から読み込むのに scanf("%f", &a) などと,お勧めできない書法で書かれているが,最小限の手当で,scanf("%lf", &a) にする。さらに,出力書式もより広範囲の数値を表現できるよう %g を使う
正しい答えはいうまでもなく,x = -1, -1e-10 である。ここまでの改変後のプログラムでは
2次方程式の定数を入力してください。
a = 1
b = 1.0000000000000001
c = 0.0000000000000001
2次方程式の解: x = -1.11022e-16, -1
が得られる。絶対値の小さい方の解の誤差がひどい。原因は先に述べたように,2つの解を解の公式をそのまま使って求めているからである。大きい方の解をまず求めて,絶対値の小さい方の解は解と係数の関係から求めるようにしなければならない。
    if (b < 0) {
      kai1 = ( -b + sqrt(discriminant) ) / (2*a);
    }
    else {
      kai1 = ( -b - sqrt(discriminant) ) / (2*a);
    }
    kai2 = (c/a)/kai1;
これでやっと,正しい答えが得られる。
2次方程式の定数を入力してください。
a = 1
b = 1.0000000000000001
c = 0.0000000000000001
2次方程式の解: x = -1, -1e-16

このプログラムではちゃんと虚数解も求めるように作られている。

2次方程式の定数を入力してください。
a = 1    
b = 1
c = 1
2次方程式の解: -0.5+0.866025i, -0.5-0.866025i
前にも書いたが,解は printf で書かれている。実部と虚部を文字列としてつなぎ合わせているだけである
printf("2次方程式の解: %.2f+%.2fi, %.2f-%.2fi\n", real, imag, real, imag);
普通,二次方程式の解は別の計算のために使われるのだから,文字列でコンソールに表示するだけでは意義は半減する。C は複素数も扱えるのだから,複素数にすればよいのだ。もっとも,その場合に結果を出力するときは,creal, cimag で実部と虚部を取り出して出力しなければならないが。
#include
:
  double _Complex x1, x2;
:
    x1 = real+I*imag;
    x2 = conj(x1);

●● https://na-inet.jp/research/equation_cpp.pdf

は,ここに挙げたものの中では,一番優れていると思う。
C++ で,高校生向けの授業用資料であるが,「「桁落ち」現象と呼びます。これを解決する手段はいくつかありますが,ここでは最後の節で,倍精度浮動小数点数よりずっと多い桁数で計算する「多倍長浮動小数点数」というものを使って強引に解決していきます」と勇ましい(^_^;)。そこまで行かなくても,解の公式と解と係数の関係を使えばよいとは思う
また,虚数解にも対処している。解を虚数型変数で求めている。
complex<double> a, b, c;
complex<double> d, sol[2];
cout << "Solutions: " <
虚数解も,
Solutions: (1,1.41421), (1,-1.41421)
のように自然に表示される。

●● https://python-programming-iot.com/2018/08/26/61/

は,Python によるもの。sqrt(x) を x**0.5 とするなど,のけぞるが,Python では,x が負の場合には(自動的に)虚数になる。つまり,つねに虚数型を使って解の公式により解を求める(これはユニーク)。x が正の場合には x**0.5 は実数であるので,実数型の解が求まるが,桁落ちを考えていないので小さい方の解の精度が低くなることがある。

●● https://books.google.co.jp/books?id=bpYhLwff3zgC&pg=PA58&lpg=PA58&dq=%E4%BA%8C%E6%AC%A1%E6%96%B9%E7%A8%8B%E5%BC%8F+%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%A0&source=bl&ots=8B4TIiQMiT&sig=ACfU3U3i-YLWrRxyy1T0PtBi6wZnosrGqA&hl=ja&sa=X&ved=2ahUKEwjTpdiRtuPjAhWtUN4KHT5nC3w4PBDoATASegQIBxAB
書名: Javaで学ぶ数値解析,著者: 和光システム研究所

は,Java によるもの。double を使う,解の公式と解と係数の関係を使う。虚数解は実部と虚部を表示する。

●● https://blog.masuqat.net/2014/04/28/cancellation-of-significant-digits-of-quadratic-formula-with-csharp/

は,Java によるもの。float を使う。解の公式のみを使うのでは桁落ちがあるので,解と係数の関係を使うべしとの指摘をしている。虚数解の場合には,実部と虚部を表示する。

●● http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/letsc/1-body.html

は,double を使う,解の公式と解と係数の関係を使う。虚数解の場合は「実数解はありません」と表示する。

●● http://www.ritsumei.ac.jp/se/~osaka/rejime/suuti/equation.pdf

は,桁落ちの起きる原因(float を使う,解の公式をそのまま使う)を説明し,それを避ける方法(float のままであったとしても,解の公式と解と係数の関係を使う)を示している。

●● http://www.crl.nitech.ac.jp/~ida/education/MaterialsDesign/1453.html

は,プログラムの提示はなく,アルゴリズムを示すもの。解の公式のみを使うのではなく,解と係数の公式を使うべきであると述べている

●● http://www.ne.phen.mie-u.ac.jp/misc/equation.pdf

では,b^2 - 4ac が「正確に0であるかどうか」に注意するようにアドバイスをしている。
> d = b^2 - 4acが0かどうかをを正しく判断するためには,fabs(d/(b*b+fabs(4.0*a*c)))

しかし,2つの解を共に解の公式をそのままあてはめたり,虚数解を扱うとき複素数型を使わないという欠点は最初のプログラムと同じである。変数としては double を使い,出力にも %lg を使うことを勧めているのは当たり前のことではあるが,よい。

●● https://takuyab.com/archives/2613

は,「プログラムに解かせる為にはおなじみ“解の公式”を使います」と堂々と書いてある,float を使っている,虚数解は出力として形式的表示。

●● https://qiita.com/watarumohawk/items/c6fa8dfc9e2f2a23a53b

は,変数の型は double であるが,2つの解を共に解の公式をそのままあてはめたり,虚数解を扱うとき複素数型を使わないという欠点は最初のプログラムと同じである。

●● https://seaotter.cite.tohoku.ac.jp/coda/clang/c-3-branching.html

は,float を使い,2つの解を共に解の公式をそのままあてはめている。虚数解は考慮されていない(「実数解はない」と表示するのみ)。

●● https://algorithm.joho.info/programming/c-language/equation-solution-c/

は,関数形式のプログラムであるが,変数の型は double であるが,2つの解を共に解の公式をそのままあてはめたり,虚数解を扱うとき複素数型を使わない。

●● https://keisan.casio.jp/exec/system/1161228770

は,さすが CASIO だけあって,a = 1,b = 1.0000000000000001,c = 0.0000000000000001 のような場合であっても,誤差のない計算結果を返します。

●● http://pc-physics.com/secondequation1.html

は,「特にコンピューターでは虚数が扱えないので、虚数解になる場合注意が必要です」という,不穏当な注意書きがある。変数の型は double であるが,2つの解を共に解の公式をそのままあてはめたり,虚数解を扱うとき複素数型を使わない。

●● https://www.kkaneko.jp/cc/program/equation.html

は,解の公式と解と係数の公式を使うプログラムを示しています。しかし,虚数解の場合は,
    printf("x1 = 虚数解\n");
    printf("x2 = 虚数解\n");
という,ビックリするような出力をする。

●● http://roko.eng.hokudai.ac.jp/studentadm/chiba_data/pro/jp1.pdf

は,C,C++,FORTRAN の場合を示している。しかし,基本的に float を使う,解の公式のみを使う,虚数解の場合は「実数解はない」と表示するのみ。

●● http://kentiku-kouzou.jp/keisan-2.html

は,前述の CASIO の計算サイトと同じく,計算サービスのみでソースプログラムは提示していない。しかし,
a = 1,b = 1.0000000000000001,c = 0.0000000000000001 の場合,-1.1102230246251565e-16,-0.9999999999999999 を返すので,2つの解とも解の公式で求めているのだろう。「建築学生が学ぶ構造力学」というページなので,ちょっと不安になる

●● https://okayan08.hatenadiary.org/entry/20081219/1229677932

でギョッとさせられるのは,dの符号変換を d *=-1 とやっているような所(^_^;)この人は -b というような書き方ができないと思っているらしい。いつも -1*b+d と書いている。単項演算の - の存在を知らない(こういう人は案外多い)。
解の公式のみを使う,虚数解は実部と虚部を表示する。

●● http://www.u.tsukuba.ac.jp/~watanabe.norio.fw/H24-2FSubject.pdf

は,C によるもの。関数型で,double を使い,解の公式をそのまま使う,虚数解は実部と虚部を返すという仕様。「解答例」と称するにはちょっと不十分。

●● http://www.hino.meisei-u.ac.jp/is/iga/hit-u/No9/C3.pdf

は,float を使う,解の公式をそのまま使う,虚数解の場合は「実数解はありません」とのみ表示。

●● http://masudahp.web.fc2.com/jfbas/kiso/mon44.html
●● http://masudahp.web.fc2.com/n88basic/kiso/mon33.html

は,珍しく BASIC。解の公式をそのまま使う,虚数解の場合は「解なし」とのみ表示。

●● https://riptutorial.com/ja/fortran/example/3115/二次方程式

は,珍しく FORTRAN。real(float) を使う,解の公式そのままを使う,虚数解の場合は "No real roots." とのみ表示。

●● http://nalab.mind.meiji.ac.jp/~mk/labo/studying-C/Programing-in-C.pdf

は,C によるもの。double を使うが,解の公式そのままを使う,虚数解は実部と虚部を表示する。

●● https://wwws.kobe-c.ac.jp/deguchi/c/prob/eq.html

は,C によるもの。float を使う,解の公式そのままを使う,虚数解は「解なし」と表示する。

●● https://www.is.kochi-u.ac.jp/kyoko/edu/ip2010/python.html

は Python。必然的に double を使う。しかし,解の公式そのままを使う,虚数解は「実数解はありません」と表示する。

●● http://blog.livedoor.jp/add20/archives/897767.html

は,C, Perl, JavaScript。基本的にアルゴリズムが同じなら,同じような結果が表示される。解の公式そのままを使う,虚数解は「実数解は存在しない」と表示する。

●● http://macroscope.world.coocan.jp/ja/edu/computer/sabun/nizi.html

は,BASIC, Pascal, FORTRAN。基本的にアルゴリズムが同じなら,同じような結果が表示される。改定版として,double を使い,解の公式と解と係数の公式を使う方法が書かれている。虚数解の場合には実部,虚部を書くだけで,虚数型は使っていない。

●● http://nlp.dse.ibaraki.ac.jp/~shinnou/zemi2009/php/php-saito-1106.pdf

は,珍しい PHP でのプログラム。しかし,解の公式そのままを使う,虚数解は「実数解を持たない」と表示する。

●● https://excelmath.atelierkobato.com/quadratic-macro/

は,VBA によるもの。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● http://stroll.hatenablog.com/entry/2015/08/16/212016

は,Java によるもの。ここでも,-1 * b + Math.sqrt(d) のように,単項演算子 "-" を知らない表記をしているので,あきれる。いろいろゴタゴタ飾り立て,119 行のプログラムになっている。double を使っているものの,解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● http://www.kodama-lab.com/seminar/javalang/summer/ex001.html

は,Java によるもの。解の公式そのままを使う,虚数解は「None」と表示する。

●● https://www.kitasato-u.ac.jp/sci/resea/buturi/hisenkei/sogo/poly4.pdf

は,Python によるもの。
import sys, math, cmath
I=cmath.sqrt(-1)
omega1=(-1+cmath.sqrt(-3))/2
などの記述がある。I はこのように定義せずとも 1j でよい。
s1*s1*s1 などが頻繁に出てくるが s1**3 などの方が良いだろう。
3次方程式で,間接的ではあるが虚数型を使って解の公式のみを使っている

●● http://newtral.blog.jp/archives/1010749637.html

は,PHP によるもの。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● https://codeday.me/jp/qa/20190425/696210.html

は珍しい Prolog によるもの。解の公式のみを使う。虚数解の場合は「空リスト」を返す。

●● http://vipprog.net/wiki/解答例/二次方程式の解/Excel%20VBA.html

は,VBA(Visual BASIC) によるもの。Double を使う。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● https://www.e-bridge.jp/eb/tcontents/plain_basic/ChapC21.html

は,Plain Basic によるもの(行番号がある)。あれは double だったっけ?解の公式のみを使う。虚数解の場合は「根がありません」と表示する。Plain Basic といい,「根」というだけあって,古い

●● http://fullbasic.web.fc2.com/510_two_degree.html

は,full Basic によるもの(行番号がある)。穴埋め問題。「センター試験数学」なのでやむを得ないが,解の公式のみを使う。虚数解の場合は「解なし」と表示する。

●● https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=78&ved=2ahUKEwjTpdiRtuPjAhWtUN4KHT5nC3w4PBAWMBF6BAgEEAI&url=http%3A%2F%2Fwww.ocw.titech.ac.jp%2Findex.php%3Fmodule%3DGeneral%26action%3DDownLoad%26file%3D201017211-90-0-50.pdf%26type%3Dcal%26JWC%3D201017211&usg=AOvVaw2tSEKr0ecdF7qKxKG59MlB

は,C プログラムによるもの。double を使い,解の公式のみを使う。虚数解は実部と虚部を表示する。

●● https://teenaka.at.webry.info/201701/article_5.html

は,Python によるもの。解の公式のみを使う。虚数解は実部と虚部を表示する。

●● http://www.akita-pu.ac.jp/system/elect/ins/kusakari/japanese/teaching/Old/Programming/2005/note/5/Slide03.html

は,Java によるもの。解の公式のみを使う。虚数解の場合は「実数解はない」と表示する。

●● http://www.asahi-net.or.jp/~YY8A-IMI/20040913/x68/c_1.htm

は,C によるもの。解の公式のみを使う。虚数解は実部と虚部を表示する。

●● http://www.avis.ne.jp/~ammonite/javatest1.htm

は,Java によるもの。解の公式のみを使う。虚数解の場合は「回答なし」。

 

コメント    この記事についてブログを書く
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« 紙は何回折れるか? | トップ | Mac でのファイル名(UTF-8-M... »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

ブログラミング」カテゴリの最新記事