理想的２次方程式の解法を考えよう

お兄ちゃんの宿題スレ
http://piza.2ch.net/test/read.cgi?bbs=tech&key=982853418
の８９からの流れがあまりに面白いので

究極の２次方程式の解を得るコードを考えてみましょう

言語はC++
入力出力は実数の範囲

>>1
> 入力出力は実数の範囲
虚数解はどーすんの？

また、ax^2 + bx + c = 0がa = 0, b = 0, c != 0の場合、
解なしでいいの？

よく数値計算の本の最初の方に書いてある話。
x=(-b±√(b^2-4ac))/2aのように解の公式そのまま使うと
係数によっては桁落ちが深刻になるので
x=2c/(-b±√(b^2-4ac))のように変形してつかうべしということ。
ここまで考えてください。

>>3
その算法も桁落ちを回避できていません(a*cが小さいとき、根の一方は
bからbに非常に近い値を引いて求めるため)。

q = -( b + sgn(b)√(b^2-4ac) )/2 を求めて
x1 = q/a x2 = c/q とするのが良いとのことです。
(参考文献:Numerical Recipes in C,pp.158,ISBN4-87408-560-1)

>>2
複素数解は入力も複素数を許す場合でなければ意味がないと思われます
また複素数の２次式は実用的にはあまり使わないでしょう

実用的には実数入力で、複解が要求されない場合が多いと思います
例えば点列を与えて空間的最小２乗法で直線を推定する時とか

だから仕様としては
　複数解を同時に得られるか解の範囲を与えて答えを出すもの

というのではどうでしょうか？

>>3
　正確には一旦求めた絶対値の大きい方の解を使って それがx2なら　x1=c/(a*x2)
　ではなかったか？

>>4 すまんかぶった

>>4さんの方法ですが、

精度の向上をやるなら　重解がある時は　ニュートン法を１回通すのはどうでしょう？

f(x)=a*x^2+b*x+c
f'(x)=2ax+b
x+f(x)/f'(x)= 　x+(a*x^2+b*x+c)/(2*a*x+b)

class　SecondCurve　{
double　a,b,c;　//y=a*x*x+b*x+c
double　getD(void){return　b*b-4*a*c;　};
double　newton2(double x){ return x+(a*x*x+b*x+c)/(2*a*x+b) ;}
double　getX(int f){
　　　　　　　double　d=getD();
　　　　　　　 if(d< 0.0) throw err();
　　　　　　　 if(d==0.0) return -b/a/2;
　　　　　　　　else return　newton2((-b+f*sqrt(d))/a/2);
　　　　　};
public:
double　getN(void){ double D=getD();
　　　　　　　　　　if(D>0) return 2;
　　　　　　　　　　 else if(D<0) return 1;
　　　　　　　　　　 else return 0
　　　　　;};

double　getX1(void){ return getX(1);};
double　getX2(void){ return getX(-1);};
SecondCurve(double　aa,double　bb,double　cc){　a=aa;b=bb;c=cc;};
};

使い方は getNで解の個数を得てから、getX1,getX2で

この場合 a=0のチェックは
　getX()内部でやるべきか？
　それともコンストラクタで？

そもそもクラスにしない方がいいのか

一応作ってみました。なげーよ。
要件は満たしているでしょ。
http://chiba.cool.ne.jp/summoner/quad.txt

>>10
動作チェックしてますか? なんか動かないような気がするなぁ。テストドライバも見せてください。
ロジックには目を通してません。

明らかに不適切と思われる点をざっと指摘すると、
・constメンバ関数にするべきところがなっていない。
　　LinearEquation::getA()など
・throwするインスタンスはnewで確保するべきでない。
　　throw new IndefiniteException() → throw IndefiniteException()
　　根拠は「More Effective C++」参照
・非組み込み型を引数で渡すときはconst参照にするべき。
　　complex<F> QuadraticEquation::getAnotherComplexSolution(complex<F> one)
　　complex<F> QuadraticEquation::getAnotherComplexSolution( const complex<F>& one)
　　(Fはテンプレート型なので、Fの値渡しもconst参照にするべきか?)

以下は議論が分かれるところかな
・templateにする必要があるのか?
　　誤差の見積もりなんかが難しくなると思う。
・例外でエラーを通知すべきか?
　　解が存在しない係数の組み合わせは、想定していない例外的な振る舞い(メモリ不足など)では
　　ないので、例外で表現すべきではないのでは。
・コメントを打とう。
・浮動小数値を == 比較しているが、誤差範囲を考慮しなくて良いのか?

