C言語でべき乗表現

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

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

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

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

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

    他の言語のべき乗2項演算子とC言語のべき乗数学関数の違いが認識出来ていないのかなぁ、と思いました。

    Fortran (参考 www.ile.osaka-u.ac.jp/research/cmp/download/fortran90-text-v4.pdf#page=10)

     x ** y

    R (参考 ds.k.kyoto-u.ac.jp/wp-content/uploads/2021/10/slide-a02.pdf#page=15)

     x ^ y

    いろいろまとめ

    hydrocul.github.io/wiki/programming_languages_diff/number/pow.html

    こぼればなし

    blog.media.teu.ac.jp/2021/08/post-89817a.html

    数学関数としてのべき乗としては、ググってみたのですが、以下のウェブページを参照してみてはどうかな、と思いました。(CC-RLライブラリでの計算方法はともかくとして、浮動小数点計算をするものである、単純に x を y 回だけ掛け算するものでは無い、というイメージとしてです。)

    hiroyukichishiro.com/power-and-exponentiation-in-c-language/

    自作関数でべき乗の計算


  • NoMaYさん、 KatoNaganoさん、Sugachanceさん
    皆さんありがとうございます。NAKAです。

    元々、MicrochipのSAM71V(intは32bit)での疑問で、RL78で確認してもおかしかったので、
    お聞き致しました。

    元々、回路のハード設計出身でソフトは見様見真似でしっかり
    勉強してないことを反省してます。

    こういう場があるのは本当に助かります。

    どうも最近ルネサスのエディタの調子が?です。

    アルファベット打ってもカタカナになるし、ヘンテコに入力されるし...?((+_+))

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

    それで、どういうことだったのかというのは理解出来た、ですかね?

  • NoMaYさん、NAKAです。

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

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

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

    諸悪の原因はEXCELかな?


  • 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
     

Reply
  • 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
     

Children
  • 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
     

  • NoMaYさん、こんばんは。Sugachanceです。

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

    >これらの「渡す」

    ご指摘どおり変な感じで考えておりました...
    イメージとしては
    double pow (double x, double y);
    という関数に対して
    float ansf = (float)pow(2.0,15.0);
    のようにキャストして受け取るようなものを思い浮かべて書いていたのですが
    確かに変ですね。doubleのpow(2.0,15.0)をfloatにキャストしても
    ご指摘のあった通り32768になりますものね。

    CC-RLにおいては、コンパイルオプションのdouble型/long~のところの違いは
    ```
    倍精度計算版 pow(2.0, 15.0) の戻り値の倍精度浮動小数点数を、単精度浮動小数点変数に代入する
    ```
    から来ていたのかと思っていましたが、
    double型/long~をはいにした場合は
    単精度計算版 powf(2.0, 15.0) の戻り値の単精度浮動小数点数を単精度浮動小数点変数に代入しているようですね。


    >単精度計算版 powf(2.0, 15.0) と倍精度計算版 pow(2.0, 15.0) とでは、異なる「ちょっとずれた値」
    ここが本題の答えになるんですかね?
    →NAKAさんが最初に示された例では表記上はpow()だが、内部的にはpowf()が使用されていた。
    →小さめになるちょっとずれた値がpowf()から返されていた
    →32768を期待したが32767に丸められてしまった

    ちょっと個人的整理
    -------
    (1) pow(x,y)は数学的なxのy乗の答えとは異なる
    →計算アルゴリズムが異なる

    (2) pow(2.0, 15.0)とpowf(2.0, 15.0)は計算アルゴリズム上の理論値とは
    異なる「32768からちょっとずれた値」が返る
    ------

    (2)はdoubleのpowの場合
    double pow(double x, double y)
    {
        //この中の計算においてdoubleの表現の限界値が現れることがある
       // ln(x)をmathの関数から持ってきているのであればそれこそdouble型??
     //その影響で理論値からちょっとずれた値になる場合がある
     return ans;
    }
    という感じでしょうかね

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

    頂いたリプライの他の部分は後で見ますけれども、以下の部分に関して補足します。

    > CC-RLにおいては、コンパイルオプションのdouble型/long~のところの違いは ``` 倍精度計算版 pow(2.0, 15.0) の戻り値の倍精度浮動小数点数を、単精度浮動小数点変数に代入する ``` から来ていたのかと思っていましたが、 double型/long~をはいにした場合は 単精度計算版 powf(2.0, 15.0) の戻り値の単精度浮動小数点数を単精度浮動小数点変数に代入しているようですね。

    CS+ V8.09のCC-RLのヘルプには以下の通りに書かれていて、倍精度浮動小数点変数型サイズと倍精度浮動小数点計算数学ライブラリの組み合わせは、自由に出来そうですけれども、ヘッダファイル math.hにカラクリが入っているのです。


    7.2 ライブラリ・ファイル命名規則

    標準ライブラリを,アプリケーション内で使用するときは,関連するヘッダ・ファイルをインクルードして,ライブラリ関数を使用します。
    ランタイム・ライブラリ関数は,浮動小数点演算や整数演算を行うときに,CC-RLが自動的に呼び出すルーチンです。
    また,リンク・オプション(-library)で,これらのライブラリを参照するようにします。同一のプロジェクト内で使用するライブラリ・ファイルの種類は統一しなければなりません。

    ライブラリ・ファイルの命名規則は,次のようになっています。
    【V1.03以降】
    malloc系のライブラリ関数は通常用とセキュリティ機能用とで機能が異なるため,他の標準ライブラリとは別のライブラリとなります。malloc系ライブラリはRL78-S1,RL78-S2,RL78-S3共通です。
    rl78<muldiv><memory_model><float><standard/runtime><lang>.lib
    malloc_<secure>.lib (malloc系ライブラリ) 【V1.03以降】

    <muldiv>
    n :拡張命令なし,演算器なし(RL78-S1コア/RL78-S2コア用)
    c :乗除・積和演算器使用(RL78-S2コア用)
    e :乗除算拡張命令使用(RL78-S3コア用)

    <memory_model>
    m :スモール/ミディアム・モデル用

    <float>
    4 :単精度浮動小数点数
    8 :倍精度浮動小数点数(乗除算拡張命令ありデバイスのみ対応)

    <standard/runtime>
    s :標準ライブラリ
    r :ランタイム・ライブラリ

    <secure>
    n :通常用
    s :セキュリティ機能用 【Professional版のみ】

    <lang>
    なし :C90用
    99 :C99用 【V1.07以降】


    math.h

    #if defined(__DBL4)
    #define acos(_x)             acosf((_x))
    #define asin(_x)             asinf((_x))
    #define atan(_x)             atanf((_x))
    #define atan2(_y, _x)        atan2f((_y), (_x))
    #define cos(_x)              cosf((_x))
    #define sin(_x)              sinf((_x))
    #define tan(_x)              tanf((_x))
    #define cosh(_x)             coshf((_x))
    #define sinh(_x)             sinhf((_x))
    #define tanh(_x)             tanhf((_x))
    #define exp(_x)              expf((_x))
    #define frexp(_value, _exp)  frexpf((_value), (_exp))
    #define ldexp(_x, _exp)      ldexpf((_x), (_exp))
    #define log(_x)              logf((_x))
    #define log10(_x)            log10f((_x))
    #define modf(_value, _iptr)  modff((_value), (_iptr))
    #define pow(_x, _y)          powf((_x), (_y))
    #define sqrt(_x)             sqrtf((_x))
    #define ceil(_x)             ceilf((_x))
    #define fabs(_x)             fabsf((_x))
    #define floor(_x)            floorf((_x))
    #define fmod(_x, _y)         fmodf((_x), (_y))
    #if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)  /* C99 */
    #define acosh(_x)            acoshf((_x))
    #define asinh(_x)            asinhf((_x))
    #define atanh(_x)            atanhf((_x))
    #define log1p(_x)            log1pf((_x))
    #define scalbn(_x, _n)       scalbnf((_x), (_n))
    #define scalbln(_x, _n)      scalblnf((_x), (_n))
    #define nearbyint(_x)        nearbyintf((_x))
    #define rint(_x)             rintf((_x))
    #define lrint(_x)            lrintf((_x))
    #define llrint(_x)           llrintf((_x))
    #define round(_x)            roundf((_x))
    #define lround(_x)           lroundf((_x))
    #define llround(_x)          llroundf((_x))
    #define trunc(_x)            truncf((_x))
    #define copysign(_x, _y)     copysignf((_x), (_y))
    #define nan(_tagp)           nanf((_tagp))
    #define fdim(_x, _y)         fdimf((_x), (_y))
    #define fmax(_x, _y)         fmaxf((_x), (_y))
    #define fmin(_x, _y)         fminf((_x), (_y))
    #endif  /* C99 */
    #endif  /* __DBL4 */

     
    [追記]

    なお、コンパイルオプションの説明にも以下の通りに書かれています。


    -dbl_size

    double型とlong double型の解釈を変更します。

    [指定形式]

    -dbl_size={4|8}

    - 省略時解釈
    double型,およびlong double型を共に単精度浮動小数点型として扱います(-dbl_size=4オプションの指定と同じです)。

    [詳細説明]

    - double型とlong double型の解釈を変更します。
    - パラメータに4を指定した場合,double型,およびlong double型を単精度浮動小数点型として扱います。
    - パラメータに8を指定した場合,double型,およびlong double型を倍精度浮動小数点型として扱います。
    - パラメータに4,8を指定しない場合は,エラーとなります。
    - -cpu=S1オプションと-dbl_size=8オプションを同時に指定した場合は,エラーとなります。
    - -cpu=S2オプションと-dbl_size=8オプションを同時に指定した場合は,エラーとなります。
    - -dbl_size=4オプション指定時にdouble型用の標準ライブラリ関数を呼んだ場合,float型用の標準ライブラリ関数に置き換えます。
    - 本オプションは定義済みマクロに影響を与えます。


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

    解説ありがとうございました。

    >- -dbl_size=4オプション指定時にdouble型用の標準ライブラリ関数を呼んだ場合,float型用の標準ライブラリ関数に置き換えます。

    これで最初の方のpowfの値になったのですね。