sin(x) & cos(x) for x = 10²² in radians

[SlideRule]

An interesting table comparing 'computed' vs. 'mathematical' value of sin(x) & cos(x) for x = 10²² {radians} on various current computing systems … Note that many computing systems return 0.0 or an error message when x is large. The worst behavior that the author has come across was on the Casio fx-8100. It produces a fatal error computing sin(x) when x > 26.
[attachment=12607]
BEST!
SlideRule

source: Argument Reduction for Huge arguments: Good to the Last Bit, SunPro, K. C. Ng, July 13, 1992

[dm319]

Interesting!

Any ideas why this is is difficult to answer?

Presumably relates to finding where 10^22 falls in the way of full turns plus radians?

On plus42, I find that 10^22 radians has 1591'549'430'918'953'357'688 full rotations, leaving behind 5.26300791462 radians.

Sin(5.26..) = -0.8522 etc
Cos(5.26..) = 0.5232

which gives similar answers to Sin(1e22) and cos(1e22). They are different by the order of 10^-13.

[dm319]

Out of curiosity, I thought I'd run this on a HP-45 emulator (HP-1973). I don't have any retro HPs to try it on in real life, but I think this is the earliest with Radian support.

HP-45:
Got sin(1e22) as -0.944... and cos(1e22) as -0.2953...

[lrdheat]

My CASIO’s (991CW and CLASSPAD 400) can do sin(10^8), but not sin(10^9). My TI TI-30X Pro MathPrint can do sin(10^9), but not sin(10^10). My HP 42, WP 34S, HP 10b11+ all get the correct answer for sin(10^22). The CLASSPAD surprised me with it’s fail.

[John Keith]

(10-02-2023 04:24 PM)dm319 Wrote:  HP-45:
Got sin(1e22) as -0.944... and cos(1e22) as -0.2953...

That seems off by a bit. On the HP 50 I get -.852200849762, where the last digit should be a 7. I would guess that any Saturn-based HP would get the same result. It would be interesting to see how various HP models improved over time.

[lrdheat]

The TI Nspire CX 11 (non CAS) can do sin(10^12), but not sin(10^13).

[KeithB]

(10-02-2023 05:13 PM)John Keith Wrote:  

(10-02-2023 04:24 PM)dm319 Wrote:  HP-45:
Got sin(1e22) as -0.944... and cos(1e22) as -0.2953...

That seems off by a bit. On the HP 50 I get -.852200849762, where the last digit should be a 7. I would guess that any Saturn-based HP would get the same result. It would be interesting to see how various HP models improved over time.

HP 71 is the same as what you got. Of course, it is a Saturn CPU. 8^)

[dm319]

Tried a couple I have at home:

Casio FX-6300G : Ma ERROR

Casio FX-115W : Ma ERROR

Ah well.

[brouhaha]

quote='dm319' pid='177906' dateline='1696263025']
Any ideas why this is is difficult to answer?
[/quote]

When asked to evaluate a trig function with an argument of 10^22, the first question is whether you mean EXACTLY 10^22, or (on a calculator) 10^22 within the limited floating point resolution of the calculator. In calculators with 10-digit mantissa precision, a display of 1.000000000E22 could be anywhere from 0.999999995x10^22 to 1.0000000004999...x10^22. Obviously the sine or cosine could then be any real value in [-1.0, 1.0]. That's the motivation for giving an error. Even if you can do perfectly accurate range reduction on an EXACT 10^22, a normal calculator has no way to know that the argument is exact.

If we did somehow know that the argument is exact, to do correct argument range reduction, we need an approximation of pi good to more than 22 digits for this case, and more than 99 digits for an exponent range to 99.

In comparison, if we ask a 10-digit calculator for sin(999999999.0), it still doesn't know that the argument is exact, but the maximum argument error is 0.05 which will result in inaccuracy of the result, but because sin and cos are well-behaved, we can still get three digits of accuracy in the result, vs zero digits for 10^10 or more.

In general, for any real world problem, if you need sin(10^22), you're doing something very wrong.

For a good background on issues with argument range reduction, see Elementary Functions: Algorithms and Implementation by Muller.

[dm319]

And on 'bc':

s(1^22) : .84147098480789650665

c(1^22) : .54030230586813971740

on 'R':

sin(1^22) : 0.841471

cos(1^22) : 0.5403023

on Python:

math.sin(1^22) : -0.8462204041751706

math.cos(1^22) : -0.5328330203333975

on Libreoffice Calc:

=sin(1^22) : 0.841470984807896

=cos(1^22) : 0.54030230586814

[dm319]

I keep finding more!

Galculator:

sin 1e22 : -0.852200849767
cos 1e22 : 0.5232147854

Quite impressed with this, it must mean that galculator implements a different maths library to R/python etc., I guess in BCD because it can do 0.3 - 0.2 without trouble. Probably Intel BCD I'm guessing, or else much higher precision floating point binary.

Android Calculator:

cos(1^22) : 0.5403023058681397174009...

sin(1^22) : 0.8414709848078078965066...

I wonder where it's going wrong for some of these... the numbers aren't far out, but the sign is often wrong.

[Maximilian Hohmann]

Hello,

I entered it in Google and it came up with the correct answer. I have no idea what Google uses for their calculations (maybe Wolfram Alpha?) but obviously they made a good choice.
Anyway, like brouhaha already wrote, if you ever enter such a problem into your calculator you should consider another year at college maybe…

Regards
Max

[Paul Dale]

