當前位置: 華文問答 > 科學

這個求 sin 的程式碼,感覺無懈可擊啊,為什麽會不對呢?

2022-12-02科學

試著用人話把這個大坑理個提綱【

「浮點數」

也就是定長科學計數法(再具體點,定「有效位數」),會「 大數吃小數 」,因為每步計算要「修約」到固定長度的有效數位,其中的出現數位最高位越高( 越大 ),被 丟掉 的小數位數就越高。

而泰勒展開在自變量大的時候,x^n初期增長就會比下面的n!快,雖然最終會被趕超,但是小數已經損失了,啥都補不回來。

部份解決的方法是pi reduction,就是整除2pi取余數後再算,眾所周知sin(x)=sin(2pi+x)。(還有更細,最後是正負pi/4)

然而這一步也有坑,就算你覺得只是輸入我們想要的有效位數,當輸入非常接近pi的倍數時,也就是輸出的值很接近0,它的小數位數會探到一個比你預料的還小的地方,最終會撞到無理數倍數在什麽時候接近有理數的問題:

The Difference Between x87 Instructions FSIN, FCOS, FSINCOS, and...

在這裏就算intel也作了妥協,如果要絕對精確(到固定位數),這種情況得提升內建的pi的精度到至少兩倍有效數位長,而且要提升到多高還不好預測(不嚴謹地說就是如果存在一個n使得n*pi的精確值在53位(如果long double就64位元)二進制有效位數小數尾數後可能有很多個0然後又繼續blablabla,那麽這個內建pi的精度就得達到53+這個「很多」再+53位才行,不然前面的pi減掉了剩下的0.0000000xxxxxx很小的數的有效數位就不對了)。

「為什麽sin(pi)會得到一個小數,把這個小數跟電腦裏的pi加起來會得到一個兩倍精度的pi?」

