CLNライブラリの精度を上げたいのですが・・・
#include <cln/cln.h>
using namespace std;
using namespace cln;
int main()
{
float_format_t precision = float_format(60);
cl_F x = "0.1_60";
cl_F j[52];
cl_F a,scale;
cl_F Fk = "50.0_60";
int n = 10;
int k = 50;
j[51] = j[50] = "1.0_60";
while(x<=40.0){
while(k>0){
a = 2.0 * Fk + 1.0;
j[k-1] = (a/x) * j[k] - j[k+1];
k = k - 1;
Fk = Fk - 1.0;
}
scale = ((sin(x))/x)/j[0];
cout << j[n]*scale << endl;
x = x + 0.1;
}
}
以上自分が組んでみたプログラムなんですけど、
一応これでも計算結果をプロットするとベッセル関数の波の形になるので間違ってはないと思うのですが、
_60とやっても60桁表示されません、せいぜい14桁くらいしか表示されなくて困っています。
どのようにしたら60桁の計算結果を出すことが出来るでしょうか?
float_format_t precision = float_format(60);
cl_F x = "0.1234567890123456789012345678901234567890123456789";
cout << x << endl;
の結果はどうなるんでしょうか?
これで出力桁が期待ほど長くないなら、出力に問題あんじゃないかと推測できますけども。
0.1234567890123456789012345678901234567890123456789L0
と表示されます。
ので、x=
0.1234567890123456789012345678901234567890123456789
と同じ表示がされています。
多倍長クラスというのを初めて聞いたのですが、少し興味があるので想像で申し訳ないですけれど横槍を入れさせてもらいます。
次の2つの行で1.0と2.0の精度はdoubleになりますので、ここで精度が落ちているということはないでしょうか?
a = 2.0 * Fk + 1.0;
Fk = Fk - 1.0;
例えば
a = (cl_F)2 * Fk + (cl_F)1;
Fk = Fk - (cl_F)1;
のように整数型とキャストを使ってみるなど
俺も想像だけで発言。違ったらごめん。
sin()もdoubleだからそこで精度落ちてない?
あとFPUは80bitだから、それ以上はどうやっても精度上がらないのでは?
FPUは使わずにint(整数レジスタ)のソフトウェアエミュレーションとかでやってるなら別だけど。
> sin()もdoubleだからそこで精度落ちてない?
については、ネットで仕様を確認してる限り一般的な関数はサポートしてるみたいですよ。
しっかりとnamespaceも使ってるので問題はないのかと思います。
> FPUは使わずにint(整数レジスタ)のソフトウェアエミュレーションとかでやってるなら別だけど。
CLNというのがそういうもののようです。
なので少し興味があったんですけれど、ぽこぴんさんはもう見てないのかな・・・
なるほど、そうでしたか。やはり嘘書いててごめんなさい。
ちなみにMSDNで検索してもヒットしてないんでまたもや想像ですけど、多倍長が128bitだと仮数部が100bit強ぐらいでしょうから、10進の60桁にはやはり足りない気がしますね。256bitなら足りるのかな?(電卓でも桁が足りなかったので計算はしてませんが…)
あ、ひょっとして無限に桁数(bit数)は増やせるのかな?ソフトウェアエミュレーションならそれも可能だな。
そのライブラリのマニュアルってどこにあるんですかね?
検索してて…
どうもマルチの上に、別の場所で厳しい突込みをくらい、質問を削除して逃亡した様ですね^^;
ttp://search.yahoo.co.jp/search?p=CLN%A5%E9%A5%A4%A5%D6%A5%E9%A5%EA&fr=top_v2&tid=top_v2&ei=euc-jp&search.x=1
私が見つけたのはこちらのPDFの文章ですが、中ほどにCLNについてふれています。
初期化するときに有効桁数を任意で決められるようで、理論上では有効桁数を無限に増やせそうです。
ですが、他のサイトにπを求めるプログラムをCLNで試した結果を載せているページがありましたが、桁数を増やすと急激にパフォーマンスが落ちるようで何百何千桁もの計算を行うのは現実的ではなさそうです。
ちなみに内部では浮動少数の演算プロセッサも使用されているようです。
何にしても質問主が逃亡したというのは残念です。
# ところでリンク先が「ttp://〜」で始まっていますが、「http://〜」から書き出すのはタブーなのでしょうか?
コメントに反応。
> # ところでリンク先が「ttp://〜」で始まっていますが、「http://〜」から書き出すのはタブーなのでしょうか?
ttpなるスキームが正しいのでしょう。
ttp://とhttp://は別物ですから。
# まさかプログラマがそれを区別しないわけがない。
当然ながら,httpスキームのURIはタブーでも何でもありません。
「ttp://」に深い意味はないです。ただ最近広告対策にhを抜かないと受け付けてくれないBBSが多いため、つい癖で抜いただけです^^;
多数の投稿ありがとうございました。
ヤフーで質問もしたのですが、解答を得られず
コメントを頂いたのですが、ある程度解決し、そのソースを載せようと思ったのですが、
二つも同じスレをたてるのも、どうかと思い削除してしまい
色々批判を受けて両方削除しました。
決して逃亡したわけでも、ムカついたから削除した等ではないのでご理解ください。
そこでCLNのソースなんですが、
#include <cln/cln.h>
using namespace std;
using namespace cln;
int main()
{
float_format_t precision = float_format(100);
int n = 0; //n次の値を出すためのn
double xf = 10; //xをどこまで計算させるか
double xia = 0.1; //x start2 // ●
cl_F x = "0.1_60"; //x start // ●
double delta = 1.0E-2; //x plus // ■
cl_F xd = "1.0E-2_60"; //x plus // ■
double xi;
cl_F j[52],scale,sum; //k=50まで計算させる配列,scale
cl_F Fk = "50.0_60"; //kと同期させるFk
cl_F Fk2 = "50.0_60"; //Fkの初期化用
cl_F b = "1.0_60"; // 数値1
cl_F c = "2.0_60"; // 数値2
j[51] = j[50] = "1.0_60";
cout << "x" << '\t' << "y" << endl;
for(xi=xia;xi<=xf;xi=xi+delta){ //xは0.1→40.0まで
Fk = Fk2; // Fk=50で初期化している
for(int k=50;k>0;k--){ // kは50→0.0まで計算
j[k-1] = (((c * Fk + b)/x) * j[k]) - j[k+1]; // 漸化式
Fk = Fk - b; // Fk と k を同期させている
}
scale = ((sin(x))/x)/j[0]; // j0=sinx/x
sum = j[n]*scale;
cout << x << '\t' << sum << endl;
x = x + xd;
}
}
このプログラムをすっきりさせる方法はありますでしょうか?
Fkなど二つあるのは、CLNの数値は1.0など加算してしまうと
桁落ちが発生してしまうのを防ぐためにあります。
そこで質問なのですが、y→0 のxを求めたいのですが(xの零点)
どのようにしたらいいのでしょうか?
if(sum<xd && sum >0.0){cout << x << endl;}
の文を考えてみたのですが、いまいちきれいな形とはいえないと思いますので・・・
他に何か助言がありましたらお願いします。
戻ってきたよ^^;
>決して逃亡したわけでも、ムカついたから削除した等ではない
本人がそう言っても真実か嘘か分からないしね…
Q&A形式のBBSで質問自体を削除するなんて最もやってはいけないことだと思うけど(後から来る人に何の役にも立たない)、それをする人に答える気はなりなりました。まあ誰か他の人が答えてくれるでしょ。
もともと興味があった話題だからマルチでも調べて答えようと思っていただけに残念。
ツイート | ![]() |