matlab - 为大 x 计算 4^x mod 2π

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 时，我的明确目标之一就是能够计算像 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 当然能够计算三角函数。但正如我们上面所展示的，我们应该只相信结果的前 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 类型的计算？因此，计算所需的功率，同时在每一步取模数。因此，以 double 完成:

X = 1;
for i = 1:1000
  X = mod(X*4,2*pi);
end
sin(X)
ans =
         0.955296299215251

啊，但请记住，真正的答案是-0.19033458127208318385994396068455455709388...

所以基本上没有什么重要的了。我们在那个计算中丢失了所有信息。正如我所说，重要的是要小心。

发生的事情是在该循环的每一步之后，我们在模数计算中都发生了微小的损失。但是然后我们将答案乘以 4，这导致误差增加了 4 倍，然后又增加了 4 倍，等等。当然，每一步之后，结果在数字末尾都会丢失一点点.最终的结果是完整的克拉波拉。

让我们看一下较小功率的操作，只是为了说服自己发生了什么。例如，在这里尝试 20 次方。使用 double ，

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

事实上，虽然 double 的结果与显示的最后一位数字一致，但那些 double 结果实际上都是错误的，超过了第 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 吗？不!

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