calc: Sine Example
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.
Info Catalog
calc: Bit Counting Example
calc: Example Definitions