(10-02-2023 04:10 PM)dm319 Wrote:  Any ideas why this is is difficult to answer?

It's the range reduction as you surmised. You need enough digits of pi to get a "good" result. In this case 22 + the number of digits you want the result to plus a few.

The WP 34S handles this case and worse just fine.

Pauli

[dm319]

(10-02-2023 06:34 PM)brouhaha Wrote:  In general, for any real world problem, if you need sin(10^22), you're doing something very wrong.

Thanks for the explanation. Yes, this is very true, but it is also quite fun pushing these calculators to their limits and seeing what fails.

(10-02-2023 09:57 PM)Paul Dale Wrote:  It's the range reduction as you surmised. You need enough digits of pi to get a "good" result. In this case 22 + the number of digits you want the result to plus a few.

Makes sense. The trig functions are such a rabbit hole. Still curious as to what exactly goes wrong with R/python/etc..

[Albert Chan]

(10-02-2023 04:10 PM)dm319 Wrote:  Interesting!

Any ideas why this is is difficult to answer?

Presumably relates to finding where 10^22 falls in the way of full turns plus radians?

Problem is cost of storing the constants, and the extra code to make it work.
Benefit may not be forthcoming, since float may represent a range (error ≤ 1/2 ulp)

To reduce costly arbitrary precision math, constants may be precomputed.
Lots of digits of pi itself doesn't work, if we can't efficiently do mod (2 pi)

This is how intel-decimal library does this (used in Free42/Plus42), turning mod (2 pi) to dot-products.

static BID_UINT384 bid_decimal128_moduli Wrote:// Values of (10^a / (2 pi)) mod 1 for -35 <= a <= 6111
// Each one is a 384-bit binary fraction. This may be a bit too much!
// My rough guideline is a bit more than 3x the working precision
...

384 bits = 48 bytes
Array length = 6111 - (-35) + 1 = 6147      → table size of almost 300K !

I recently downloaded Python 3.8.9 for my Win7 laptop. It had fixed the issue.

p3> from cmath import *
p3> exp(1e22j) # = cis(1e22)
(0.523214785395139-0.8522008497671888j)
p3> log(_).imag # reduced angle
-1.020177392559087

Older python work too (below is Python 2.6.6), if we use gmpy2

p2> from gmpy2 import *
p2> exp(1e22j)
mpc('0.52321478539513899-0.85220084976718879j')
p2> log(_).imag
mpfr('-1.0201773925590869')

Chez Scheme 9.5 fixed this too (newer version here)

scheme> (exp 0+1e22i)
0.523214785395139-0.8522008497671888i
scheme> (imag-part (log last))
-1.020177392559087

Another way to get arg(cis(θ)) is 2*atan(tan(θ/2))
RHS form is cheaper to compute.

Free42/Plus42:
Note: angle reduction code honors DEG/RAD/GRAD mode. I am not sure to count this as feature or bug.

Code:

00 { 13-Byte Prgm }
01▸LBL "A↓"
02 2
03 ÷
04 TAN
05 ATAN
06 ENTER
07 +
08 END

XEQ "RAD" 1e22 XEQ "A↓"      → -1.020177392559086973318201985281164

A convergent of pi (53 bits numerator) = 5706674932067741/1816491048114374

5706674932067741 XEQ "A↓" → 4.237546464512562184164292621941626e-16

[John Keith]

(10-02-2023 10:59 PM)Albert Chan Wrote:  scheme> (exp 0+1e22i)
0.523214785395139-0.8522008497671888i

Again, HP 50 approximate mode: (0, 1.E22) EXP
( .523214785403, -.852200849762 )
Real part looks bad but only off by 8 ULP's.

On the HP11C (no complex, so no EXP).
1E22 SIN -> -.984807753
1E22 COS -> .1736481777

There seems to have been a major improvement between the Nut era and the Saturn era.

[Albert Chan]

(10-02-2023 11:49 PM)John Keith Wrote:  Again, HP 50 approximate mode: (0, 1.E22) EXP
( .523214785403, -.852200849762 )
Real part looks bad but only off by 8 ULP's.

Actually, reduced angle is off by only 1 ULP (more precisely, 0.9087 ULP)
It look bad because reduced angle start with 1.0, making relative error big.
(this is more a problem with decimal vs binary machine, not angle reduction code)

I wonder how big is lookup table to get this good ...

>x = -1.020177392559087 ! = round(arg(cis(1e22)), -16)
>x
-1.02017739256
>exp((0,x))
 (.523214785394,-.852200849768)
>exp((0,neighbor(x,0)))
 (.523214785403,-.852200849762)

[lrdheat]

Curious concerning the problem of a calculator with 10 display digits (3 internally?) such as the HP 42S and HP 10b11+…how do they come up with such and accurate answer for sin(10^22)? sin(10^22)=-.852200849762 on the HP 42S, -.852200850 on the HP 10b11+ ?

[AndiGer]

(10-02-2023 04:24 PM)dm319 Wrote:  Out of curiosity, I thought I'd run this on a HP-45 emulator (HP-1973). I don't have any retro HPs to try it on in real life, but I think this is the earliest with Radian support.

HP-45:
Got sin(1e22) as -0.944... and cos(1e22) as -0.2953...

My 45 gives the same on sin but -0.329593335 for cos
Serial 1349S59…

[Joe Horn]

DM42 and DM32 are impressive: they get sin(1E22) correct to all 34 digits, and cos(1E22) off by only one ULP.