10次の多項式を求めるプログラムについて

質問

ARMコア(32bit)を使ったCPUで多項式を求める演算をさせています(C)。

Y軸0-4095、X軸0-4095の範囲にプロットされる点が（整数値）11点で、非線形の点の集りです。多項式近似を使ってこの式を求め、求まった式に対してXに0～4095の値を代入し、4096個のY値を求めます。

絶対条件として、0≦x≦4095の範囲では、接線の傾きは必ず0よりも大です。

曲線式:
y=a + bx + cx^2 + dx^3 + ex^4 + fx^5 + gx^6 + hx^7 + ix^8 + jx^9 + kx^10

算出方法は、11点の(x,y)をfloat型の配列変数を使って行列を作成し、ガウスの消去法→後退代入させて求めてます。

この中でa～kの係数は、Windows上のCで同じソースで処理させたところ、ほぼ同一の結果が出ます。

Xの値を0から4095までひとつづつ代入して求めていくY値に対して、思った精度で結果が出ません。X値が大きい所で妙な計算結果が出ます。

係数a～kを使い、エクセル上でYを計算した結果、期待通りとなります。

例）
X値 エクセル結果 CPUでの結果 期待値
2962 2525.434 2526 2525
2963 2526.543 2528 2527
2964 2527.652 2527 2527
※2963⇒2964で接線の傾きがマイナスになる　→　おかしい

べき乗の計算は、掛算で対処しています。

float型で10乗するような計算処理をすること自体、無謀なのでしょうか？　

識者の方の意見をお伺いしたく思います。

回答No.6

rinkun · ベストアンサー

No.1は数値計算で展開式を計算するときの標準的手法です。
これは演算量を削減するとともにオーダーの大きく異なる数の加減算による情報落ち(情報欠落)を防止するものでもあります。

質問の曲線式を使った場合、10次の項はx=1000のとき10^30オーダーになります。各係数のオーダーが同等だと仮定すると定数項～7次項くらいまでは情報落ちでないも同然です。
No.1の式ではオーダー差が大きくならないうちに加減算することで情報落ちを少なくしています。

実のところ係数の符号が揃っていたりするとNo.1の式でも情報落ちは免れないのですが、実際は正負の係数が適当に出てきて中間結果が(最終結果に対して)あまり大きくならないようにできてます。

参考URL：

http://ja.wikipedia.org/wiki/%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%95%B0

お礼

実際にやってみましたが、↑の補足で書いたことはまったく的違いでした。情報落ちに対し、今回の10次の計算をしようとしたら、結構複雑な処理をさせないと実現しないのですね。　和算での積み残しをちょっとずつ拾って、最後に加える操作がいるようです・・・

きっかけを見出したので、もう少し調べてみたいと思います。ありがとうございました。

補足

確かに、いろんな係数を求めて、結果を出したところ、まったく情報落ちがなくなるというわけでなく、4096データのうち数パーセント程度はでてくるようです。ちょっと試しにやってみたら、20個ほど情報落ちしたものがありました。

そこで、参考書にもかいていたのですが、式の各項bx , cx^2 , dx^3などを10個のfloat型変数に格納して、値の小さいものから順次足し合わせていくほうがもっと精度があがるのではないかと考えましたが、いかがでしょうか？

回答No.5

Trick--o--

0.123
0.000123
どちらも、有効桁は３桁です。
これらを加算すると
0.123123
になりますが、有効桁が3桁なので
0.123
になってしまいます。

これが、＃４の方の言う
「大きな数値と小さな数値を加減算しようとしたことによる情報落ち」
だと思いました。

http://ja.wikipedia.org/wiki/%E8%AA%A4%E5%B7%AE#.E8.A8.88.E7.AE.97.E8.AA.A4.E5.B7.AE.E3.81.AE.E7.A8.AE.E9.A1.9E

お礼

情報落ちでいろいろ身の回りの参考書をみてますと、たくさんのってました。　クリアに理解できました。ありがとうございました。

回答No.4

sakusaker7

#2です。
有効桁を勘違いされていませんか?

たとえば有効桁が7桁ということはたとえ10^100 を表すことができたとしても、
1.234567*10^100 のように意味のある数字が7桁しかないということです。
つまり、12345678 と 12345679 を区別することができません。

#1の方のおっしゃるように式を変形することで解決したのであれば、
わたしの最初考えていたオーバーフロー/アンダーフローではなく、
大きな数値と小さな数値を加減算しようとしたことによる情報落ちでは
ないでしょうか。

お礼

ご回答ありがとうございます。
言われてみれば、そうですね。　ぼやっと理解してましたが、ご説明でクリアになりました。

回答No.3

麻野 なぎ（@AsanoNagi）

float の上下限ですが、表現形式が、IEEE 754 で、単精度（32bit）であれば、10^38 程度（2^127 程度）の表現能力しかありません。
倍精度（64bit）であれば、10^300 を超えます（2^1023 程度）

お礼

ご回答ありがとうございます。

回答No.2

sakusaker7

> float型で10乗するような計算処理をすること自体、無謀なのでしょうか？　

float だと十進6～7桁しか有効桁がありませんので、オーバーフローもしくは
アンダーフローが起きているということはありませんか?

お礼

ご回答ありがとうございます。
有効桁7桁というのは承知していましたが、指数表示されるはずで、float型は実質的に100乗あたりまでの数値を扱えたと、認識しています。

回答No.1

rinkun

> y=a + bx + cx^2 + dx^3 + ex^4 + fx^5 + gx^6 + hx^7 + ix^8 + jx^9 + kx^10

曲線の式をこれで計算しているなら精度は望めません。次の式に変形して計算してください。

y=a+x*(b+x*(c+x*(d+x*(e+x*(f+x*(g+x*(h+x*(i+x*(j+x*k)))))))))

お礼

ありがとうございます！
上記式であれば、先に述べた奇妙な計算結果にならず、増加傾向の値が算出されました。

でもどうしてなんでしょうか・・・
こういうやり方でやるのだ、という方法論的なものでしょうか。　あるいはもっと理論的なことがあるんでしょうか。　もし、そこらへんをご存知であれば、ご説明いただきたく思います。