我认为将其转化为类似于sin(n*π+z)的形式可能是一种解决方案。
需要注意的是,精度会有所损失。sin函数是周期性的,但4^1000是一个很大的数字。因此,为了将参数移到区间[0,2*pi),我们实际上要减去2*pi的倍数。
4^1000大约是1e600,是一个非常大的数字。因此,我将使用MATLAB中的高精度浮点工具HPF进行计算。(事实上,我编写HPF的一个显式目标之一就是能够计算像sin(1e400)这样的数字。即使您只是为了好玩而做某件事,正确地做也很有意义。)在这种情况下,由于我知道我们感兴趣的幂大约是1e600,因此我将以超过600位的精度进行计算,预计通过减法消除的位数约为600位。这是一个巨大的减法消除问题。想想看吧。那个模操作实际上是两个数字之间的差,它们的前600位左右是相同的!
X = hpf(4,1000);
X^1000
ans =
114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029376
twopi = 2*hpf('pi',1000);
twopi*floor(X^1000/twopi)
ans = 114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029372.6669043995793459614134256945369645075601351114240611660953769955068077703667306957296141306508448454625087552917109594896080531977700026110164492454168360842816021326434091264082935824243423723923797225539436621445702083718252029147608535630355342037150034246754736376698525786226858661984354538762888998045417518871508690623462425811535266975472894356742618714099283198893793280003764002738670747
如您所见,前600位数字是相同的。现在,当我们减去这两个数字时,
X^1000 - twopi*floor(X^1000/twopi)
ans =
3.333095600420654038586574305463035492439864888575938833904623004493192229633269304270385869349155154537491244708289040510391946802229997388983550754583163915718397867356590873591706417575657627607620277446056337855429791628174797085239146436964465796284996575324526362330147421377314133801564546123711100195458248112849130937653757418846473302452710564325738128590071680110620671999623599726132925263826
sin(X^1000)
ans =
-0.1903345812720831838599439606845545570938837404109863917294376841894712513865023424095542391769688083234673471544860353291299342362176199653705319268544933406487071446348974733627946491118519242322925266014312897692338851129959945710407032269306021895848758484213914397204873580776582665985136229328001258364005927758343416222346964077953970335574414341993543060039082045405589175008978144047447822552228622246373827700900275324736372481560928339463344332977892008702220160335415291421081700744044783839286957735438564512465095046421806677102961093487708088908698531980424016458534629166108853012535493022540352439740116731784303190082954669140297192942872076015028260408231321604825270343945928445589223610185565384195863513901089662882903491956506613967241725877276022863187800632706503317201234223359028987534885835397133761207714290279709429427673410881392869598191090443394014959206395112705966050737703851465772573657470968976925223745019446303227806333289071966161759485260639499431164004196825
所以,我是正确的,我们不能信任所有这些数字吗?我将进行相同的计算,一次在1000位数字精度下,然后第二次在2000位数字精度下。计算绝对差异,然后取log10。2000位数字的结果将作为我们的参考,与1000位数字的结果相比几乎完全精确。
double(log10(abs(sin(hpf(4,[1000 0])^1000) - sin(hpf(4,[2000 0])^1000))))
ans =
-397.45
X = 1;
for i = 1:1000
X = mod(X*4,2*pi);
end
sin(X)
ans =
0.955296299215251
mod(4^20,2*pi)
ans =
3.55938555711037
X = 1;
for i = 1:20
X = mod(X*4,2*pi);
end
X
X =
3.55938555711037
但这是正确的值吗?再次使用hpf计算正确的值,并显示该数字的前20位数字。(由于我已经在50个总数字中进行了计算,我绝对信任其中的前20个数字。)
mod(hpf(4,[20,30])^20,2*hpf('pi',[20,30]))
ans =
3.5593426962577983146
512 + 256 + 128 + 64 + 32 + 8
ans =
1000
4^1000 = 4^8 * 4^32 * 4^64 * 4^128 * 4^256 * 4^512
然而,可以通过反复平方4并在每次操作后取模来实现这一点。然而,这种方法失败了,因为模运算只能去除2 * pi的整数倍。毕竟,模运算确实是设计用于整数的。所以看看会发生什么。我们可以将4^2表示为:
4^2 = 16 = 3.43362938564083 + 2*(2*pi)
mod(3.43362938564083^2,2*pi)
ans =
5.50662545075664
mod(4^4,2*pi)
ans =
4.67258771281655
4^4 = (4^2)^2 = (3.43362938564083 + 2*(2*pi))^2
如果你移除2*pi的整数倍,会得到什么?你需要理解为什么直接循环可以移除2*pi的整数倍,但上面的平方操作却不行。当然,直接循环也因为数值问题而失败了。
sin(4^4) ~ -0.9992
。sin(4^1000) ~ -0.19033
,我想是这样的。 - DSM在Pavel的提示下,我在mathwors.com上找到了一个高次幂的mod函数。bigmod(number,power,modulo)不能计算4^4000 mod 2π。因为它只能处理整数模数而不是小数。
这个语句已经不正确了:sin(4^x)
是sin(bidmod(4,x,2*pi))
。
sin(2*pi*bigmod(4, 100, 2*pi))
,结果不为零。我相当确定在Matlab中也是如此。 - db12344^100 mod 2pi
不是整数。 - Danicasin(bigmod(4*2*pi, 100, 2*pi))
,这会得到零,正如它应该做的那样。 - db1234bigmod(4*2*pi + 0.001, 1.123, 2*pi)
的结果与mod((4*2*pi + 0.001)^1.123, 2*pi)
不同。即使对于小的x值也是如此。我想像这样的bigmod(4, 100.132, 2*pi)
可能由于这个观察结果而不准确。 - db1234