6

我需要sin(4^x)在 Matlab 中用 x > 1000 进行计算，基本上是sin(4^x mod 2π)因为 sin 函数内部的值变得非常大，Matlab 为4^1000. 我怎样才能有效地计算这个？我更喜欢避免使用大数据类型。

我认为转换为类似的东西sin(n*π+z)可能是一个可能的解决方案。

4

3 回答 3

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 当然可以计算三角函数。但正如我们上面所展示的，我们应该只信任结果的大约前 400 位数字。（在某些问题上，sin 函数的局部形状可能会导致我们丢失更多的数字。）

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

啊。因此，在我们一开始的 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 吗？不！

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

于 2012-12-02T23:37:36.297 回答
5

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

使用一些数学技巧：

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

因此，您可以自己乘以 4 并在每个阶段计算 mod 2pi。当然，真正要问的问题是这件事的精确度是多少。这需要仔细的数学分析。它可能完全是废话，也可能不是。

于 2012-12-02T22:31:09.343 回答
3

按照 Pavel 对 mod 的提示，我在 mathwors.com 上找到了一个用于高功率的 mod 函数。 bigmod(number,power,modulo)不能计算 4^4000 mod 2π。因为它只使用整数作为模数而不是小数。

这种说法不再正确：sin(4^x)is sin(bidmod(4,x,2*pi)).

于 2012-12-02T22:57:28.673 回答