整数演算で立方根を求めるには?

解決


山田 由美子  2005-03-16 02:36:03  No: 56671

整数演算だけで 3√-300 ⇒ -669432 の様に
立方根(n乗根でも可)の値を10万倍とか100万倍等の整数(int)で
返す方法を教えて下さい。

一応整数部(-6)はニュートン法で 出来たけど
そのあとの小数部(69432)の求め方が分かりません。
宜しくお願い致します。(C言語)


REE  2005-03-16 03:43:57  No: 56672

元の数字を10万の3乗倍したものに対して3乗根を求めれば、
元の数字の3乗根の10万倍が得られます。


tetrapod  2005-03-16 17:53:58  No: 56673

http://www.nikonet.or.jp/spring/tikuji/tikuji.htm
とか。


山田 由美子  2005-03-17 02:28:25  No: 56674

REEさん>
>元の数字を10万の3乗倍したものに対して3乗根を求めれば
ありがとうございます。ですが、10万の3乗倍の時点で
オーバーフローしてしまいます。使える数値は32bitの
符号付整数なので 先に大きな値を掛け算することは出来ないのです。

tetrapodさん>
ありがとうございます。
その方法を そのまま整数演算だけで行った場合、
整数部だけしか求められないのではないでしょうか? 1の次は0なので。
私の質問は、その方法で з√10 の 収束値2(整数部)を得た後、
8の余りの2を どういう風に計算して小数部の1544を出すか?です。

目的は整数演算の立方根の関数を 出来るだけ小さなサイズの
マシン語にしたいのです。多少であれば、計算回数が多くなっても
マシン語サイズが小さくなるのであれば構いません。

その様な算法(ニュートン法の改とか)、または 手計算(開立法等)で
ソースが短く書けるのであれば そのソースを教えて下さい。


tetrapod  2005-03-17 03:28:47  No: 56675

計算する以前の段階の話になるわけですが
2.1544 をどのように整数で表現するつもりでしょうか?
それが決まらんことにはなんとも言いようがない。


tetrapod  2005-03-17 03:59:50  No: 56676

ああ、結果だけ100万倍したい、ってことね。
単純に 64bit 演算すればいいような気がしますけど。

Pentium 等 CPU を対象にするのであれば浮動小数点数で演算するほうが
コードが短くなりますけど、それぢゃダメですか?
fyl2xp1 fdiv f2xm1 だけでほとんど計算終わりですけど...


tetrapod  2005-03-17 19:19:44  No: 56677

f2xm1 が癖のある命令なので、結構めんどくさいことが判明。要 fscale ですね。
double ftroot(int x) {
    __asm fld1;
    __asm fld1;
    __asm fld1;
    __asm fld1;
    __asm fadd;
    __asm fadd;
    __asm fdiv;
    __asm fild [x];
    __asm fyl2x;
    __asm fld st(0);
    __asm frndint;
    __asm fsubr st(1),st;
    __asm fxch;
    __asm fchs;
    __asm f2xm1;
    __asm fld1;
    __asm faddp st(1),st;
    __asm fscale;
    __asm fstp st(1);
}
が多分最短 (結果は double だけど) 。
[0..2147483647] でしか答えが得られないのでアレげですが。
符号セットの命令を前後に追加しておけば負数も計算できますが宿題ということで。

# ここは VC++ 掲示板なのでこれでよいのですよね。


山田 由美子  2005-03-17 20:13:17  No: 56678

tetrapodさん わざわざ ありがとうございます。
まだ解決しておりませんが、tetrapodさんの親切に大変感謝しております。

>単純に 64bit 演算すればいいような気がしますけど。
マシン語は 32bit符号付整数しか扱えないインタリプタの
プログラミング言語に 読み込ませ それを実行させ、
高速な演算をさせる目的のものです。

64bit、浮動小数、型変換、や 命令、関数(インラインのabsのみOK)は、
他言語でマシン語を実行させるための制限から使うことは出来ません。

使うことが出来るのは 32bit符号付整数と演算子です。
すべて自前で計算を行う必要があります。

># ここは VC++ 掲示板なのでこれでよいのですよね。
すみません。言語は関係なくアルゴリズムだけ知りたかったのです。
ここが数学とプログラミング両方に関係していると思ったので
よく分からずに こちらで質問してまいました。
場違いでしたら適した場を教えて頂いたら そちらに移動します。

この問題は、他所でも何度か質問した事がありますが、未だ誰にも
回答が得られていないものです。簡単そうで難問かもしれません。


山田 由美子  2005-03-17 20:26:42  No: 56679

