计算大x下的4^x mod 2π。

Question

计算大x下的4^x mod 2π。

matlabmathtrigonometrymodulolargenumber

6

我需要在Matlab中计算x>1000时的sin(4^x)，即基本上是sin(4^x mod 2π)。由于sin函数内部的值变得非常大，Matlab会返回无穷大的结果，如4^1000。我该如何高效地计算它？我希望避免使用大数据类型。

我认为将其转化为类似于sin(n*π+z)的形式可能是一种解决方案。

- Bene

1

你可能也想查看 mathoverflow 上关于数学的问题 - 尽管这个问题与编程有关。 - Jeff

3

我有一种预感，任何方法都会因为双精度算术而导致垃圾数据：例如，在Matlab（或C）中计算sin（2pi10^i），其中i = 1:100，会得到垃圾数据，同样地，在计算mod（2pi10^i + 0.0001，2*pi）时，其中i = 1:100，也会得到垃圾数据。 - db1234

这也是我的担忧。因此，在使用循环之前，我通常会尝试找到数学解决方案。谢谢。 - Bene

@Jeff，这个问题对于数学溢出（mathoverflow）来说是不相关的，因为它是针对“研究级别的数学问题”，也就是在撰写或阅读文章或研究生级别书籍时遇到的那些问题。你可能在想[math.se]，它是为“任何级别的数学学习者和相关领域的专业人士”而设立的。 - AakashM

@AakashM 哦，是的，那就是我想的。我的错... - Jeff

3个回答

5

我首先要重新定义问题：计算4^1000模2pi的余数。因此我们将问题分为两部分。

使用一些数学技巧：

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

因此，你可以将4乘以自身多次，并在每个阶段计算模2pi。当然，真正需要问的问题是这个东西的精度是多少。这需要仔细的数学分析。它可能是完全没用的。

- Pavel Radzivilovsky

我有一个预感，这种方法或任何方法都会因为双精度算术而导致垃圾结果。例如，对于i = 1:100计算sin(2pi10^i)会给出垃圾结果，同样地，计算mod(2pi10^i + 0.0001, 2*pi)也是如此。 - db1234

是的，这是可能的。但需要小心处理机器算术。 - Pavel Radzivilovsky

1

我找到了一个处理幂和大数字的模函数： http://www.mathworks.com/matlabcentral/fileexchange/7908-big-modulo-function/content/bigmod.m这可能是一个可行的解决方案。 - Bene

看我的回答。我展示了这里建议的过程会导致垃圾结果。 - user85109

3

在Pavel的提示下，我在mathwors.com上找到了一个高次幂的mod函数。bigmod(number,power,modulo)不能计算4^4000 mod 2π。因为它只能处理整数模数而不是小数。

这个语句已经不正确了：sin(4^x)是sin(bidmod(4,x,2*pi))。

- Bene

我不信任这个。在Octave（Matlab的开源版本）中，我尝试使用bigmod.m例程，输入sin(2*pi*bigmod(4, 100, 2*pi))，结果不为零。我相当确定在Matlab中也是如此。 - db1234

1

除非我大错特错，它不应该是零... 4^100 mod 2pi 不是整数。 - Danica

糟糕...犯了愚蠢的错误，谢谢@Dougal，我想测试应该是sin(bigmod(4*2*pi, 100, 2*pi))，这会得到零，正如它应该做的那样。 - db1234

1

@Bene，不确定您是否对非整数x感兴趣，但请注意bigmod(a，x，m)似乎无法处理非整数x：bigmod(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

1

顺便说一句，我最终解释了（在我的答案中）为什么大模数计算对于非整数参数会失败。 - user85109

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user85109 · Accepted Answer

需要注意的是，精度会有所损失。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

什么是小于该数字但最接近2 * π的倍数？我们可以通过简单的运算得到这个数字。

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倍，依此类推。当然，在每个步骤之后，结果在数字末尾失去了一点点。最终结果是完全无用的。

让我们看看更小幂次的操作，只是为了让自己相信发生了什么。例如，尝试20次方。使用双精度，

mod(4^20,2*pi)
ans =
          3.55938555711037

现在，在幂模计算中使用循环，在每一步之后进行取模。本质上，这会在每一步之后丢弃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位的精度。

最后，要彻底解决这个问题，我们可以问一下是否可以进行更好的幂模运算。也就是说，我们知道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的整数倍。毕竟，模运算确实是设计用于整数的。所以看看会发生什么。我们可以将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的整数倍，但上面的平方操作却不行。当然，直接循环也因为数值问题而失败了。