微分積分むずいんじゃｺﾞﾙｧ！！

最近C++してないから、いろいろ忘れてました。
一応動いたんだけど、まともなテストはしていません。つかそんな暇ない。
cout << QuadraticEquation<>(1, 2, 1).getSolution(); とかその程度。

テンプレートは、C++は自作の数値型すらありうる世界なので、あったほうがイイかと。

解を返す関数で、解がないときってどうするのがイイのでしょうか？
個数も一緒に返す奴は、個数0でやっているけど、解を一個返す関数は？
例外を投げる以外の良いエラー通知法は、ないような。

誤差の見積もり・・・　ちょっと考えたけどね・・・

ニュートン法だけでは解けないかな？
そうでないと折角テンプレートにしてもsqrtが必要になるし。

片方が求まれば　片方は 解と係数の関係
x1*x2 = c / a
x1+x2 = -b / a
から求まるでしょ？

>>7-8 のニュートン法の式は符号が間違いでした

>>14
　sqrtもニュートン法を使えばいいんじゃないかなと

x <- x-f(x)/f'(x)= 　x-(x*x-D))/(2x) = x-(x-D/x)/2;

double NewtonSqrt(double D)
{ double x=D/2;
for(int i=0;i<8;i++) x = x-(x-D/x) / 2 ;
return x = x-(x-D/x) / 2;
}

>>14 15
うーん、ニュートン法は f'(x)≒0のときに発散してしまうので、あまり
万能視しないほうがいいと思うよ。15でもD=0で0除算するよね。

とりあえず、>>11とニュートン法を採用してみました。
http://chiba.cool.ne.jp/summoner/quad.txt

>>17
おお、すごいなあ

BCC32でエラーが出たので
return d < 0 ? complex<F>(b, r) : s;
は
return d < 0 ?(complex<F>(b, r)) : (complex<F>(s,0));
かな？

計算後の精度向上にニュートン法を使えばと思ったけど、sqrtをニュートン法で出し
た場合、最後の１ビットの改善（かどうか判らない微妙な差）しかなかったです

コメント入れた。
http://chiba.cool.ne.jp/summoner/quad.txt

>>19
＞無限ループ防止の為、収束計算は64回で打切り
＞ 2^64桁以上の精度の処理系の人はごめんなさい

ってあるけど、実験してみたらニュートン法の収束は１回ごとに有効桁が
倍になる急峻なもの指数部さえあっていれば後５回も実行すれば１００桁
以上の精度が出る。

ただ、初期値が離れてる＝有効桁が無い状態　では　１回の演算で１ビット近
づくだけ（指数部が１加減されるだけ）となる。　この期間が遅い　から、

template<class F　> F calcurateSquareroot(const F &x) {
　F f=1,last;
　if(f<x)
　 { F oldf=f; f=2;F ff=4;
　　while(ff<x) { oldf=f; f=ff;ff=f*f;}; //２乗で近づくので高速
　　 f=oldf+(x-f)*(f-oldf)/(ff-f); // oldf〜fの間で直線補間
　 }else {
　　　F oldf=f; f=f/2;F ff=f*f;
　　while(ff>x) { oldf=f; f=ff;ff=f*f;};
　　 f=oldf+(x-f)*(f-oldf)/(ff-f);
　 }
　for (int i = 0; i < 10; i++) f = (x / f + f) / 2;
　return f;
}

麻衣ちゃんも　お局さんもやってこないなあ・・・・

>>20
今ちょっと試しているんだけど、引数を1.0+e200でやってみると・・・
まだまだ、お互い改良の余地ありですよ・・・

あ、そうか、　オーバーフローだ

型が固定小数点の場合もあるし　不動小数点数でも・・

オーバフローすると掛算したら前より値が小さくなる事を利用するとかしないとダメだな
チョッと考え中

オーバーフローも心配だけど、この場合はそうでなくて収束しきらないのですわ。
私のバージョンだとfor文のループを数百回、20でも、数十回まわさないとダメで。
回数さえ増やせば、√1e+200 ＝ 1e+100 は計算できるのだけど。

不動かよ。

あ　浮動　・・　スマン

なるほど>>24