所以嵌入式爺傳統娛樂計畫之一是人均自擼三角函式【具體可以是查表、查表加擬合、cordic等等【不用弧度,而用整數2^n為圓周,這樣輸入角度的時候已經自動取模(實際是把modpi這事交給了使用者去做

已知輸入接近pi的倍數(或者就是接近0了),又想求精確值,那才是自己泰勒展開的時候。

另一類小數運算的風格是約定「保留多少位小數」(而非有效數位),也就是定點數,比如一個linux程式(最嚴格說是posix定義,然後常用gnu實作)「bc」,指定小數位數,不管整數部份飛起來多少位,小數點後保留多少位都是固定的(可用參數調,計算過程中固定),於是直接用這個公式也可以算出一個「尾數最多錯一兩位(此處省略一本數值計算.pdf)」結果。

這樣某種意義上的計算量就「不常數」的。

浮點數之「定長」,追求的是某種程度的「固定計算量、固定計演算法(不管你這數數量級多大)」,便於硬體實作,至少加減乘除平方根上還是符合這一認知的。

但對於超越函式們(雖然至今好像也就intel和某奇妙的stc32有硬體弧度sin),就在這個對pi取模的時候,如果輸入非常接近pi了,又要絕對精確,也難以固定(代價比不上收益)。所以實際上,它們計算的結果也是精確到「某個小數位」,而非總是精確到「某個數量有效位數」,intel整個文章就是在解釋這個事。

「來點我看得懂的」

CatastrophicCancellationVisualizer

js的「toFixed」是個好東西(其實也不是完全的好)(c也可以插一句也能看個大概

printf ( "@[email protected] \n " , term , sum );

原版程式碼計算sin(50)的過程打印展示出來大概是這樣:

term sum 50.00000000000000??? 50.00000000000000??? 20833.33333333333?????? -20783.33333333333?????? 2604166.666666666???????? 2583383.333333333???????? 155009920.6349206?????????? -152426537.3015872?????????? 5382288910.934743??????????? 5229862373.633156??????????? 122324747975.7896????????????? -117094885602.1564????????????? 1960332499612.013?????????????? 1843237614009.856?????????????? 23337291662047.77??????????????? -21494054048037.92??????????????? 214497166011468.5???????????????? 193003111963430.6???????????????? 1567961739849916.????????????????? -1374958627886485.????????????????? 9333105594344740.????????????????? 7958146966458254.????????????????? 4611218179024080?.????????????????? -3815403482378255?.????????????????? 1921340907926700??.????????????????? 1539800559688874??.????????????????? 6842382150736111??.????????????????? -5302581591047237??.????????????????? 2106644750842398???.????????????????? 1576386591737675???.????????????????? 5663023523769889???.????????????????? -4086636932032214???.????????????????? 1340677917559159????.????????????????? 9320142243559378???.????????????????? 2816550246973024????.????????????????? -1884536022617086????.????????????????? 5286318031105526????.????????????????? 3401782008488439????.????????????????? 8917540538302169????.????????????????? -5515758529813729????.????????????????? 1359381179619233?????.????????????????? 8078053266378601????.????????????????? 1881756893160621?????.?????????????????-1073951566522761?????.????????????????? 2375955673182602?????.????????????????? 1302004106659841?????.????????????????? 2747404802477570?????.?????????????????-1445400695817729?????.????????????????? 2920285716919186?????.????????????????? 1474885021101457?????.????????????????? 2863025212665868?????.?????????????????-1388140191564411?????.????????????????? 2597083828615628?????.????????????????? 1208943637051216?????.????????????????? 2186097498834703?????.????????????????? -9771538617834874????.????????????????? 1712169093698859?????.????????????????? 7350152319153720????.????????????????? 1250854101182685?????.????????????????? -5158388692673131????.????????????????? 8544085390592113????.????????????????? 3385696697918981????.????????????????? 5468564638115792????.????????????????? -2082867940196810????.????????????????? 3286397018098432????.????????????????? 1203529077901622????.????????????????? 1857981127373605????.????????????????? -6544520494719832???.????????????????? 9899728939543932???.????????????????? 3355208444824100???.????????????????? 4979742927335982???.????????????????? -1624534482511881???.????????????????? 2368599185376704???.????????????????? 7440647028648227??.????????????????? 1066936569989506???.????????????????? -3228718671246835??.????????????????? 4557999700912108??.????????????????? 1329281029665273??.????????????????? 1849237139285990??.????????????????? -5199561096207168?.????????????????? 7134402543541630?.????????????????? 1934841447334462?.????????????????? 2620629791192194?.????????????????? -6857883438577320.????????????????? 9175874619020288.????????????????? 2317991180442968.????????????????? 3065983232765399.????????????????? -747992052322431.5???????????????? 978671869498659.2???????????????? 230679817176227.7???????????????? 298739886904352.6???????????????? -68060069728124.87??????????????? 87289588272660.31??????????????? 19229518544535.43??????????????? 24437174768381.94??????????????? -5207656223846.507?????????????? 6560667624672.988?????????????? 1353011400826.480?????????????? 1690545151688.566?????????????? -337533750862.0859????????????? 418451770219.9421????????????? 80918019357.85620???????????? 99574474162.36962???????????? -18656454804.51342???????????? 22796353974.90147???????????? 4139899170.388042??????????? 5024765027.089901??????????? -884865856.7018594?????????? 1067100965.657896??????????? 182235108.9560366?????????? 218489141.2075954?????????? -36254032.25155875????????? 43159201.40794789????????? 6905169.156389147???????? 8230206.218144145???????? -1325037.061754998???????? 1516026.786425019???????? 190989.7246700208??????? 269909.3409815232??????? -78919.61631150235?????? 46471.99397064794?????? -32447.62234085441?????? 7742.235434267615????? -40189.85777512202?????? 1248.747650688324????? -38941.11012443370?????? 195.0924338658175???? -39136.20255829951?????? 29.53797751117634??? -39106.66458078833?????? 4.336168160771630?? -39111.00074894911?????? 0.6174766690549712? -39110.38327228005?????? 0.08533397858692250 -39110.46860625864?????? 0.01144992198729639 -39110.45715633666?????? 0.00149227426588682 -39110.45864861092?????? 0.00018899116842538 -39110.45845961976?????? 0.00002326789722562 -39110.45848288765?????? 0.00000278590723487 -39110.45848010174?????? 0.00000032451626536 -39110.45848042626?????? 0.00000003678989041 -39110.45848038947?????? 0.00000000406069431 -39110.45848039353?????? 0.00000000043652115 -39110.45848039309?????? 0.00000000004571860 -39110.45848039314?????? 0.00000000000466669 -39110.45848039313?????? 0.00000000000046440 -39110.45848039313?????? 0.00000000000004507 -39110.45848039313?????? 0.00000000000000427 -39110.45848039313?????? 0.00000000000000039 -39110.45848039313?????? 0.00000000000000004 -39110.45848039313??????

什麽叫「定長科學計數法」,就是有效數位長度一定(二進制53位,約等於十進制17位),有效數位之後不知道是啥了,就可以理解為「??????????」,只要加和的時候碰到過這個位數,那一位就廢了。