当前位置: 华文问答 > 科学

这个求 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位),有效数字之后不知道是啥了,就可以理解为「??????????」,只要加和的时候碰到过这个位数,那一位就废了。