>>20でやっても、まだ誤差が大きいんですね。
>>20で１回だけやってる簡易ニュートン法を
繰り返して、挟み込む方式ががいいかも・・・後で実験します

frexp()相当を作れば、オーダーは簡単に揃えられますね。
つーわけで、ためしに作った。引数が2以上でないと動かない奴だけど。

template<class F>
F tmp(const F &x, F &n, int &exp) {
F nn = n * n, cn = n;
if (nn >= x) return x / n;
int cexp = exp;
exp *= 2;
F t = tmp(x, nn, exp);
F cnn = n * nn;
if (cnn < x) { n = cnn; exp += cexp; return t / cn; }
n = nn;
return t;
}

template<class F>
F frExp(const F &x, int &exp) {
F f = 2;
exp = 1;
f = tmp(x, f, exp);
exp++;
return f / 2;
}

あ、あとpowかexp2がいるね。

そっちが王道みたい >>27

>>26　は　簡易ニュートン法だと誤差が大きすぎると全然収束しないからダメ

2のベキを求める関数は

template<class F > F Power2(int n) {
F r=1;
　switch(n){
　case -2: return r/4;
　case -1: return r/2;
　case 0: return r;
　case 1: return r*2;
　case 2: return r*4;
　default:
　 return Power2(m)*Power2(n-m);
　 }
}

だけど、ハテ？　これどうやって使うんだろ？　帰り値じゃ区別出来ないよね？

あ、
int m=n/2;

が抜けてる

正負で２の対数の概算がえられるように工夫してみたが・・・
なんか格好悪いなあ

template<class F　> int AboutLog2(const F &x) {
if(x==0) return 0;
if(x==1) return 1;
　F　y=1;
int r=0;
　F　 oldy;
int oldr;
if(x>y){ //1以上である
　 y=y*2;r++;
　 if(x<y) return 0;
　 if(x==y) return 1;
　　F　 oldy;
　　int oldr;
　 do{
　　 oldy=y;
　　 oldr=r;
　　　y=y*y;
　　　r=r+r;
　　　if(y<oldy) break;//オーバフローチェック
　　 }while(x>y); //xはoldy以上である
　return　 AboutLog2(x/oldy)+oldr;
}else{
　 y/=2;r--;
　 if(x>y) return　0;
　 if(x==y) return -1;
　 do{
　　 oldy=y;
　　 oldr=r;
　　　y=y*y;
　　　r=r+r;
　　　if(y==0) break;
　　 }while(x<y);
　return　 AboutLog2(x/oldy)+oldr;
}
}

ああ、最初の２行は間違い
if(x==0) エラーを投げる
if(x==1) return 0;

２次方程式ぐらい自分で解け！！何のために習ったんだ。

荒い対数を求めるの、少しだけスマートになった

　ところで関数内クラス定義って　テンプレート関数だと出来ないのかな？


template<class F>
struct _AboutLog{
　 F x; int r;
　　 _AboutLog(F x0){x=x0;r=0;}; //コンストラクタ
　
　int plus(const int &step , const F & pw){
　　　int step2=step+step;
　　　 F pw2 =pw*pw;
　　　 if((pw2>pw) //オーバフローチェック
　　　　&&(pw2<x ) ) plus(step2,pw2); //再帰を使って対数２のベキを求める
　　　　　while(x>pw)　{ x=x/pw; r=r+step;} //あとはバイナリ処理で
　　　　return r;
　　　};

　int minus(const int &step , const F & pw){
　　　int step2=step+step;
　　　 F pw2 =pw*pw;
　　　 if((pw2<pw) //アンダーフローチェック
　　　　&&(pw2>x ) ) minus(step2,pw2);
　　　　　while(x<pw)　{ x=x/pw; r=r+step;}
　　　　cout <<"("<< step <<" "<<pw<<　" "<< r <<")";
　　　return r;
　　　};
};
template<class F　> int AboutLog2(const F &x) {
_AboutLog <F> w(x);
F y=1.0;
if(x<=0){/*例外*/};
if(x>y){ //1以上である
　　　return w.plus(1,y*2);
}else{
　　　return w.minus(-1,y/2);
}
}

>>34
　いや自分で解こうとしてるんだけど？

汎用にしようとすると、色々あって思ったより難しいよ

いや、プログラムで解くんじゃなくて自分で計算しろってことだ。
プログラムにまかせるまでもないと思うんだが。

