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

[Paul Dale]

(10-03-2023 04:10 AM)Joe Horn Wrote:  DM42 and DM32 are impressive: they get sin(1E22) correct to all 34 digits, and cos(1E22) off by only one ULP.

Intel's sloppy maths library at work ?


Pauli

[rprosperi]

(10-03-2023 10:37 AM)Paul Dale Wrote:  

(10-03-2023 04:10 AM)Joe Horn Wrote:  DM42 and DM32 are impressive: they get sin(1E22) correct to all 34 digits, and cos(1E22) off by only one ULP.

Intel's sloppy maths library at work ?


Pauli

Damn those guys!! It's a good thing it's free!

[Ajaja]

(10-02-2023 06:58 PM)dm319 Wrote:  And on 'bc':
s(1^22) : .84147098480789650665
c(1^22) : .54030230586813971740

on 'R':
sin(1^22) : 0.841471
cos(1^22) : 0.5403023

on Libreoffice Calc:
=sin(1^22) : 0.841470984807896
=cos(1^22) : 0.54030230586814

Android Calculator:
cos(1^22) : 0.5403023058681397174009...
sin(1^22) : 0.8414709848078078965066...

It's just sin(1) and cos(1).

[dm319]

(10-03-2023 12:40 PM)Ajaja Wrote:  It's just sin(1) and cos(1).

So it is! Good spot. Now I'm trying to work out what's going on here. I'm guessing something to do with a limited precision pi, maybe multiplied up a lot or some sort of modulo calculatoin, maybe taking the integer part of something that doesn't have a fractional part... Hmmm.

[Joe Horn]

(10-03-2023 12:40 PM)Ajaja Wrote:  

(10-02-2023 06:58 PM)dm319 Wrote:  And on 'bc':
s(1^22) : .84147098480789650665
c(1^22) : .54030230586813971740

on 'R':
sin(1^22) : 0.841471
cos(1^22) : 0.5403023

on Libreoffice Calc:
=sin(1^22) : 0.841470984807896
=cos(1^22) : 0.54030230586814

Android Calculator:
cos(1^22) : 0.5403023058681397174009...
sin(1^22) : 0.8414709848078078965066...

It's just sin(1) and cos(1).

Maybe because 1^22 = 1?

[dm319]

(10-03-2023 02:13 PM)Joe Horn Wrote:  Maybe because 1^22 = 1?

<facepalm>! Serves me right for trying such a silly calculation on serious bits of software! I think it was because I was inputting 1e22. Ok I'll need to run those again to at least set the record straight.

[Jlouis]

Gaxio: +/- 0,6789876543215687654328

With double precision on.

Precise result, if it's not raining.

[Jlouis]

(10-03-2023 02:13 PM)Joe Horn Wrote:  

(10-03-2023 12:40 PM)Ajaja Wrote:  It's just sin(1) and cos(1).

Maybe because 1^22 = 1?

[dm319]

Ok lets try again not muddling up e for ^!

R:
Sin(10^22) : -0.8522008
Cos(10^22) : 0.5232148

Python 3.10.12:
math.sin(10^22) : 0.27090578830786904
math.cos(10^22) : -0.9626058663135666

Julia 1.9.3:
sin(10^22) : -0.7138413366055608
cos(10^22) : 0.7003074654415633

Libreoffice Calc:
VALUE! and VALUE!

Android Calculator:
sin(10^22) : -0.8522008497671888017727058937530...
cos(10^22) : 0.523214785395138945497594473384709...

BC:
s(10^22) : -.852200849767188801772705893753029368261762150410043656256509326025\
91031199209620153543628018037908962775444734135871466289192960541752\
16043583440915435596194706342456081425546376901913864520332009761575\
60796633380693540220007456689428988629977079633384981377028139856825\
53803450531758922619115674748231842464566971379516418439401013378169\
06152024705255524515959251984144296572956924819102578669162060547699\
99938433149951023354338851794739740671689887537040538589009732269418\
65066978450714962805064022980941788591903233291613477370182803451002\
71881558355654347085787733978013640483250684988058946829605911021898\
46890926169796552029258387512882452076937215435193053898062745371324\
36650823819125824317591259681988156401020352128082469004610666201995\
59828390637847945525492738363284282602995651681901027000201538098623\
52143653195089980176678047361861374760966870159859856968855103691674\
50318018031949140264007794360143089385888218113790071739783432417512\
91760137987646938423978141058291541637613409005287

