Sine Table
Challenge
🙜
July 24, 2021
Approach:
• Use trig identities:
sin(a+1) = sin(a)*cos(1) + cos(a)*sin(1)
cos(a+1) = cos(a)*cos(1) - sin(a)*sin(1)
• Compute 180 degrees of the table incrementally
• Negate for the second half
• Use single 64-bit cells to compute sine
• Use double 128-bit cells to check
0.0000000000000000000000000000000000000 2constant ZERO
0.5877852522924731291687059546390727686 2constant SIN36
0.9510565162951535721164393333793821434 2constant SIN72
SIN36 DNEGATE 2constant -SIN36
SIN72 DNEGATE 2constant -SIN72
create ref-sines
ZERO 2, SIN36 2, SIN72 2, SIN72 2, SIN36 2,
ZERO 2, -SIN36 2, -SIN72 2, -SIN72 2, -SIN36 2, ZERO 2,
: ref-sin ( n -- d ) 2* cells ref-sines + 2@ ;
1.0000000000000000000000000000000000000 2constant ONE
3.1415926535897932384626433832795028842 2constant dPI
: deg>rad ( n -- d ) dPI rot 1 * 180 m*/ ;
1000000000000000000 constant sONE
0017452406437283513 constant SIN1
0999847695156391239 constant COS1
: SIN1* ( n -- n ) SIN1 sONE */ ;
: COS1* ( n -- n ) COS1 sONE */ ;
( sin[a+1] = sin[a]*cos[1]+cos[a]*sin[1] )
( cos[a+1] = cos[a]*cos[1]-sin[a]*sin[1] )
: sines-step ( sin cos -- sin+1 cos+1 )
2dup COS1* swap SIN1* - >r SIN1* swap COS1* + r> ;
: fill-sines 0 sONE 179 for over , sines-step next ;
create sines fill-sines
: sin ( n -- d )
180 /mod swap cells sines + @
swap 2 mod if negate then ONE rot sONE m*/ ;
: deg. ( n -- ) s>d <# # # #s #> type ;
: un. ( d -- ) <# 36 for # next [char] . hold #s #> type ;
: n. ( d -- ) 2dup d0< if [char] - emit dnegate else space then un. ;
: sinerow. ( n -- )
dup 36 * deg. ." ~ "
dup 36 * deg>rad n. ." : "
dup ref-sin n. ." "
dup ref-sin rot 36 * sin d- dabs n. cr
dup 49 spaces 36 * sin n. cr ;
: dashes ( n -- ) 1- for [char] - emit next ;
: center ( a n w -- ) over - 2 /mod dup spaces + >r type r> spaces ;
: gap 4 spaces ;
: column ( a n ) gap 39 center ; : dcolumn gap 39 dashes ;
: split 3 dashes dcolumn dcolumn dcolumn cr ;
: sines. ." DEG" s" RAD" column s" REF-SIN vs COMP-SIN" column
s" DIFF" column cr split
0 10 for dup sinerow. 1+ next drop ;
sines. bye
DEG RAD REF-SIN vs COMP-SIN DIFF
--- --------------------------------------- --------------------------------------- ---------------------------------------
000 ~ 0.0000000000000000000000000000000000000 : 0.0000000000000000000000000000000000000 0.0000000000000000000000000000000000000
0.0000000000000000000000000000000000000
036 ~ 0.6283185307179586476925286766559005768 : 0.5877852522924731291687059546390727686 0.0000000000000000331687059546390727686
0.5877852522924730960000000000000000000
072 ~ 1.2566370614359172953850573533118011536 : 0.9510565162951535721164393333793821434 0.0000000000000000681164393333793821434
0.9510565162951535040000000000000000000
108 ~ 1.8849555921538759430775860299677017305 : 0.9510565162951535721164393333793821434 0.0000000000000000881164393333793821434
0.9510565162951534840000000000000000000
144 ~ 2.5132741228718345907701147066236023073 : 0.5877852522924731291687059546390727686 0.0000000000000000791687059546390727686
0.5877852522924730500000000000000000000
180 ~ 3.1415926535897932384626433832795028842 : 0.0000000000000000000000000000000000000 0.0000000000000000000000000000000000000
0.0000000000000000000000000000000000000
216 ~ 3.7699111843077518861551720599354034610 : -0.5877852522924731291687059546390727686 0.0000000000000000331687059546390727686
-0.5877852522924730960000000000000000000
252 ~ 4.3982297150257105338477007365913040378 : -0.9510565162951535721164393333793821434 0.0000000000000000681164393333793821434
-0.9510565162951535040000000000000000000
288 ~ 5.0265482457436691815402294132472046147 : -0.9510565162951535721164393333793821434 0.0000000000000000881164393333793821434
-0.9510565162951534840000000000000000000
324 ~ 5.6548667764616278292327580899031051915 : -0.5877852522924731291687059546390727686 0.0000000000000000791687059546390727686
-0.5877852522924730500000000000000000000
360 ~ 6.2831853071795864769252867665590057684 : 0.0000000000000000000000000000000000000 0.0000000000000000000000000000000000000
0.0000000000000000000000000000000000000