で、結局、　関数テンプレートやめて、テンプレートクラスにして挿入した。

これで　平方根の本体は簡単になった


template<class F >
struct _Power2 {F Power2(int n) {
F r=1;
switch(n){
case -2: return r/4;
case -1: return r/2;
case 0: return r;
case 1: return r*2;
case 2: return r*4;
default:
int m=n/2;
　return　Power2(m)*Power2(n-m);
　}
}
};
// ニュートン法にて、平方根を求めて返すテンプレート関数。引数 x は正のこと。
template<class F　> F calcurateSquareroot(const F &x) {
if(x<=0) return 0;
　_Power2<F> pw2;
　F f=pw2.Power2( AboutLog2(x)/2);
　for (int i = 0; i <8; i++) {
　　f = (x / f + f) / 2;
　}
　return f;
}

桁がほぼ合ってる状態からだから、
実際は５回でもOKと思う。ただ汎用だから８回にした

固定回数　ループするのは勿体ないようだけど、
８回くらいならまあいいかなと思う

>>37
ごめん見落としていた・・・
ウーン。
テンプレートで汎用のプログラムにするって事は、
必要な精度の実数クラス（加減乗除）さえ作れば
数１００桁の精度でも計算出来るって事になる訳で、
これを　自分で計算は無理だろうと思うのだが？

だったら２次方程式みたいに解の公式があるやつを扱っても
アピールしないのでは？ｎ次方程式の複素数解をちゃんと
求めるテンプレートがあればありがたいけど（ベアストー法とか
なんとかよりいいやつ）。

やっぱ乱数で入れるというのがもっとも早い

但し運良く正解が入った場合だけ・・・

ｎ次方程式の複素数解ですか・・・
　２次方程式は運動方程式とか特殊な最小2乗解に
　タマに使うけど、ｎ次は要求された事がないのですが
　どういう用途があるのでしょう？

ちなみに複素数そのものは使うのですが、２次方程式の
解としての複素数はやはり要求された事はありませんで
した。　これも用途があるなら知りたいです。

非線型の電気回路の解析とか，光学関係の設計とかに使う。
あとは減衰振動を扱う時は複素数解がいるね。

それはn次じゃなくてｎ階の方程式（いわゆる微分方程式）では？

ああ、微分方程式をラプラス変換で解くって奴ですか？

議論があらぬ方向に行っていますが、すこし考えをまとめました。

frexp と ldexp 相当が存在すれば、√x は、

int n;
F y = frexp(x, &n);
if (n % 2 != 0) y *= 2;

√x == ldexp(√y, n / 2)

で、求められるわけで、この場合の √y の計算では、0.5 <= y < 2 が保証されているため、
ニュートン法の初期値は、y < 1 なら 1、y >= 1 なら y で問題がない。
しかもこの場合、初期値が √y 以上である事が保証されるため、
ニュートン法の各過程では、計算値は減少する事が保障される。したがって、

F f = y < 1 ? 1 : y, last;
do { last = f; f = (y / f + f) / 2; } while (last > f);

で、収束する。多分数回のループで収束。
永久に減少しつづけることはないはずから、回数は指定しなくていいはず。

>>47
　うん正統だと思う。

DSPなんかでsqrtを実装する時は
その方法で処理した後、
３次くらいの級数展開
√(1+x)=1+x/2-x^2/8+x^3/16
を使って、精度を7ビットくらい出してから　必要な回数
（１〜３回）ニュートン法を使う事が多いみたい
割算の速度に応じて級数展開の次数と　ニュートン法の回数を調整する

でも、汎用となると、実数表現には固定小数点もあるわけで
必ずしもfrexpが効率的でない事があるのが難点かも

汎用でも効率を求めるなら、多倍長計算でも、
一旦 doubleで概算を計算してからニュートン法
が一番いいでしょうね。

特にハードウエアでsqrtが実装されていればね。

doubleへの代入くらいは普通用意するでしょ？

>49
ひとつだけ問題があります。
多倍長の場合、いや long double でさえ、指数範囲は double より大きい（事が多い）。
つまり、double に変換した時点でオーバー・アンダーフローする恐れがある。

そろそろネタが尽きたかなあ

やってみて結構勉強になった。特に10さんありがとう。

ん〜、とりあえず今までの議論を反映しつつ、STL風にしつつ、書き直してます。
んで、次はn次方程式に拡張する予定。

高次にしてなにするんですか？