C言語でべき乗表現

こんにちは、NAKAといいます。

C言語基礎を確認させてください。

pow()っていう関数を使うと思うのですが

例えば2^15って0x8000だと思うのですが、実際やってみると0x7FFFになってしまいます。

Parents Reply
  • NoMaYさん、NAKAです。

    数学的本質みたいなのは理解できてませんが、教えていただいたページは興味深く見させて頂きました。
    0.5乗はルートだし、そもそも^はEXORですものね!

    VBでは^がべき乗みたいですが、今更なにを!って感じですかね!

    ちょっと前まで、switch~caseをselect~caseって打って怒られてました (*_*;

    諸悪の原因はEXCELかな?


Children
  • NAKAさん、こんにちは。NoMaYです。

    > 諸悪の原因はEXCELかな?

    たぶん、NAKAさんの場合、まずはそこでしょうかね。ちなみに、ググってみたのですが、EXCELのべき乗演算(x^y)は、xにもyにも浮動小数点数が使えるようですね。

    [Excel] べき乗・累乗の計算 - Excel VBA 数学教室
    excelmath.atelierkobato.com/power/

    [関連リンク]

    浮動小数点数の内部表現に関しては以下の検索結果を時間のある時に見ていくと知識が溜まるかな、と思います。あとは、xもyも浮動小数点数である場合の、べき乗の数値計算ライブラリの計算精度、の話のウェブページがあれば良いのかなぁ。(その前に、xもyも浮動小数点数である場合の、べき乗の数値計算ライブラリの計算方法、の話の別のウェブページかなぁ。)

    Google検索: IEEE浮動小数点数 site:qiita.com
    www.google.com/search?q=IEEE浮動小数点数+site:qiita.com

    float型では123456789すらも表現できない話
    qiita.com/y-yoshinari/items/76260f6359d5b4418b33
     

  • NAKAさん、こんにちは。NoMaYです。

    > 浮動小数点自体が仮数部で表現できるのに限界があるため、誤差が発生し

    むしろ、べき乗の計算方法がNAKAさんが思っているのとは全く別物、という点に思い至っていない、のかなと私は思うのです。


    > CC-RLライブラリでの計算方法はともかくとして、浮動小数点計算をするものである、単純に x を y 回だけ掛け算するものでは無い、というイメージとしてです。

    > その前に、xもyも浮動小数点数である場合の、べき乗の数値計算ライブラリの計算方法、の話の別のウェブページかなぁ。

    SAMマイコンでも、べき乗を計算する数値計算ライブラリは、コンパイラがGCCとかClangとかなら、newlibと呼ばれるオープンソースライブラリが使われているのではないかな、と、思いますけれども、そのソースコードは以下のようなものになります。繰り返しになりますけれども、単純に x を y 回だけ掛け算するものでは無い、のです。

    github.com/eblot/newlib/blob/master/newlib/libm/math/e_pow.c#L16-L23

         *                    n
         * Method:  Let x =  2   * (1+f)
         *  1. Compute and return log2(x) in two pieces:
         *          log2(x) = w1 + w2,
         *     where w1 has 53-24 = 29 bit trailing zeros.
         *  2. Perform y*log2(x) = n+y' by simulating multi-precision
         *     arithmetic, where |y'|<=0.5.
         *  3. Return x**y = 2**n*exp(y'*log2)

     

  • NoMaYさん、こんにちは。Sugachanceです。

    たとえば序盤に出てきた2^15 = 32768が
    floatで受けると 3.276798E+004
    になっていた事象ですが、
    これはNAKAさんの仰るような「 浮動小数点自体が仮数部で表現できるのに限界が~
    によるものではなく、計算アルゴリズムの時点で既に誤差が発生しているという見解なのでしょうか?

    doubleで受けても確かに3.276800000000004E+004 なので、
    計算アルゴリズムの時点で既に誤差というのは何となく理解できますが、

    計算アルゴリズムの時点で既に誤差 -> floatに入れることによりさらに表現限界の誤差 -> intでさらに…

    というのが最初の事象として正しいような気がするのですがいかがでしょうか?
    (関数内部の計算アルゴリズムの時点で登場する変数が
    正確に全てを表現できる理想的な変数で、それが返ってくるという仕様のものであれば、
    double,floatの表現限界って話にまとまるのでしょうけども...)

  • Sugachanceさん、こんにちは。NoMaYです。

    >  floatに入れることによりさらに表現限界の誤差

    以下のQiitaの投稿で知ったのですけれども、ウェブ上で浮動小数点数の内部表現を確認出来るサイトがありました。そのサイトで、32768(とか32767とか32769とか)の単精度浮動小数点数の内部表現を確認すると、ちゃんと数値全体が収まっていますよ。なので、これらの数値ぴったりの数値が単精度浮動小数点数の内部表現限界の外にあるわけではないのですよ。

    浮動小数点数内部表現シミュレータ
    tools.m-bsys.com/calculators/ieee754.php


    0.1の浮動小数点数は0.1より大きいのに、10回足すと1.0より小さいのはなぜか【前編】
    qiita.com/yucatio/items/85a8e8d381d300092295

    0.1の浮動小数点数は0.1より大きいのに、10回足すと1.0より小さいのはなぜか【後編】
    qiita.com/yucatio/items/b5ed14bbd1e9a751b703
     

  • NAKAです。

    NAKA:浮動小数点自体が仮数部で表現できるのに限界が~


    Qiita:0.1を浮動小数点数にすると0.1よりわずかに大きくなるのに、10回足すと1.0より小さいのはなぜか、の答えは、2つの浮動小数点数を足す際に、仮数に入りきらない数を丸めるから、でした。

    同じ事のような気が........

    いっそ今回のような場合は自分で累乗の関数を作ったほうが早かったですね!

  • NAKAさん、こんにちは。NoMaYです。

    > の答えは、 2つの浮動小数点数を足す際に、仮数に入りきらない数を丸めるから、でした。 同じ事のような気が........

    いやいや、今回は、以下なのですよ。32768は、数値自体は、内部表現的には、丸める必要など無い数値、なのですよ。

    > 以下のQiitaの投稿で知ったのですけれども、ウェブ上で浮動小数点数の内部表現を確認出来るサイトがありました。そのサイトで、32768(とか32767とか32769とか)の単精度浮動小数点数の内部表現を確認すると、ちゃんと数値全体が収まっていますよ。なので、これらの数値ぴったりの数値が単精度浮動小数点数の内部表現限界の外にあるわけではないのですよ。

    そして、丸める必要が無い、という点では、以下の数値もそうである、なのですよ。

    2、4、8、16、32、64、128、256、、、16384

  • NoMaYさん NAKAです。

    >今回は、以下なのですよ。32768は、数値自体は、内部表現的には、丸める必要など無い数値、なのですよ

    ⇒確かに!

    単純なNAKAは上記のような関数を作ったほうが早かったなあ!
    そもそも、依頼元の意図が汲み取れれば2^15を0x8000に変更すればよかっただけだし
    でも、自信が無いから、大本に出来るだけ忠実にしたほうが....と考えちゃう!(;_;)

    でも、ちょっと盛り上がりましたね!(^_-)-☆

  • NoMaYさん、こんにちは。Sugachanceです。

    解説ありがとうございます。(そのシミュレータはたまに利用させてもらってます!便利ですよね)

    32768そのものをfloatに入れる分に関してはおっしゃる通り問題ないと思うのですが
    CのPow()関数がこれまでに出てきたアルゴリズムで実装されているとして
    ```
    pow関数の中身のアルゴリズムの時点で(^の表現を使ってしまいますが)
    pow(2,15) = e^(y*log(x))  = 32768とはちょっとずれた値
    になってしまっている。
    これをfloatに渡すと近い値の3.276798E+004
    doubleに渡すと近い値の3.276800000000004E+004

    ```
    という形"ではない"という事ですかね?

    Windowsの電卓でたたいてみましたが(おそらくdecimal型?)
    ln(2) = 0.69314718055994530941723212145818‬
    15 * ln(2) = 10.397207708399179641258481821873
    e^(15 * ln(2) ) = 32768
    ときれいに計算できました。


    一方で内部表現シミュレータに15*ln(2)の結果を入れると
    10.397207708399179641258481821873
    10.397207708399179641258481821
    どちらもdouble型では0x4024cb5ecf0a9651
    となりました。
    Windowsの電卓では
    e^10.397207708399179641258481821 = 32,767.999999999999999999999971405‬
    となったので、コンパイラ依存なのか計算機依存なのかは分かりませんが
    こういう内部の差が出たのではないかと考えていたのですが…


  • >でも、ちょっと盛り上がりましたね!

    皆様の回答を含めて、大変勉強になります。
    興味深い話題提供、ありがとうございます!

  • Sugachanceさん、こんにちは。NoMaYです。

    > これをfloatに渡すと近い値の3.276798E+004
    > doubleに渡すと近い値の3.276800000000004E+004

    これらの「渡す」という言い回しが私にはピンと来なかったのですけれども、私なら「計算する」という言い回しを使うのですけれども。

    単精度計算コンパイルオプション版 pow(2.0, 15.0) を、そのコンパイラの単精度浮動小数点ライブラリで、計算したら戻り値は○○○
    ↑ doubleが32bitで、floatも32bit、の意

    倍精度計算コンパイルオプション版 pow(2.0, 15.0) を、そのコンパイラの倍精度浮動小数点ライブラリで、計算したら戻り値は□□□
    ↑ doubleが64bitで、floatは32bit、の意

     
    他方で、単精度計算版と倍精度計算版で、powf()とpow()というように関数が異なって、かつ、floatが32bitで、doubleが64bit、という場合であれば、以下の言い回しも私は使いますけれども。

    単精度計算版 powf(2.0, 15.0) の戻り値の単精度浮動小数点数を、倍精度浮動小数点変数に代入する/渡す
    倍精度計算版 pow(2.0, 15.0) の戻り値の倍精度浮動小数点数を、単精度浮動小数点変数に代入する/渡す

     
    でも、そういう話をされているようには思えなくて、何と言うか、pow(2.0, 15.0) が単精度浮動小数点数でも倍精度浮動小数点数でもない「型」の戻り値をもっているような、そういうことをイメージされているのでしょうか?


    そして、

    > pow関数の中身のアルゴリズムの時点で(^の表現を使ってしまいますが)
    > pow(2,15) = e^(y*log(x))  = 32768とはちょっとずれた値
    > になってしまっている。

    数学的には、上の式の後半部分の途中までの「e^(y*log(x)) = 32768」は厳密に正しいです。他方で、 pow(2.0, 15.0) の浮動小数点ライブラリの実装上の計算結果の戻り値は、「pow(2.0, 15.0) = 32768とはちょっとずれた値」になってしまいます。また、「ちょっとずれた値」というのも、単精度計算版 powf(2.0, 15.0) と倍精度計算版 pow(2.0, 15.0) とでは、異なる「ちょっとずれた値」になっている、のです。


    あと、ググってみたのですけれども、Windowsの電卓アプリケーションは、「そこそこ多量」から「めちゃくちゃ多量」といった計算量の科学計算(それと3Dゲームでの計算なども)をするものではないから、1回の計算に0.01秒とか掛かっても構わないけど人間社会の常識に従った計算をして欲しい、という要請に対応する為に、もはや咄嗟にどういう内部表現になっているのか分からない、というようものになっているのかなぁ?という印象を持ちました(以下の記事の斜め読みで)。なので、IEEE単精度浮動小数点数やIEEE倍精度浮動小数点数の話に、それを持ち出して大丈夫かなぁ?と思ったりします。

    Windows10 の新しい電卓は計算結果が演算誤差により正しくない場合がある上に指数表示になって使い物にならないしもとに戻すこともできない
    qiita.com/Q11Q/items/a724e7f51b22beeecac6

    あと、Windowsの電卓アプリケーションでは、これに類することも関係するかも知れません。

    0.1は浮動小数点数で正確に表せないのに、printしたときに0.1と表示されるのはなぜか
    qiita.com/yucatio/items/c03def631bc5eaaa9887