Ben*_*ene 6 math matlab largenumber trigonometry modulo
我需要sin(4^x)在Matlab中用x> 1000 进行计算,基本上是sin(4^x mod 2?)因为sin函数内部的值变得非常大,所以Matlab返回无穷大4^1000.我怎样才能有效地计算出来？我更喜欢避免大数据类型.
我认为转变为类似的东西sin(n*?+z)可能是一种可能的解决方案.
小智 13
你需要小心,因为精度会下降.sin函数是周期性的,但是4 ^ 1000是一个很大的数字.因此,我们有效地减去2*pi的倍数以将参数移动到区间[0,2*pi]中.
4 ^ 1000大约是1e600,这是一个非常大的数字.所以我将在MATLAB中使用我的高精度浮点工具进行计算.(事实上,当我编写HPF时,我明确的目标之一是能够计算像sin(1e400)这样的数字.即使你正在做一些有趣的事情,做正确的事情仍然有意义.)在这种情况下,因为我知道我们感兴趣的功率大约是1e600,那么我将以超过600位的精度进行计算,期望通过减法消除我将丢失600位数.这是一个巨大的减法取消问题.想一想.该模数操作实际上是两个数字之间的差异,对于前600个数字左右相同!
X = hpf(4,1000);
X^1000
ans =
114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029376
2*pi的最接近倍数是多少,不超过这个数字？我们可以通过简单的操作来实现.
twopi = 2*hpf('pi',1000);
twopi*floor(X^1000/twopi)
ans = 114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029372.6669043995793459614134256945369645075601351114240611660953769955068077703667306957296141306508448454625087552917109594896080531977700026110164492454168360842816021326434091264082935824243423723923797225539436621445702083718252029147608535630355342037150034246754736376698525786226858661984354538762888998045417518871508690623462425811535266975472894356742618714099283198893793280003764002738670747
如您所见,前600位数字相同.现在,当我们减去这两个数字时,
X^1000 - twopi*floor(X^1000/twopi)
ans =
3.333095600420654038586574305463035492439864888575938833904623004493192229633269304270385869349155154537491244708289040510391946802229997388983550754583163915718397867356590873591706417575657627607620277446056337855429791628174797085239146436964465796284996575324526362330147421377314133801564546123711100195458248112849130937653757418846473302452710564325738128590071680110620671999623599726132925263826
这就是我将其称为大量减法取消问题的原因.对于许多数字,这两个数字是相同的.即使携带1000位精度,我们也失去了很多数字.当你减去这两个数字时,即使我们携带1000位数的结果,现在只有最高位400位数才有意义.
HPF当然能够计算trig功能.但正如我们上面所示,我们应该只信任结果的前400位数字.(在某些问题上,sin函数的局部形状可能会导致我们失去更多数字.)
sin(X^1000)
ans =
-0.1903345812720831838599439606845545570938837404109863917294376841894712513865023424095542391769688083234673471544860353291299342362176199653705319268544933406487071446348974733627946491118519242322925266014312897692338851129959945710407032269306021895848758484213914397204873580776582665985136229328001258364005927758343416222346964077953970335574414341993543060039082045405589175008978144047447822552228622246373827700900275324736372481560928339463344332977892008702220160335415291421081700744044783839286957735438564512465095046421806677102961093487708088908698531980424016458534629166108853012535493022540352439740116731784303190082954669140297192942872076015028260408231321604825270343945928445589223610185565384195863513901089662882903491956506613967241725877276022863187800632706503317201234223359028987534885835397133761207714290279709429427673410881392869598191090443394014959206395112705966050737703851465772573657470968976925223745019446303227806333289071966161759485260639499431164004196825
我是对的,我们不能相信所有这些数字？我将进行相同的计算,一次是1000位精度,然后是2000位数的第二次.计算绝对差值,然后取log10.与1000位数结果相比,2000位数结果将是我们的参考.
double(log10(abs(sin(hpf(4,[1000 0])^1000) - sin(hpf(4,[2000 0])^1000))))
ans =
-397.45
啊.因此,在我们开始使用的1000位数精度中,我们丢失了602位数.结果中的最后602个数字为非零,但仍然是完全垃圾.这是我所期待的.仅仅因为您的计算机报告的精度很高,您需要知道何时不信任它.
我们可以在不使用高精度工具的情况下进行计算吗？小心.例如,假设我们使用powermod类型的计算？因此,计算所需的功率,同时在每一步取模数.因此,以双精度完成:
X = 1;
for i = 1:1000
X = mod(X*4,2*pi);
end
sin(X)
ans =
0.955296299215251
啊,但请记住,真正的答案是-0.19033458127208318385994396068455455709388 ...
所以基本上没有任何重要意义.我们在该计算中丢失了所有信息.正如我所说,小心谨慎是很重要的.
发生的事情是在该循环中的每个步骤之后,我们在模数计算中引起了微小的损失.但后来我们将答案乘以4,这导致误差增长4倍,然后是4倍等等.当然,在每一步之后,结果在数字末尾会丢失一点点.最终结果是完整的crapola.
让我们看一下较小功率的操作,只是为了说服自己发生了什么.例如,尝试20次幂.使用双精度,
mod(4^20,2*pi)
ans =
3.55938555711037
现在,在powermod计算中使用一个循环,在每一步之后取mod.基本上,这会在每一步之后丢弃2*pi的倍数.
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
实际上,虽然双精度的结果与所示的最后一位数相符,但这些双重结果在第5位有效数字之后实际上都是错误的.事实证明,我们仍然需要为此循环提供超过600位的精度,以产生任何重要的结果.
最后,为了完全杀死这匹死马,我们可能会问是否可以进行更好的powermod计算.也就是说,我们知道1000可以分解成二进制形式(使用dec2bin):
512 + 256 + 128 + 64 + 32 + 8
ans =
1000
我们是否可以使用重复的平方方案来扩展具有较少乘法的大功率,从而减少累积误差？从本质上讲,我们可能会尝试计算
4^1000 = 4^8 * 4^32 * 4^64 * 4^128 * 4^256 * 4^512
但是,通过重复平方4,然后在每次操作后取出mod来做到这一点.但是,这会失败,因为模运算只会删除2*pi的整数倍.毕竟,mod真的是专门用于整数.那么看看会发生什么.我们可以表达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的INTEGER倍数时你会得到什么？你需要理解为什么直接循环允许我删除2*pi的整数倍,但上面的平方操作却没有.当然,由于数值问题,直接循环也失败了.
我首先将问题重新定义如下:计算4 ^ 1000模2pi.所以我们将问题分成两部分.
使用一些数学技巧:
(a+2pi*K)*(b+2piL) = ab + 2pi*(garbage)
因此,你可以自己多次乘以4并在每个阶段计算mod 2pi.当然,要问的真正问题是这件事的精确性.这需要仔细的数学分析.它可能是也可能不是完全废话.