[過去ﾛｸﾞ] プログラミングのお題スレ Part21 (1002ﾚｽ)
上下前次1-新
抽出解除 ﾚｽ栞
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
915(14): デフォルトの名無しさん [sage] 2023/07/11(火) 15:24:57.95 ID:zBgA5NaS(1) AAS
お題：神経質なipow10()
整数を引数とし、その10のべき乗を浮動小数点数で返す関数を書け。
結果が浮動小数点数で厳密に表現出来ない場合(10の-1乗=0.1とか)、
浮動小数点数として、(結果を)厳密に最近値に丸めた値を返すこと。
917: 915 [sage] 2023/07/12(水) 01:10:01.78 ID:OtxkddAG(1) AAS
正の整数の範囲だけなら、CPythonのmath.pow()もそう悪くないんだ…
(10**23と10**210が外れ)。
外部ﾘﾝｸ:sagecell.sagemath.org
919(1): 915 [sage] 2023/07/12(水) 15:28:21.53 ID:6VP/LOOe(1) AAS
"1e23" とかの文字列を構成して、Cのstrtod()やsscand()的な関数を使って変換する分には誤差が出ないっぽいのですけど
(これの誤差が出るライブラリは駆逐されたらしい)、
これはこれで「負け」た気分になるのが…ｗ
920: 915 [sage] 2023/07/12(水) 15:35:55.97 ID:9l9n30C9(1) AAS
>>919
sscanfですたorz
923: 915 [sage] 2023/07/13(木) 05:18:24.50 ID:102DBCHI(1/2) AAS
何で「負け」なのかと云うと、
IEEE754的な意味での「再現性」を担保する為に
(binary64の数値を有効桁17桁で10進文字列に変換して、
それをstrtod()とかするとbit perfectで元に戻るってやつ)、
strtod()とかは内部で多倍長演算をやってるので…。
(glibcの実装(GNU MPに丸投げ？)と、ルーセントの中の人
が書いたやつ(そのほか)の2系統あり)。
…なんとかズル出来ないのか？ｗ
924: 915 [sage] 2023/07/13(木) 05:22:25.37 ID:102DBCHI(2/2) AAS
要するに牛刀割鶏感がでかい
いまさら気にしても仕方がなくもないけど
937: 915 [sage] 2023/07/14(金) 19:28:34.12 ID:Q8y8QZ/7(1) AAS
仕事で書いたら怒られるやつｗ
(プラットホームや、ライブラリのバージョンとかに依存する可能性があるから)
外部ﾘﾝｸ:sagecell.sagemath.org
944: 915 [sage] 2023/07/15(土) 10:11:42.89 ID:HHeu+C79(1/3) AAS
今からIEEE754に準拠してない処理系が制作される事は無いでしょうから
(ホビーや学習用途は別ね)、
今日の処理系で、2のべき乗での割り算で仮数部が変わるって事は
無いと思われます(アンダーフローやオーバーフローは除外)。
…IBM hexadecimal floating-pointって基数が16だったな。
(仮数部が変わらないのは16のべき乗での乗除算のみ)
あれを常用してる人ってまだ居るのか？ｗ
>>943943(1): デフォルトの名無しさん [] 2023/07/15(土) 08:05:25.03 ID:DpclFvGX(1) AAS
>>925
ツッコミどころが多すぎる
正しくはこうでは?
10^(-30)
=1.010_0010_0100_0010_0101_1111_...×2^(-100)
さらに、...部分は1111...と続くので、単精度にするなら切り上げになって、
1.010_0010_0100_0010_0110_0000×2^(-100)
よって、単精度での内部表示は
0_00011011_010_0010_0100_0010_0110_0000
整数0xa24260をセットして2^123で割ればいける
2のべき乗での割り算すら信用できないなら、ldexpあたりを使うのはどうか。pythonならこちら
(pythonは単精度が無いので結果は少しずれる)
from math import ldexp
print(ldexp(0xa24260, -123))
もしくは、バイト列経由で直接変換もあり。pythonならこちら
from struct import pack, unpack
print(unpack('f', pack('I', 0b0_00011011_010_0010_0100_0010_0110_0000))[0])
FORTRANならEQUIVALENCE文で…
948: 915 [sage] 2023/07/15(土) 18:15:18.86 ID:HHeu+C79(2/3) AAS
>>946946(1): デフォルトの名無しさん [] 2023/07/15(土) 16:07:59.58 ID:STkOcJIm(2/2) AAS
そもそも浮動小数点ライブラリなら
x^y = exp( y * log 2 )
で実装してるんじゃないかな
そもそものexpやlogが“最近値”を与えてくれることなど規格にはないやろ
それは数値計算の理論勉強してたらすぐわかる、メチャクチャ難しくてほとんど実装不能になってしまう
だから規格上は“±その周辺の値で表示できる値の中では最良”、すなわち“表現可能な上界の最小値または表現可能な下界の最大値”、それならギリギリ実用可能な速度くらいは出せるけど、それとて難しい、そこまで今のライブラリが保証してくれてるかあやしい
まぁその手のライブラリ使わずに“最近値”を出すことはメチャクチャ手間かかるけどできなくはない、しかしそれが限界、その値を変数に確実にセットする方法がない
規格で正確な値(を最近値に丸めたもの)を求められてるのは
加減乗除とsqrt()とfma()でしたっけ？
pow()は、glibcが誤差0.523ulpとか(0.5ulpが限界値)、
イイ線行ってるそうですね。
外部ﾘﾝｸ[pdf]:members.loria.fr
950: 915 [sage] 2023/07/15(土) 23:31:05.73 ID:HHeu+C79(3/3) AAS
>>949llvmの多くの項目がNA(non available)なのは、
実装されてないからではないかと。
この中では、OpenLibmがfdlibmの「血」を残してる方なのかな？
外部ﾘﾝｸ:www.netlib.org
952: 915 [sage] 2023/07/16(日) 01:29:25.23 ID:34tdZpnq(1/5) AAS
>>951951(1): デフォルトの名無しさん [] 2023/07/16(日) 00:22:02.04 ID:Cy0F/L8f(1) AAS
うん、多分実装されてない分は自分でやれと
実装されてるところは軒並み0.5
まぁこれも本当に確率100%で必ず最近値出すんではなくて外れる確率が0.000001%とかなんだろうな
例えば単精度演算は毎回倍精度にキャストして計算してから丸めればほぼほぼ最近値にならない確率0.0000001%とかになるからな
逆に言えば他のライブラリはその手のキャストも何もしないで“丸め誤差上等”でやってるんやろな
ちゃんとやろうとすると、CRlibmみたいに、部分的にでも
6倍精度とかで計算しないとならないっぽいです。
glibcのexp()は768bitで計算する場合があるそうで。
…世間的には、それ程の精度を求めない用途もあって、
実際、libmのhypot()(じゃなくてもいいのですが)が遅いので自分で
実装しますたと云う報告がネット上でも散見されます
(もちろん精度は落ちる)
954(1): 915 [sage] 2023/07/16(日) 10:57:33.80 ID:34tdZpnq(2/5) AAS
>>953953(1): デフォルトの名無しさん [] 2023/07/16(日) 08:49:44.74 ID:Q2/lNXhi(1) AAS
まじですか
768bitって何をどう計算してるんやろ？
結果が24bitで768bitなら求める答えの72倍のアキュムレーター使うのかな？
そんな事ないよな、アキュムレータは48bitくらいで36回計算するからのべ768bitとかいう事なのかな？
まぁそれだけの事やって「末尾の最後の1ビットの正確性を求める意味あんのか問題」は発生するわな
”みんなが使う汎用ライブラリ”だとどうしてもそういう“誰も求めてない精度”に無意味にこだわってしまう部分があるんかも
丸めて切り上げ/切り捨てのぎりぎりのとこを精査するのに
768bitで計算する場合もあるそうです(そういう
ぎりぎりのとこでない場合はやらないので、それなりに速い)。
956: 915 [sage] 2023/07/16(日) 12:36:42.33 ID:34tdZpnq(3/5) AAS
>>955955(1): デフォルトの名無しさん [] 2023/07/16(日) 11:52:54.00 ID:YiXq44MV(1) AAS
>>954
まじですか？
そんな最後の1ビット正確に決定するために768bitもの無駄な計算するくらいなら誤差±2^(-23)でいいからとっとと値返してくれた方がいいのに
なので当然通常“最後の1ビット”の正確性までは規格に入れてないんだよな、それでもライブラリ作成者は自己満のために誰も求めない“最後の1ビット”にこだわる
これは高校くらいまでの近似計算の「小数第××位まで求めよ」のノリが抜けてない事の現れでもあるんだよな
もう計算論の近似計算の世界ではゲームチェンジが起こってる事に気づけてない
誤差なし(±0.5ulp以内)を追及する人達(CORE-MATH)も居れば、
double-double演算とかの「飛び道具」を使うのをよしとしない人達も居て、
えーっと、みんな違ってみんないい(こなみ)
958: 915 [sage] 2023/07/16(日) 13:09:49.90 ID:34tdZpnq(4/5) AAS
>>957957(1): デフォルトの名無しさん [] 2023/07/16(日) 13:03:43.49 ID:nweoM3+S(1) AAS
まぁしかし“払うコスト”に対する“リターン”が少なすぎる気はする
例えば内部表現の精度が2^24まである単精度の計算をする場合、もちろん理想は“誤差±1/2×2^(-24)”で返してくれるとありがたい
しかし現実できるか？
例えばアキュムレータをを32ビット用意して計算する、マクローリン展開で求めるとして100回で打ち切るとする、打ち切り誤差は入力のサイズによるけど十分小さい、丸め誤差は2^(-32)×100で誤差トータルは2^(-25)程度、返すのが24bitだから丸められる2^8の可能性のうち真値のポジションから最大100ずれたところでウロチョロしてるわけでその“ウロチョロ”が二項分布、真値の分布が256個の箱の一様分布として外れ値になるの確率がどのくらいやろ？暗算ではできないけど言うほどない、この言うほどない誤差を返さないただそれだけのために768bitも計算繰り返すとかどうなんって感じ
768bitの話は、倍精度(53bit精度)のexp()の過去のバージョンでした。
glibcでも問題になったんでしょう(たぶん)
960: 915 [sage] 2023/07/16(日) 14:30:41.46 ID:34tdZpnq(5/5) AAS
>>959959(1): デフォルトの名無しさん [] 2023/07/16(日) 14:04:25.82 ID:JCFIIFPR(1) AAS
でしょうね
流石に768bitはない
せめてそこまで行けば完全に最近値が決定できるならともかくそこまで行っても最近値を100%決めるのは無理なんじゃないでしょうか？
24ビットの値返すのに768bitまで計算するなら744bit、2^744の可能性につきあってそこまで行っても両端の100通りの可能性は残り、最近値を返せない可能性は0にはなってない、まぁほとんど0だけど
多分llvmの実装は単精度でも最初からアキュムレータに64bit(レジスタ2個分)で計算するんやろな、あまりのビットが40bitあって100/2^40とか1000/2^40とかはほとんど0だからほぼ確率1で最近値返しますよ、それで十分でしょって実装なのではなかろかと、実際それで十分
最近のcpuはなんかレジスタ2個分に分けて計算してもレジスタ1個で計算するのと時間大して変わらないという話もききますしね
以前は最近値だったそうです(その代わり遅い)
外部ﾘﾝｸ:sourceware.org
上下前次1-新書関写板覧索設栞歴
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 0.195s*