计算4 ^在x mod2π为大的x

Question

我需要计算sin(4^x)与在Matlab X> 1000，用的是基本上sin(4^x mod 2π)由于sin函数变得非常大内的值，MATLAB返回无限为4^1000。我如何有效地计算这个？ 我宁愿避免大数据类型。

我认为对sin(n*π+z)这样的转变可能是一种可能的解决方案。

Answer 1

你需要小心，因为会有精度损失。 sin函数是周期性的，但4^1000是一个很大的数字。所以有效地，我们减去2 * pi的倍数以将参数移动到区间[0,2 * pi）。

4^1000大概是1e600，是一个非常大的数字。所以我会用我的high precision floating point tool in MATLAB来做我的计算。 （实际上，当我写HPF时，我的一个明确目标就是能够计算出一个类似于罪的数字（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个数字是相同的。现在，当我们减去两个数字时，这就是为什么我将它称为大量的减法取消问题的原因。这两个数字对于许多数字都是相同的。即使携带1000个数字的准确性，我们也失去了许多数字。当您扣除这两个数字时，即使我们携带1000位数的结果，现在只有最高位的400位数才有意义。

当然，HPF能够计算trig函数。但正如我们上面所显示的那样，我们只应该大致相信结果的前400位数字。 （在一些问题上，sin函数的局部形状可能会导致我们失去更多数字。）

sin(X^1000) 
ans = 
-0.19033458127208318385994396068455455709388374041098639172943768418947125138650234240955423917696880832346734715448603532912993423621761996537053192685449334064870714463489747336279464911185192423229252660143128976923388511299599457104070322693060218958487584842139143972048735807765826659851362293280012583640059277583434162223469640779539703355744143419935430600390820454055891750089781440474478225522286222463738277009002753247363724815609283394633443329778920087022201603354152914210817007440447838392869577354385645124650950464218066771029610934877080889086985319804240164585346291661088530125354930225403524397401167317843031900829546691402971929428720760150282604082313216048252703439459284455892236101855653841958635139010896628829034919565066139672417258772760228631878006327065033172

所以我是对的，我们不能相信所有这些数字？我将进行相同的计算，一次以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 

事实上，虽然在双精度结果同意所示，这些双最后一位结果在第五位有效数字后实际上都是错误的。事实证明，我们仍然需要为这个循环携带超过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来完成此操作，然后在每次操作之后取模。然而，这失败了，因为模操作将仅移除2 * pi的整数倍。毕竟，mod确实是设计用于整数。所以看看会发生什么。我们可以将4^2表示为：

4^2 = 16 = 3.43362938564083 + 2*(2*pi) 

但是，我们可以将剩下的部分平方，然后再次采用mod？没有！

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的整数倍，但上面的平方运算不能。当然，由于数字问题，直接循环也失败了。

Answer 2

我首先将问题重新定义如下：计算4^1000模2pi。所以我们把问题分成两部分。

使用一些数学挂羊头卖狗肉：

(a+2pi*K)*(b+2piL) = ab + 2pi*(garbage)

因此，你可以乘4多次受到自身和计算MOD二皮音乐的每一个阶段。真正要问的问题当然是这件事的精确程度。这需要仔细的数学分析。它可能或可能不是一个完全废话。

Answer 3

继Pavel的提示与MOD之后，我在mathwors.com上发现了一个用于高功率的mod函数。 bigmod(number,power,modulo)无法计算4^4000 mod2π。因为它只用整数作为模而不用小数。

此声明不再正确：sin(4^x)为sin(bidmod(4,x,2*pi))。