c(10^22) : .5232147853951389454975944733847094921409199724393879535272113921042\
98247376710623283422632663065703715385969475522571594595880745174680\
99321764045804940209665373863421408519470724298211356174350312925233\
98021007499038858472332797896752075063671442470524436453359121516556\
92059065599678840901524160704752550092663164693824519608203841467373\
80803528574364758776522462732448745942524057325543154555034884711181\
47501402012988472497468939637543401786831366568442456574644912702060\
79372629404872529281963153892600854826686134815968789505080864866007\
14174348461193099621694555233629201922602257999209709395680292160284\
11734789771712082866554711231837258300506268315731006465207372434666\
59195270736258730243593424153868297088887975644725197732877520987861\
92993200753697676582734988415314876979196271838701902857681809783167\
72286687053782571235210822799722647879031199806540161655728783633866\
55272107859476026100879998571094124337377148333098768674226179154655\
1774338870167170459282266282375581076913475851709

Ok, I'm kinda surprised that R, Python and Julia all have different answers. I wonder if android calculator has some sort of arbitrary precision thing going on. Of course the oldest program here has to be the one that can show off the most. I've no idea if it is at all accurate though.

[Ajaja]

(10-03-2023 06:29 PM)dm319 Wrote:  Python 3.10.12:
math.sin(10^22) : 0.27090578830786904
math.cos(10^22) : -0.9626058663135666

In Python operator ^ means XOR. You should use here 10**22 or just 1e22.

[Nigel (UK)]

(10-03-2023 06:29 PM)dm319 Wrote:  Ok, I'm kinda surprised that R, Python and Julia all have different answers. I wonder if android calculator has some sort of arbitrary precision thing going on. Of course the oldest program here has to be the one that can show off the most. I've no idea if it is at all accurate though.

The Android calculator literally has precision limited only by memory. Just swipe left on the displayed result to get as many significant figures as you like. All displayed digits guaranteed correct! You can read a little about it here and more in the articles linked to from there. I'm sure it's been discussed on this forum previously too.

Nigel (UK)

[dm319]

(10-03-2023 07:56 PM)Ajaja Wrote:  In Python operator ^ means XOR. You should use here 10**22 or just 1e22.

Darn, I really need to stop now. Someone more competent can take over!

[dm319]

(10-02-2023 10:59 PM)Albert Chan Wrote:  turning mod (2 pi) to dot-products.

This is fascinating and makes sense. I guess this isn't an issue for calculating in degrees. Do you think the hp-20s uses a similar look up table?

[Albert Chan]

(10-03-2023 07:56 PM)Ajaja Wrote:  In Python operator ^ means XOR. You should use here 10**22 or just 1e22.

10^22 = 0x1.0f0cf064dd592p+73

It is the biggest pow-of-10 that can be exactly store in IEEE double. (10^23 need 54 bits)
Perhaps this is the reason for picking it, to compare both decimal and binary machines.

(10-05-2023 01:33 PM)dm319 Wrote:  Do you think the hp-20s uses a similar look up table?

Very likely.

Of course, there are ways to do this without table.
We don't need pi ... we just keep throwing precisions to it

lua> qd = require'qd' -- quad double, precision = 4*53 bits
lua> t = 1e22 / qd.pow(5,40) -- about 1e-6
lua> tt = t*t
lua> s = t*(1 - tt/6*(1 - tt/20*(1 - tt/42*(1 - tt/72*(1 - tt/110)))))
lua> s -- ≈ sin(10^22 / 5^40)
1.0995116277757784620007025274123330676478201219099885281501104800e-06

Now, map 40 times, sin(x) --> sin(5x)

lua> m5 = fn'x: x*(5-x*x*(4+16*(1-x)*(1+x)))'
lua> for i=1,40 do s = m5(s) end
lua> s -- ≈ sin(10^22)
-8.5220084976718880177270589375302936826176863529279677310847219854e-01

[Ajaja]

(10-05-2023 03:08 PM)Albert Chan Wrote:  10^22 = 0x1.0f0cf064dd592p+73
It is the biggest pow-of-10 that can be exactly store in IEEE double. (10^23 need 54 bits)

Yes, and this whole comparison of sin(10^22) is a bit misleading.
In Python for example:

Code:

