Previous: Bit Counting Example, Up: Example Definitions

##### 18.5.5.2 The Sine Function

A somewhat limited sine function could be defined as follows, using the well-known Taylor series expansion for ‘`sin(x)`’:

         (defmath mysin ((float (anglep x)))
           (interactive 1 "mysn")
           (setq x (to-radians x))    ; Convert from current angular mode.
           (let ((sum x)              ; Initial term of Taylor expansion of sin.
                 newsum
                 (nfact 1)            ; "nfact" equals "n" factorial at all times.
                 (xnegsqr :"-(x^2)")) ; "xnegsqr" equals -x^2.
             (for ((n 3 100 2))       ; Upper limit of 100 is a good precaution.
               (working "mysin" sum)  ; Display "Working" message, if enabled.
               (setq nfact (* nfact (1- n) n)
                     x (* x xnegsqr)
                     newsum (+ sum (/ x nfact)))
               (if (~= newsum sum)    ; If newsum is "nearly equal to" sum,
                   (break))           ;  then we are done.
               (setq sum newsum))
             sum))

The actual `sin` function in Calc works by first reducing the problem to a sine or cosine of a nonnegative number less than ‘`pi/4`’. This ensures that the Taylor series will converge quickly. Also, the calculation is carried out with two extra digits of precision to guard against cumulative round-off in ‘`sum`’. Finally, complex arguments are allowed and handled by a separate algorithm.

         (defmath mysin ((float (scalarp x)))
           (interactive 1 "mysn")
           (setq x (to-radians x))    ; Convert from current angular mode.
           (with-extra-prec 2         ; Evaluate with extra precision.
             (cond ((complexp x)
                    (mysin-complex x))
                   ((< x 0)
                    (- (mysin-raw (- x)))    ; Always call mysin-raw with x >= 0.
                   (t (mysin-raw x))))))

         (defmath mysin-raw (x)
           (cond ((>= x 7)
                  (mysin-raw (% x (two-pi))))     ; Now x < 7.
                 ((> x (pi-over-2))
                  (- (mysin-raw (- x (pi)))))     ; Now -pi/2 <= x <= pi/2.
                 ((> x (pi-over-4))
                  (mycos-raw (- x (pi-over-2))))  ; Now -pi/2 <= x <= pi/4.
                 ((< x (- (pi-over-4)))
                  (- (mycos-raw (+ x (pi-over-2)))))  ; Now -pi/4 <= x <= pi/4,
                 (t (mysin-series x))))           ; so the series will be efficient.

where `mysin-complex` is an appropriate function to handle complex numbers, `mysin-series` is the routine to compute the sine Taylor series as before, and `mycos-raw` is a function analogous to `mysin-raw` for cosines.

The strategy is to ensure that ‘`x`’ is nonnegative before calling `mysin-raw`. This function then recursively reduces its argument to a suitable range, namely, plus-or-minus ‘`pi/4`’. Note that each test, and particularly the first comparison against 7, is designed so that small roundoff errors cannot produce an infinite loop. (Suppose we compared with ‘`(two-pi)`’ instead; if due to roundoff problems the modulo operator ever returned ‘`(two-pi)`’ exactly, an infinite recursion could result!) We use modulo only for arguments that will clearly get reduced, knowing that the next rule will catch any reductions that this rule misses.

If a program is being written for general use, it is important to code it carefully as shown in this second example. For quick-and-dirty programs, when you know that your own use of the sine function will never encounter a large argument, a simpler program like the first one shown is fine.