自己レス まちがい訂正
>インタリプタ
インタプリタ(^^;

>関数(インラインのabsのみOK)は、
>他言語でマシン語を実行させるための制限から使うことは出来ません。
正確にはポインタでの呼び出しは可能。この場合意味がありませんが。


tetrapod  2005-03-17 21:44:15  No: 56680

いやだからアルゴリズムは提示済みです。あとは応用してください。
例えば32ビットの演算が出来れば64ビットの演算ができることは自明です。

三乗根の有効数字としてどの程度必要なのでしょうか?
小数点以下1桁あればいい程度の話であれば、説明するまでもなく簡単に整数演算できます。
わざわざ掲示板に書くほどの内容ではないと思ったのですけど。


ryou  2005-03-18 00:28:25  No: 56681

僕も考えてみましたが これは難しい問題かもしれませんね。
立方根でググってみましたが 山田さんの なさりたい事は見つけられませんでした。

僕が時々ROMってる「HSPラウンジ」というのがあるのですが
http://dream.freespace.jp/perl-bin/puma/lng/joyful.cgi
言語はHSPというものなんですが、無理難題OK!と書かかれてあるので
たぶん Cであっても 管理人さんが 教えてくれると思います。
そこで聞かれれば 解決率が99%以上?なので 適切な回答が得られるかもしれません。
ここでは理解すらしてもらえないと思いますので ダメ元で 
そちらで質問してみては?僕もどうやるのか知りたいので

>小数点以下1桁あればいい程度の話であれば、
10万倍とか100万倍等の整数(int)で返す方法を教えて下さい。という質問なので


tetrapod  2005-03-18 02:23:23  No: 56682

立方根の100万倍の値を得て何に使うのかがさっぱり判らんので答えようが無いんですよ。
後で演算する目的にも使えんし。
原理的に不可能なことはなにをどうやろうと不可能。
・元数値の有効数字が10桁しかないのにそれ以上の精度を求めることに意味が無い。
・途中の演算における有効数字に 32bit しか使えないのに、
  結果は 32bit 以上の有効数字を求めているし。

表示するだけ?

簡単な対策:
元数値を8倍して演算し、結果を 1/2 すれば2進数での小数点以下1桁が得られる。
元数値を64倍すれば結果を 1/4 すれば2進数小数点以下2桁が得られる。


REE  2005-03-18 03:11:24  No: 56683

・元数値の有効数字が10桁しかないのにそれ以上の精度を求めることに意味が無い。
・途中の演算における有効数字に 32bit しか使えないのに、
  結果は 32bit 以上の有効数字を求めているし。

求められているのは、小数点以下5〜6桁の精度ですよ
32bit符号付整数の最大値の立方根の100万倍は32bit符号付整数で表せますよ。

#  実際に計算できるかどうかは別問題ですが。


REE  2005-03-18 03:12:11  No: 56684

引用記号を入れ忘れました。  最初の3行が引用です。


ryou  2005-03-18 03:41:21  No: 56685

REEさんの言う通りなんですが、補足。
tetrapodさん VC++から頭を切り離して もう一度山田さんのを
最初から読み直せば 意味を理解出来ると思いますが。

フリーのプログラミング言語などで32bit整数(int)しか対応していないもの
いくつかありますが その場合 3√7*12345 を計算しようとしても
intでは 3√7は 1にしかならないので 3√7*12345=12345というとんでもない
誤差が出てしまう訳ですよ。それを 3√7をたとえば 1万倍した値が
求められると 3√7は19129になって 19129*12345/10000 で 2360という
そこそこの精度が得られ、ゲームなどには十分使える精度のものになります。
そういうことだと思います。


WIZ  2005-03-18 18:45:36  No: 56686

tetrapod さんへ。

> 原理的に不可能なことはなにをどうやろうと不可能。
> ・元数値の有効数字が10桁しかないのにそれ以上の精度を求めることに意味が無い。
> ・途中の演算における有効数字に 32bit しか使えないのに、
>  結果は 32bit 以上の有効数字を求めているし。

ryou さんも既に仰られていることなのですが、質問の意味を理解していないようです。
実際、コンピュータで円周率は何万桁まで計算されていますよね。でも計算に使われて
いるコンピュータは演算の有効数字が10進数で何万桁もある訳ではなく、精々数桁から
数十桁でしょう。

私も立方根(n乗根)を求める数学的アルゴリズムの存在有無については興味があります。
開立については http://yosshy.sansu.org/cbr.htm があり、C(C++) プログラムに変更
してみたのですが、上記ホームページ記載のアルゴリズムが間違っているのが、私の
コードが間違っているのか満足な結果は得られていません。


外野  2005-03-18 19:28:34  No: 56687

>実際、コンピュータで円周率は何万桁まで計算されていますよね。

これを整数演算だけでやってるとは思いません

>でも計算に使われているコンピュータは演算の有効数字が10進数で何万桁もある訳ではなく、>精々数桁から数十桁でしょう。

その言語でも動くように整数配列か文字列で BigInteger を作ればいい
後は何万桁の演算でも時間さえかければ思いのままです


REE  2005-03-18 19:54:53  No: 56688

開立法(立方根を筆算)は下記のページにもありますが、
http://www004.upp.so-net.ne.jp/s_honma/root.htm
第二副運算の桁数が、答えの2倍になっているので、32bitでは不足しそうです。


tetrapod  2005-03-18 20:13:28  No: 56689

判ってないといわれるのははなはだ心外なので:
でも手抜き。説明も無し。32bit 全ての範囲で結果が得られもしません。精度も無し。
int qroot(int x) {
    // An+1=(2/3)An + (1/3)(A/An^2)
    int y;
    int z;
    int n=0;
    if (x>10000) z=100*100; else z=1*100;
    do {
        y=z;
        z=(y*2/3) + ((x*100/y)*100/y)*100/3;
//      std::cout << "x=" << x << " y=" << y << " z=" << z << std::endl;
    } while (++n<20);
    return z;
}


ryou  2005-03-18 22:59:00  No: 56690

>判ってないといわれるのははなはだ心外なので:

>   if (x>10000) z=100*100; else z=1*100;
だからぁ。そういうのは最初にREEさんが答えて
そういう問題では無いと言ってるでしょ。
平方根を10倍、100倍のレベルであれば 中学生でも思いつきますよ
これは立方根で100万倍が出来ないかという問題。

intだけの立方根計算で最初に1000倍や1000000倍したら
いくつまで求められますか?返値が12までしか求められません。
では役に立ちませんよ。


REE  2005-03-19 00:18:17  No: 56691

ryouさんご紹介の「HSPラウンジ」で解決したらしいです。
山田さんには具体的な解決方法を報告して欲しいですね


WIZ  2005-03-19 04:06:03  No: 56692

外野さんへ。
# 既に解決しているらしいので、ゴミレスですが。

> これを整数演算だけでやってるとは思いません
> その言語でも動くように整数配列か文字列で BigInteger を作ればいい
> 後は何万桁の演算でも時間さえかければ思いのままです

貴方も質問の意味や、私の書き込みの意味を理解できていないようです

私の書き込みでは整数かどうかは別として、tetrapod さんの意見に対する
反論です。要約すれば、高い精度の計算を有限回の低い精度の計算に分解
するアルゴリズムの例として円周率の計算があるということです。

> その言語でも動くように整数配列か文字列で BigInteger を作ればいい
> 後は何万桁の演算でも時間さえかければ思いのままです

では、C 言語で実現できるアルゴリズムを示してください。
下位の桁を保持する整数がオーバーフローして上位の桁を保持する整数へ
波及していく部分をどう実装すればよいでしょうか?
機械語でキャリーフラグを参照してなんていう方法はダメですよ。


外野  2005-03-19 04:57:59  No: 56693

>私の書き込みの意味を理解できていないようです

>高い精度の計算を有限回の低い精度の計算に分解
>するアルゴリズムの例として円周率の計算があるということです。

読解力無くて御免。最初の書き込みからこれには直接つながらない。

> 下位の桁を保持する整数がオーバーフローして上位の桁を保持する整数へ
> 波及していく部分をどう実装すればよいでしょうか?

ryou 氏の例示なら早く近似値が取れることが大事だろうと思うけど、
円周率の例だと WIZ 氏の主眼は有効桁を超える値の結果がだせる方にない?

制限付のインタプリンタ用なんだしメモリも時間も効率なんてあんま考慮せずに、
1桁1整数で確保して10 になったら隣に繰り上げていけばいいだけのはず。
C みたいにビット演算できる言語ならキャリー用に1ビット捨てるだけでもいいけど。


山田 由美子  2005-03-20 03:43:02  No: 56694

ryouさん ご紹介の「HSPラウンジ」で 失礼ながら 同じ質問をさして頂きました。
そこでの ぷまさん回答の平方根の例では 立方根では うまく出来ません
でしたので メールで コソっと伺ってみたところ いろいろな整数での
計算テクニックを 整数演算だけで360度三角関数(πの計算)を
10億倍の精度で求める方法を例に いろいろ教えて頂き大変感動致しました。

どの様に解決したかは 難しくて まだあまり解っていないのですが、
要は 私の質問の場合、1万桁を求める訳ではないので、大きな数を延々と
掛け続ける必要は全く無い(不要な精度)とか、大きな数の立方根を求める場合は、
最初に1000等で割ってもいいということでした。 その他いろいろ。
説明が3行で終わってしまった。(^^;

回答して頂いた皆様 有難うございました。


※返信する前に利用規約をご確認ください。

※Google reCAPTCHA認証からCloudflare Turnstile認証へ変更しました。






  このエントリーをはてなブックマークに追加