18

For several days, I've been wondering how it would be possible to compute the sine of huge numbers like 100000! (radians). I obviously don't use double but cpp_rational from the boost multiprecision library. But I can't simply do 100000! mod 2pi and then use the builtin function sinl (I don't need more than 10 decimal digits..) as I'd need several million digits of pi to compute this accurately.

Is there any way to achieve this?

Oliver Borchert
  • 297
  • 5
    Why? Why is this necessary? – Jared Jul 05 '16 at 09:06
  • 1
    You chould try to use the following relationship: $\Gamma(n+0.5)=\frac{2n!}{n!4^n}\sqrt{\pi}$ – MrYouMath Jul 05 '16 at 09:08
  • 3
    $100000!=3141510000099999!/31415≒pi*k$. So $sin100000!≒0.000$. – Takahiro Waki Jul 05 '16 at 09:14
  • 1
    Applying the taylor-series won't work. $1)$ The alternating signs lead to numerical instability $2)$ we need too many terms. It seems that we need a program with high precision-calculation and cannot compute , for example $sin(10^{10^{10^{10}}})$ – Peter Jul 05 '16 at 09:21
  • 1
    $(10^5)!$ has $456,574$ digits. Do you need the sine for even larger numbers ? You need less tha a million digits for $sin((10^5)!)$ – Peter Jul 05 '16 at 09:27
  • 1
    @TakahiroWaki This approximation is far too inaccurate to calculate $sin((10^5)!)$ – Peter Jul 05 '16 at 09:42
  • For those who are interested : $(10^5)!\ modulo\ 2\pi=4.33131820813936408660\ \cdot\ \cdot\ \cdot$ – Peter Jul 05 '16 at 10:10
  • So, the sine is $-0.9282669319268813746181175648\cdot\cdot\ \cdot$, coinciding with Mario's result. – Peter Jul 05 '16 at 10:15
  • @Peter : how do you know this is a good approximation ? (how do you compute $n! \bmod 2 \pi$ in general ?) – reuns Jul 06 '16 at 00:29
  • 1
    @TakahiroWaki 31415 is not nearly enough precision considering you're multiplying it by such an astronomical number... the "small" error of $.000096\ldots$ becomes amplified to way beyond the size of $\pi$. – MT_ Jul 06 '16 at 00:32
  • @user1952009 I used PARI/GP with a precision of a million digits. – Peter Jul 06 '16 at 10:04
  • @Peter : my question was : how do you know a million digit is enough for computing it with a 20 digit precision ? (and in general what is the needed precision for $\pi$ when computing $n ! \bmod 2 \pi$ with a 20 digit precision ?) – reuns Jul 06 '16 at 17:55
  • $(10^5)!$ has less than $500,000$ digits. (Since it is an integer, PARI can calculate it exactly). If we know $\pi$ with a precision of a million digits and divide $(10^5)!$ by $2\pi$, the number is a bit smaller than $(10^5)!$, so many of the digits of the fraction will be correct (It is nearly impossible that hundred thousands of nines or zeros follow). If I subtract $2\pi$ times this number from $(10^5)!$, I get the remainder and the precision will be almost the same. I see no reason why the final result should not be good to $20$ digits (In fact several hundred thousand digits) – Peter Jul 07 '16 at 22:02
  • The real challenge is to find a method to calculate $\sin(x)$ for monster numbers like $x=10^{10^{10^{10}}}$ or at least for $x=10^{10^{100}}$. Mario's answer could be a promising approach. – Peter Jul 07 '16 at 22:04
  • Problems with on chip value for pi: Intel Underestimates Error Bounds by 1.3 quintillion Intel https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ – dantopa Mar 17 '17 at 01:07

1 Answers1

15

I believe you may be able to calculate this without obscene numbers of digits of $\pi$ if you take advantage of the fact that these are factorials. To simplify the algebra, we can calculate $a_n=e^{i(n!)}$ instead (you want the imaginary part). Then $$a_{n+1}=e^{i(n!)(n+1)}=a_n^{n+1},$$ and it's perfectly reasonable to calculate $a_{100000}$ recursively with a high-precision library.

