Calculate 4 ^ x mod 2π for large x

Question

I need to calculate sin(4^x) with x> 1000 in Matlab, basically sin(4^x mod 2π) . As the values ​​inside the sin function become very large, Matlab returns infinity for 4^1000 . How can I calculate this efficiently? I prefer to avoid large data types.

I think converting to something like sin(n*π+z) may be a possible solution.

Answer 1

You need to be careful, as there will be a loss of accuracy. The sin function is periodic, but 4 ^ 1000 is a large number. So efficiently, we subtract the multiple of 2 * pi to move the argument to the interval [0.2 * pi).

4 ^ 1000 is approximately 1e600, a really large number. So I will do my calculations using a high precision floating point tool in MATLAB . (Actually, one of my explicit goals when I wrote HPF was to calculate a number like sin (1e400). Even if you are doing something for fun, doing it right still makes sense.) Q in this case, since I know that the power of interest to us is approximately 1e600, then I will do my calculations for more than 600 digits of accuracy, expecting me to lose 600 digits by subtractive cancellation. This is a serious subtraction question. I'm thinking about it. This module operation actually represents the difference between two numbers that will be identical for the first 600 digits or so!

 X = hpf(4,1000); X^1000 ans = 

What is the nearest multiple 2 * pi that does not exceed this number? We can get this with a simple operation.

 twopi = 2*hpf('pi',1000); twopi*floor(X^1000/twopi) ans = .6669043995793459614134256945369645075601351114240611660953769955068077703667306957296141306508448454625087552917109594896080531977700026110164492454168360842816021326434091264082935824243423723923797225539436621445702083718252029147608535630355342037150034246754736376698525786226858661984354538762888998045417518871508690623462425811535266975472894356742618714099283198893793280003764002738670747 

As you can see, the first 600 digits were the same. Now that we subtract two numbers,

 X^1000 - twopi*floor(X^1000/twopi) ans = 3.333095600420654038586574305463035492439864888575938833904623004493192229633269304270385869349155154537491244708289040510391946802229997388983550754583163915718397867356590873591706417575657627607620277446056337855429791628174797085239146436964465796284996575324526362330147421377314133801564546123711100195458248112849130937653757418846473302452710564325738128590071680110620671999623599726132925263826 

That's why I called it a serious subtraction problem. Two numbers were the same for many numbers. Even with 1000 precision digits, we lost a lot of digits. When you subtract two numbers, although we have a result with 1000 digits, now only 400 higher order digits matter.

HPF can, of course, calculate the trigger function. But, as we showed above, we should trust approximately the first 400 digits of the result. (For some problems, the local form of the sin function may cause us to lose more digits than this.)

 sin(X^1000) ans = -0. 

Am I right and we can’t trust all these numbers? I will do the same calculation, once every 1000 digits of accuracy, then a second time to 2000 digits. Calculate the absolute difference, then take log10. The result of the digit in 2000 will be our reference as essentially accurate compared to the result of 1000 digits.

 double(log10(abs(sin(hpf(4,[1000 0])^1000) - sin(hpf(4,[2000 0])^1000)))) ans = -397.45 

Oh. So, out of those 1000 digits of accuracy with which we started, we lost 602 digits. The last 602 digits are non-zero as a result, but still complete garbage. It was as I expected. Just because your computer reports high accuracy, you need to know when not to trust it.

Can we do the calculations without resorting to a high-precision tool? Be careful. For example, suppose we use the type of computation powermod? Thus, calculate the desired power by taking the module at each step. Thus, it is performed with double precision:

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

Ah, but remember that the true answer was -0.19033458127208318385994396068455455709388 ...

Thus, essentially nothing substantial remained. We have lost all our information in this calculation. As I said, it is important to be careful.

What happened after each step in this cycle, we suffered a tiny loss of module calculation. But then we multiplied the answer by 4, which caused the error to grow 4 times, and then another factor of 4, etc. And, of course, after each step, the result loses a tiny bit at the end of the number. The end result was a complete crapola.

Let's look at an operation for less power, just to make sure what happened. For example, try the 20th degree. Using double precision,

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

Now, use the loop in computing powermod, taking the mod after each step. Essentially, this discards a multiple of 2 * pi after each step.

 X = 1; for i = 1:20 X = mod(X*4,2*pi); end X X = 3.55938555711037 

But is there a right meaning? Again, I will use hpf to calculate the correct value, showing the first 20 digits of this number. (Since I did the calculation in 50 common digits, I fully trust the first 20 of them.)

 mod(hpf(4,[20,30])^20,2*hpf('pi',[20,30])) ans = 3.5593426962577983146 

In fact, although the results are in double precision consistent with the last figure shown, these double results were actually erroneous for the 5th significant digit. As it turned out, for this cycle we need to carry more than 600 digits of accuracy to get the result of any value.

Finally, in order to completely kill this dead horse, we may ask if a more efficient powermod calculation could be done. That is, we know that 1000 can be decomposed into binary form (use dec2bin) as:

 512 + 256 + 128 + 64 + 32 + 8 ans = 1000 

Can we use the re-squaring scheme to expand this large power with fewer multiplications and therefore cause less accumulated error? Essentially, we could try to calculate

 4^1000 = 4^8 * 4^32 * 4^64 * 4^128 * 4^256 * 4^512 

However, do this by repeating square 4, and then take the mod after each operation. However, this does not work, since the modulo operation will only remove integer multiples of 2 * pi. In the end, the mod is really designed to work with integers. See what happens. We can express 4 ^ 2 as:

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

Can we just compose the remainder and then take the mod again? NOT!

 mod(3.43362938564083^2,2*pi) ans = 5.50662545075664 mod(4^4,2*pi) ans = 4.67258771281655 

We can understand what happened when we expand this form:

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

What do you get when you remove INTEGER shorts from 2 * pi? You need to understand why the direct loop allowed me to remove the integer multiples of 2 * pi, but this squaring operation does not work. Of course, the direct cycle also failed due to numerical problems.

Answer 2

First, I would reconsider the question as follows: calculate 4 ^ 1000 modulo 2pi. So, we divided the problem into two parts.

Use some math tricks:

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

Therefore, you can simply multiply 4 times individually and calculate mod 2pi at each stage. The real question to ask, of course, is the accuracy of this thing. This requires careful mathematical analysis. This may or may not be complete shit.

Answer 3

After Paul hinted at mod, I found a mod feature for high permissions at mathwors.com. bigmod (number, power, modulo) CANNOT calculate 4 ^ 4000 mod 2π. Because it just works with integers as modulo, not decimals.

This operator is no longer suitable: sin(4^x) is sin(bidmod(4,x,2*pi)) .