dontrig.py: simplified and documented.

generate/dontrig.py

@@ -1,46 +1,48 @@
+'''dontrig.py
+Don Cross <[email protected]>
+https://github.com/cosinekitty/astronomy
+
+These are my own hand-rolled replacements for sine and cosine
+on platforms where I'm getting divergent results (order 1e-12).
+I want to find out who is correct and who is wrong!
+'''
 import math
 
 Debug = False
 
 def xsin(angle:float) -> float:
     angle = math.fmod(angle, 2*math.pi)
-    angleSquared = angle * angle
+    numerator = -(angle * angle)
     sum = 0.0
-    fact = 1.0
-    pow = angle
+    term = angle
     n = 1
     while n <= 99:
         # n = 1, 3, 5, ...
-        term = pow / fact
         prev = sum
         sum += term
         if Debug:
             print('xsin: n={:02d}  term={:24.16e}  sum={:20.16f}  fact={:g}  diff={:24.16e}'.format(n, term, sum, fact, sum-prev))
         if prev == sum:
             return sum
-        pow *= angleSquared
-        fact *= -((n+1) * (n+2))
+        term *= numerator / ((n+1) * (n+2))
         n += 2
     raise Exception('xsin({:0.16g}): failure to converge'.format(angle))
 
 
 def xcos(angle:float) -> float:
     angle = math.fmod(angle, 2*math.pi)
-    angleSquared = angle * angle
+    numerator = -(angle * angle)
     sum = 0.0
-    fact = 1.0
-    pow = 1.0
+    term = 1.0
     n = 0
     while n <= 99:
         # n = 0, 2, 4, ...
-        term = pow / fact
         prev = sum
         sum += term
         if Debug:
             print('xcos: n={:02d}  term={:24.16e}  sum={:20.16f}   fact={:g}'.format(n, term, sum, fact))
         if prev == sum:
             return sum
-        pow *= angleSquared
-        fact *= -((n+1) * (n+2))
+        term *= numerator / ((n+1) * (n+2))
         n += 2
     raise Exception('xcos({:0.16g}): failure to converge'.format(angle))