理想的２次方程式の解法を考えよう
■ このスレッドは過去ログ倉庫に格納されています
0001デフォルトの名無しさん
NGNGhttp://piza.2ch.net/test/read.cgi?bbs=tech&key=982853418
の８９からの流れがあまりに面白いので
究極の２次方程式の解を得るコードを考えてみましょう
言語はC++
入力出力は実数の範囲
0003デフォルトの名無しさん
NGNGx=(-b±√(b^2-4ac))/2aのように解の公式そのまま使うと
係数によっては桁落ちが深刻になるので
x=2c/(-b±√(b^2-4ac))のように変形してつかうべしということ。
ここまで考えてください。
0004デフォルトの名無しさん
NGNGその算法も桁落ちを回避できていません(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)
0005１
NGNG複素数解は入力も複素数を許す場合でなければ意味がないと思われます
また複素数の２次式は実用的にはあまり使わないでしょう
実用的には実数入力で、複解が要求されない場合が多いと思います
例えば点列を与えて空間的最小２乗法で直線を推定する時とか
だから仕様としては
複数解を同時に得られるか解の範囲を与えて答えを出すもの
というのではどうでしょうか？
>>3
正確には一旦求めた絶対値の大きい方の解を使って それがx2なら x1=c/(a*x2)
ではなかったか？
0006１
NGNG0007１
NGNG精度の向上をやるなら 重解がある時は ニュートン法を１回通すのはどうでしょう？
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)
0008１
NGNGdouble 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で
0009１
NGNGgetX()内部でやるべきか？
それともコンストラクタで？
そもそもクラスにしない方がいいのか
0010例のスレの122
NGNG要件は満たしているでしょ。
http://chiba.cool.ne.jp/summoner/quad.txt
0011デフォルトの名無しさん
NGNG動作チェックしてますか? なんか動かないような気がするなぁ。テストドライバも見せてください。
ロジックには目を通してません。
明らかに不適切と思われる点をざっと指摘すると、
・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にする必要があるのか?
誤差の見積もりなんかが難しくなると思う。
・例外でエラーを通知すべきか?
解が存在しない係数の組み合わせは、想定していない例外的な振る舞い(メモリ不足など)では
ないので、例外で表現すべきではないのでは。
・コメントを打とう。
・浮動小数値を == 比較しているが、誤差範囲を考慮しなくて良いのか?
0012名無しのエリー
NGNG0013>11
NGNG一応動いたんだけど、まともなテストはしていません。つかそんな暇ない。
cout << QuadraticEquation<>(1, 2, 1).getSolution(); とかその程度。
テンプレートは、C++は自作の数値型すらありうる世界なので、あったほうがイイかと。
解を返す関数で、解がないときってどうするのがイイのでしょうか？
個数も一緒に返す奴は、個数0でやっているけど、解を一個返す関数は？
例外を投げる以外の良いエラー通知法は、ないような。
誤差の見積もり・・・ ちょっと考えたけどね・・・
0014デフォルトの名無しさん
NGNGそうでないと折角テンプレートにしてもsqrtが必要になるし。
片方が求まれば 片方は 解と係数の関係
x1*x2 = c / a
x1+x2 = -b / a
から求まるでしょ？
0015１
NGNG>>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;
}
0016デフォルトの名無しさん
NGNGうーん、ニュートン法は f'(x)≒0のときに発散してしまうので、あまり
万能視しないほうがいいと思うよ。15でもD=0で0除算するよね。
001710
NGNGhttp://chiba.cool.ne.jp/summoner/quad.txt
0018１
NGNGおお、すごいなあ
BCC32でエラーが出たので
return d < 0 ? complex<F>(b, r) : s;
は
return d < 0 ?(complex<F>(b, r)) : (complex<F>(s,0));
かな？
計算後の精度向上にニュートン法を使えばと思ったけど、sqrtをニュートン法で出し
た場合、最後の１ビットの改善（かどうか判らない微妙な差）しかなかったです
001910
NGNGhttp://chiba.cool.ne.jp/summoner/quad.txt
0020１
NGNG＞無限ループ防止の為、収束計算は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;
}
0021１
NGNG002210
NGNG今ちょっと試しているんだけど、引数を1.0+e200でやってみると・・・
まだまだ、お互い改良の余地ありですよ・・・
0023１
NGNG型が固定小数点の場合もあるし 不動小数点数でも・・
オーバフローすると掛算したら前より値が小さくなる事を利用するとかしないとダメだな
チョッと考え中
002410
NGNG私のバージョンだとfor文のループを数百回、20でも、数十回まわさないとダメで。
回数さえ増やせば、√1e+200 ＝ 1e+100 は計算できるのだけど。
0025デフォルトの名無しさん
NGNG0026１
NGNGなるほど>>24
>>20でやっても、まだ誤差が大きいんですね。
>>20で１回だけやってる簡易ニュートン法を
繰り返して、挟み込む方式ががいいかも・・・後で実験します
002710
NGNGつーわけで、ためしに作った。引数が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;
}
002810
NGNG0030１
NGNGtemplate<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);
}
}
だけど、ハテ？ これどうやって使うんだろ？ 帰り値じゃ区別出来ないよね？
0031１
NGNGint m=n/2;
が抜けてる
0032１
NGNGなんか格好悪いなあ
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;
}
}
0033１
NGNGif(x==0) エラーを投げる
if(x==1) return 0;
0034デフォルトの名無しさん
NGNG0035１
NGNGところで関数内クラス定義って テンプレート関数だと出来ないのかな？
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);
}
}
0036１
NGNGいや自分で解こうとしてるんだけど？
汎用にしようとすると、色々あって思ったより難しいよ
003734
NGNGプログラムにまかせるまでもないと思うんだが。
0038１
NGNGこれで 平方根の本体は簡単になった
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;
}
0039１
NGNG桁がほぼ合ってる状態からだから、
実際は５回でもOKと思う。ただ汎用だから８回にした
固定回数 ループするのは勿体ないようだけど、
８回くらいならまあいいかなと思う
0040１
NGNGごめん見落としていた・・・
ウーン。
テンプレートで汎用のプログラムにするって事は、
必要な精度の実数クラス（加減乗除）さえ作れば
数１００桁の精度でも計算出来るって事になる訳で、
これを 自分で計算は無理だろうと思うのだが？
0041デフォルトの名無しさん
NGNGアピールしないのでは？ｎ次方程式の複素数解をちゃんと
求めるテンプレートがあればありがたいけど（ベアストー法とか
なんとかよりいいやつ）。
0042デフォルトの名無しさん
NGNG但し運良く正解が入った場合だけ・・・
0043１
NGNG２次方程式は運動方程式とか特殊な最小2乗解に
タマに使うけど、ｎ次は要求された事がないのですが
どういう用途があるのでしょう？
ちなみに複素数そのものは使うのですが、２次方程式の
解としての複素数はやはり要求された事はありませんで
した。 これも用途があるなら知りたいです。
0044デフォルトの名無しさん
NGNGあとは減衰振動を扱う時は複素数解がいるね。
0045１
NGNG0046１
NGNG004710
NGNGfrexp と 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);
で、収束する。多分数回のループで収束。
永久に減少しつづけることはないはずから、回数は指定しなくていいはず。
0048１
NGNGうん正統だと思う。
DSPなんかでsqrtを実装する時は
その方法で処理した後、
３次くらいの級数展開
√(1+x)=1+x/2-x^2/8+x^3/16
を使って、精度を7ビットくらい出してから 必要な回数
（１〜３回）ニュートン法を使う事が多いみたい
割算の速度に応じて級数展開の次数と ニュートン法の回数を調整する
でも、汎用となると、実数表現には固定小数点もあるわけで
必ずしもfrexpが効率的でない事があるのが難点かも
0049デフォルトの名無しさん
NGNG一旦 doubleで概算を計算してからニュートン法
が一番いいでしょうね。
特にハードウエアでsqrtが実装されていればね。
doubleへの代入くらいは普通用意するでしょ？
005010
NGNGひとつだけ問題があります。
多倍長の場合、いや long double でさえ、指数範囲は double より大きい（事が多い）。
つまり、double に変換した時点でオーバー・アンダーフローする恐れがある。
0051１
NGNGやってみて結構勉強になった。特に10さんありがとう。
005210
NGNGんで、次はn次方程式に拡張する予定。
■ このスレッドは過去ログ倉庫に格納されています