The downside is that to start the recursion you need a very good approximation of $e^i$, and I don't know if the error dependence works out any differently than in the $\pmod{2\pi}$ approach.

But to answer your actual question, Mathematica doesn't even break a sweat with the mere million digits needed for this:

> Block[{$MaxExtraPrecision = 1000000}, N[Sin[100000!], 10]]

-0.9282669319

takes about 15 ms on my computer.

For calculating the sine or cosine of a large arbitrary precision real number $x$, the gains of this method (which are tuned for $\sin n$ for integer $n$) are mostly lost, so I would recommend your original idea of reducing the argument $\bmod 2\pi$. As has been noted, the main bottleneck is a high-precision estimation of $\pi$. Your answer will be useless unless you can at least calculate $\frac{x}{\pi}$ to within $1$ (otherwise you may as well answer "somewhere between $-1$ and $1$"), so you need at least $\log_2(x/\pi^2)$ bits of precision for $\pi$. With $x\approx100000!$, that's about $1516701$ bits or $456572$ digits. Add to this the number $a$ of bits of precision you want in the result, so about $1516734$ digits of $\pi$ to calculate $33$ bits ($\approx 10$ digits) of $\sin x$ in the range $x\approx 100000!$.

Once you have an integer $n$ such that $y=2\pi n$ is close to $x$ (ideally $|x-2\pi n|\le2\pi$, it doesn't have to be perfectly rounded), calculate $\pi$ to precision $a+\log_2(n)$, so that $y$ is known to precision $a$, and then $x-y$ is precision $a$ and $\sin x=\sin (x-y)$ can be calculated to precision $a$ as well.

Mario Carneiro
  • 28,221
  • If I am right, the relative error on a power is the relative error on the base times that power. So you will need to evaluate $e^i$ with $n!$ decimals. –  Jul 05 '16 at 09:37
  • 3
    @YvesDaoust Actually, you only need error within $\frac1{n!}$ of the target error, which means about $\log(n!)$ decimals. $n!$ decimals would truly be obscene. – Mario Carneiro Jul 05 '16 at 09:40
  • you are right, my bad. –  Jul 05 '16 at 09:46
  • I don't actually compute the sine of a factorial but of a number as big as 100000! - but anyway: how is Mathematica actually capable of computing the sine that fast?! – Oliver Borchert Jul 05 '16 at 18:04
  • @OliverBorchert Actually, nothing about this method really depends on it being $n!$. I'm just taking $e^i$ to a very large integer power, which you could do just as well with other numbers by exponentiation by squaring. The important part is that to calculate $e^{in}$ you need about $a+\log_2(n)$ bits of precision on $e^i$ to get $a$ bits at the end. – Mario Carneiro Jul 05 '16 at 20:11
  • @OliverBorchert If you want to calculate the sine or cosine of a large real number, I would suggest breaking it into the form $n+x$ where $0\le x<1$ and then calculate $(e^i)^ne^{ix}$ using an $a+\log_2(n)$ bit approximation of $e^i$ and an $a$ bit approximation of $e^{ix}$ (which you can use standard sine and cosine math functions for). – Mario Carneiro Jul 05 '16 at 20:23
  • @MarioCarneiro your last comment could be extremely helpful (I'll try that out). But what exactly is a? (And what is n supposed to be? log_2(100000!) is still more than a million..) Would it be possible to add this approach (and some code) - I actually don't know a lot about complex numbers - to your answer? – Oliver Borchert Jul 05 '16 at 20:37
  • @OliverBorchert $a$ is the desired precision, for example you said you wanted 10 digits of precision so $a$ is about $33$ bits. I'm afraid you're on the wrong site for specific code examples, but you can calculate the product of two complex numbers via $(a+ib)(c+id)=(ac-bd)+i(ad+bc)$ - just store two real numbers of high precision and do the multiplications like that. (Many high-precision libraries already come with support for complex numbers; see if yours does.) If you have a high-precision sine and cosine, you can calculate $e^i=\cos 1+i\sin 1$. – Mario Carneiro Jul 06 '16 at 00:08