import math
a = 10**22
b = 10**22 - 1 
print(a)
print(b)
print(math.sin(a))
print(math.sin(b))

Output:

Code:

10000000000000000000000
9999999999999999999999
-0.8522008497671888
-0.8522008497671888

The math library converts both integers to double and gets exactly the same IEEE number. The mathematically correct result for 10^22 here it's just like "even a broken clock is right twice a day"

[Albert Chan]

(10-03-2023 12:56 AM)Albert Chan Wrote:  I wonder how big is [HP71B] lookup table to get this good ...

I should just try it! No table, just many digits of pi.

>sin(1e22)
-.852200849762
>sin(1e33)
-5.11830428837E-2

Using rpn.exe, HP71B stored 32 digits of pi (note, next digit of pi is 0)
Note, "pi =-31" means digits=31, but only for pi, global digits (default=36) unaffected.

C:\> rpn
pi ?
3.14159265358979323846264338327950288

1e22 pi =-31 2x % sin ?12
-0.852200850595
1e22 pi =-32 2x % sin ?12
-0.852200849762
1e22 sin ?12
-0.852200849767

1e33 pi =-31 2x % sin ?12
0.899003019303
1e33 pi =-32 2x % sin ?12
-0.0511830428837
1e33 sin ?12
0.762302394953

For comparison, Free42 Decimal "stored" thousands of digits of pi

(10-24-2021 12:46 PM)Albert Chan Wrote:  min(|cos(x)|) for integer x (also, for real x)

2344813655066356855719930664718056e1381 COS ABS → 1.030557387629248882465543827418861e-37

C:\>spigot --printf "%.34g" cos(2344813655066356855719930664718056e1381)
1.030557387629248882465543827418861e-37

[J-F Garnier]

(10-05-2023 09:20 PM)Albert Chan Wrote:  

(10-03-2023 12:56 AM)Albert Chan Wrote:  I wonder how big is [HP71B] lookup table to get this good ...

I should just try it! No table, just many digits of pi.

>sin(1e22)
-.852200849762
>sin(1e33)
-5.11830428837E-2

Using rpn.exe, HP71B stored 32 digits of pi (note, next digit of pi is 0)

Actually, the Saturn machines are using the PI/4 constant with 31 digits.

The main goal of the PI/4 double precision value was not to manage huge exact arguments, since this has no real-life interest and there will always be a limit unless storing 500 digits, but to guarantee 12-digit accuracy trig results for all "realistic" arguments, especially SIN(X) with X close to PI.

J-F

[KeithB]

"Actually, the Saturn machines are using the PI/4 constant with 31 digits."

How many significant digits does the "45" have? 8^)
(If this confuses you, look at the source in the link)

[Albert Chan]

(10-06-2023 08:40 AM)J-F Garnier Wrote:  Actually, the Saturn machines are using the PI/4 constant with 31 digits.

Thanks for the info.

I still think code somehow adjusted stored 31-digits pi/4 (error ≈ -1/4 ULP)
Otherwise, it cannot match HP71B SIN(x) signatures.

HP71B Wrote:>for i=21 to 30 @ sin(10^i) @ next i
-.667120177071
-.852200849762
 .701140639796
-.996472229026
-.744789854871
 .853455100159
 .718062857184
-.987830531786
 .999958451937
-9.10311968368E-2

32 digits of pi, it mached them all.

Quote:C:\> seq 21 30 | rpn pi =-32 2x st [ 10 $ pow dup rt / f rt x - sin ?12 k
-0.667120177071
-0.852200849762
 0.701140639796
-0.996472229026
-0.744789854871
 0.853455100159
 0.718062857184
-0.987830531786
 0.999958451937
-0.0910311968368

31 digits of pi/4 (= 32 digits of pi, error = -1 ULP), no match:

Quote:C:\> seq 21 30 | rpn pi 4/ =-31 8x st [ 10 $ pow dup rt / f rt x - sin ?12 k
-0.667120177095
-0.852200849929
 0.701140642065
-0.996472231697
-0.744789642463
 0.853456759079
 0.718085010496
-0.987780973718
 0.999982401941
-0.0592916087884

Note: "[ ..." means setup a pipe, i → round(sin(remainder(10^i, t ≈ 2*pi)), -12)

[Gilles]

NewRPL give perfect result in 32 digits précisions and beyond. For example
127 SETPREC
1E22 SIN
returns the 126 good decimals