Calculating Sines

January 12, 2010

Taylor Series

Just as we did previously when calculating logarithms, we start our solution for calculating sines by defining epsilon as the desired accuracy of the final result; we also give a definition for π:

(define epsilon 1e-7)

(define pi 3.141592654)

Since each term of the Taylor series contains a factorial in its denominator, we’ll need a fast way to compute factorials. The following iterative procedure comes from section 1.2.1 of Structure and Interpretation of Computer Programs (SICP):

(define (factorial n) (fact-iter 1 1 n))

(define (fact-iter product counter max-count) (if (> counter max-count) product (fact-iter (* counter product) (+ counter 1) max-count)))

The following procedure will compute the nth term of the Taylor sine series given n and the angle in radians:

(define (term n radians) (* (/ (expt radians (+ (* 2 n) 1)) (factorial (+ (* 2 n) 1))) (expt -1 n)))

Since the sine function is periodic, computation time for large angles can be greatly improved by reducing the angle to the range [-π, π] before computing the Taylor series.

(define (reduce x) (- x (* (round (/ x (* 2 pi))) 2 pi)))

Finally, we can define a procedure that iteratively adds terms of the Taylor series to a sum until the desired precision is reached.

(define (good-enough? current next) (< (abs (- current next)) epsilon))

(define (sine-iter radians n current next) (if (good-enough? current next) next (sine-iter radians (+ n 1) next (+ next (term (+ n 1) radians)))))

(define (taylor-sine radians) (sine-iter (reduce radians) 0 0 (term 0 (reduce radians))))

The sine-iter procedure keeps track of both the current and the next sum of the terms of the Taylor series, only halting execution when the good-enough? procedure determines that the difference between the two is smaller than the epsilon value we defined earlier. The taylor-sine procedure simply reduced the input angle and calls the iterative procedure with initial values.

Here are a few sample calculations:

> (taylor-sine 1) 0.841470984648068 > (taylor-sine (/ pi 2)) 1.0000000006627803 > (taylor-sine 10) -0.5440211113737637

Recursive formula

One solution that uses the approximation sin x ≈ x, and the triple-angle formula $\sin x = 3 \sin\frac{x}{3} - 4 \sin^3\frac{x}{3}$ can be found in SICP section 1.2.3:

(define (cube x) (* x x x)) (define (p x) (- (* 3 x) (* 4 (cube x)))) (define (sine angle) (if (not (> (abs angle) 0.1)) angle (p (sine (/ angle 3.0)))))

The only thing we need to change is the criteria for considering an angle “sufficiently small.” If we change 0.1 in the comparison on the 4th line of the code above to the epsilon value we defined earlier, we gain the needed precision.

(define (cube x) (* x x x)) (define (p x) (- (* 3 x) (* 4 (cube x)))) (define (sine angle) (if (not (> (abs angle) epsilon)) angle (p (sine (/ angle 3.0)))))

Once again, here are a few sample calculations:

> (sine 1) 0.8414709848078971 > (sine (/ pi 2)) 1.0 > (sine 10) -0.5440211108893757

You can run the program at http://programmingpraxis.codepad.org/uyaebBRj.

19 Responses to “Calculating Sines”

programmingpraxis said
January 12, 2010 at 2:07 PM
Here is my version of taylor-sine, which keeps the numerator and denominator separately and calculates each by updating the previous value rather than recalculating each from scratch at each iteration.

(define (taylor-sine x) (let loop ((n x) (d 1) (k 1) (s x) (t -1)) (let* ((next-n (* n x x)) (next-d (* d (+ k 1) (+ k 2))) (next (* t (/ next-n next-d)))) (if (< (abs next) epsilon) (exact->inexact (+ s next)) (loop next-n next-d(+ k 2) (+ s next) (* t -1))))))

I also eschew the range-reduction optimization, as http://dotancohen.com/eng/taylor-sine.php shows little improvement from the optimization.

novatech said

January 12, 2010 at 4:32 PM

my solution using c

#include <stdio.h>
#include <math.h>
#define PI 3.141592654

long double factorial(int n)
 {
  return (n<=1) ? 1 : n*factorial(n-1);
 }
double tylor_sine(double x) {
int n=0,fac;
double result=0,term=0;
do {
        fac=(2*n+1);
        result=result+term;
        term=pow(-1,n)*(pow(x,fac)/factorial(fac));
        n++;
       }while(fabs(term)>=.00000001);
 return result;
}
double cube (double x) {
    return x*x*x;
}
double sine (double angle) {
    if (angle < 0.00000001) { return angle; } 
    return 3*(sine (angle/3.0)) - 4*cube(sine (angle/3.0));
}
int main() {
    printf("tylor-series sine(10): %f\n",tylor_sine(10.0));
    printf("triple-angle sine(10): %f\n",sine(10.0));
    printf("tylor-series sine(1): %f\n",tylor_sine(1.0));
    printf("triple-angle sine(1): %f\n",sine(1.0));
    printf("tylor-series sine(PI/2): %f\n",tylor_sine(PI/2));
    printf("triple-angle sine(PI/2): %f\n",sine(PI/2));
    return 0;
}

Bill the Lizard said
January 12, 2010 at 4:56 PM
programmingpraxis,
I’d say that the page you linked to gives a strong argument in favor of range reduction. It shows that when you reduce the range, only very few terms of the Taylor series are required to accurately approximate sine(x), even for large x. Your optimization of keeping the numerator and (particularly) the denominator, and updating each instead of recalculating on each iteration helps a lot. However, as x grows, even just updating the factorial in the denominator can take quite a bit of time. For example, without using range reduction it takes my machine about 25 seconds to compute (taylor-sine 1000). With range reduction I get the answer immediately, with no significant loss of accuracy.
Programming Praxis – Calculating Sines « Bonsai Code said
January 12, 2010 at 7:26 PM
[…] Praxis – Calculating Sines By Remco Niemeijer In today’s Programming Praxis exercise we have to implement two ways of calculating sines. Let’s get […]

Remco Niemeijer said

January 12, 2010 at 7:26 PM

My Haskell solution (see http://bonsaicode.wordpress.com/2010/01/12/programming-praxis-calculating-sines/ for a version with comments):

import Data.Fixed

taylorSin :: Double -> Double
taylorSin x = sum . useful $ zipWith (/) 
    (map (\k -> (mod' x (2*pi)) ** (2*k + 1) * (-1) ** k) [0..])
    (scanl (\a k -> a * k * (k - 1)) 1 [3,5..]) where
    useful ~(a:b:c) = a : if abs (a-b) > 1e-7 then useful (b:c) else []

recSin :: Double -> Double
recSin = f . flip mod' (2 * pi) where
    f x = if abs x < 1e-7 then x else
          let s = f (x / 3) in 3 * s - 4 * s**3

programmingpraxis said
January 12, 2010 at 9:50 PM
Bill the Lizard: You are correct. My mistake.
Matías Giovannini said
January 13, 2010 at 1:35 AM
My purely iterative solution, with range reduction, in 27 lines: http://codepad.org/HgaJbSbh
Manish Mathai said
January 13, 2010 at 5:00 AM
Just a nitpicking. Shouldn’t it be lim ( sin(x) / x ) = 1 when x -> 0
norsetto said
January 13, 2010 at 1:32 PM
My first attempt to programming in haskell:

sin' :: (Double -> Double) -> Double -> Double sin' f x = reduce1 f (x - fromIntegral(truncate(x / (2*pi)))*(2*pi))
reduce1 :: (Double -> Double) -> Double -> Double reduce1 f x | x Double) -> Double -> Double reduce2 f x | x <= pi/2 = f x | x <= pi = f (pi - x) | x Double taylor x = let (_,_,_,r) = foldl sumOdd (1,1,-1,0) [1..20] where sumOdd (n,d,s,p) el | even el = (n*x,d*fromIntegral el,-s,p-s*n/d) | otherwise = (n*x,d*fromIntegral el,s,p) in r
sinIter :: Double -> Double sinIter x | x < 1e-7 = x | otherwise = 3*sinIter(x/3) - 4*sinIter(x/3)**3

And is amazingly working:

*Main> sin’ taylor (pi*0.25)
0.7071067811865475
*Main> sin’ sinIter (pi*0.25)
0.7071067811865475
*Main> sin (pi*0.25)
0.7071067811865475

And it took me only 10 times more time that it would have taken me in C …
programmingpraxis said
January 13, 2010 at 1:49 PM
Manish: Fixed. My fault, not Bill’s.
Erik Schulz said
January 14, 2010 at 12:03 AM
A version written in F#.

let epsilon = 1e-7
let pi = 3.141592654

let sin (x:float) =
let rec series sum f d v k2 =
let next = f * d / float v
if abs(next) < epsilon then
sum + next
else
series (sum + next) (-f) (d * x * x) (v * (k2 + 1) * (k2 + 2)) (k2 + 2)
series 0.0 1.0 1.0 1 1

let rec sin' x =
if abs(x) float
val sin’ : float -> float

> sin 1.0;;
val it : float = 0.8414709846
> sin’ 1.0;;
val it : float = 0.8414709848
Erik Schulz said
January 14, 2010 at 12:06 AM
the comment section ate my code …
Erik Schulz said
January 14, 2010 at 12:07 AM
Here is a version written in F#.

http://pastie.org/777280
Norsetto's Journal said
January 15, 2010 at 4:05 PM
[…] dabbling with haskell in the recent days. My first semi-serious attempt was inspired by a prompt at programming praxis. The code I concocted is as follows (I guess it would ashame any serious haskell programmer, if you […]

Mike said

April 10, 2010 at 12:34 AM

Python:


def reduce_range( x ):
    sign,x = (1, x) if x >= 0 else (-1, -x)
    
    if x >= 2*pi: x -= int(x/2/pi) * pi
    if x > pi: x = x - 2*pi
    
    return x if sign>0 else -x


def sin_taylor( x, eps=1.0e-7 ):
    x = reduce_range( x )        
    x2 = x*x
    sinx = 0
    num, den, k, sign = ( float( x ), 1.0, 3.0, 1 )
    
    while True:
        delta = num / den
        sinx += sign * delta

        if delta < eps: break

        num, den, k, sign = ( num * x2, den*k*(k-1), k+2, -sign )

    return sinx

def sin_3angle( x, eps=1e-7 ):

    def _3angle( x ):
        if x < eps:
            return x
        else:
            f = _3angle( x/3.0 )
            return f * ( 3 - 4*f*f )

    return _3angle( reduce_range( x ) )

David said

March 28, 2011 at 3:15 AM

Version in Factor. Interesting that the Taylor series is minimally accurate to the given precision (1e-7), but the recursive version is much more precise than our minimum of 1e-7. Probably the case as for values that small sin x ~= x.

USING: kernel math math.constants math.libm locals
       arrays sequences lists lists.lazy ;
IN: sines

CONSTANT: epsilon 1e-7

: normalize-theta  ( x -- x )
    [let 1 :> sign!
      dup 0 <
        [ -1 *  -1 sign! ]  when

      pi 2 *  fmod

      dup pi >
        [ 2 pi * - ] when

      sign *
    ] ;


! create the lazy list of taylor terms (not so simple :)

: numerators  ( x -- lazy-list )
    dup 2array  [ first2
                  over *
                  over *
                  -1 *
                  2array
                ] lfrom-by  [ second ] lazy-map ;

: denominators  ( -- lazy-list )
    { 1 1 }  [ first2
               [ 1 + ] dip over *
               [ 1 + ] dip over *
               2array
             ] lfrom-by  [ second ] lazy-map ;

: taylor-terms ( x -- lazy-list )
    numerators denominators lzip [ first2 / ] lazy-map ;

: taylor-sine ( x -- x )
   >float normalize-theta taylor-terms  [ abs epsilon < ] luntil  0 [ + ] foldl ;


! sin x = 3 sin (x/3) - 4 sin^3 (x/3)

: sine  ( x -- x )
    dup abs epsilon >
      [ 3.0 / sine [ 3 * ] [ dup sq * 4 * ] bi - ]          
    when ;

Session:


( scratchpad ) pi 4 / sin .
0.7071067811865475
( scratchpad ) pi 4 / sine .
0.7071067811865475
( scratchpad ) pi 4 / taylor-sine .
0.7071067811796195
( scratchpad ) 1 sin .
0.8414709848078965
( scratchpad ) 1 sine .
0.8414709848078971
( scratchpad ) 1 taylor-sine .
0.841470984648068
( scratchpad ) 100 sin .
-0.5063656411097588
( scratchpad ) 100 sine .
-0.5063656411096491
( scratchpad ) 100 taylor-sine .
-0.5063656411334029
( scratchpad ) pi 1.5 * sin .
-1.0
( scratchpad ) pi 1.5 * sine .
-1.0
( scratchpad ) pi 1.5 * taylor-sine .
-1.00000000066278
( scratchpad )

Graham said

May 21, 2011 at 3:05 AM

My sin_taylor() works quite well, but my recursive sin_limitruns into numerical trouble if epsilon is too small…

#!/usr/bin/env python

from __future__ import division
from itertools import count, takewhile
from operator import mul

def factorial(n):
    return reduce(mul, xrange(2, n + 1), 1)

def sin_taylor(x, eps):
    a = lambda k: pow(x, 2 * k + 1) / factorial(2 * k + 1)
    return sum(pow(-1, k) * a(k) for k in takewhile(lambda k: a(k) > eps,
        count()))

def sin_limit(x, eps):
    #Loses accuracy if eps < 1e-6...
    if x < eps:
        return eps
    else:
        return 3 * sin_limit(x / 3, eps) - 4 * pow(sin_limit(x / 3, eps), 3)

ardnew said

October 7, 2011 at 4:40 PM

like Matías Giovannini, i went for the purely iterative implementations

use strict;
use warnings;

our $ABSTOL = 0.0000000001;
our $PI2    = 6.28318530717958648;

#
# evaluates the taylor expansion of sin(x) by
# first performing range reduction, and then
# updating the numerator, denominator, and 
# partial sum at each iteration.
#
sub taylor_sine
{  
  my $x = shift;
  
  # sine is periodic, so map domain to [-x, x]
  $x -= int($x / $PI2) * $PI2;

  # local vars  
  my ($s, $t, $n, $d, $k) = 
     ($x,  0, $x,  1,  1);
    
  do
  {
    # save previous term
    $t = $s;
    
    # update numerator and denominator
    $n *= $x * $x;
    $d *= ++$k;
    $d *= ++$k;
    
    # divide (odd) iterator by 2 and check if odd
    if (($k >> 1) & 1) { $s -= $n / $d }
                  else { $s += $n / $d }
    
  }  
  until abs($t - $s) <= $ABSTOL;
  
  return $s;
}

#
# iteratively calculates sin(x) through the triple
# angle formula, first performing range reduction
#
# MUCH faster than the recursive implementation
#
sub triple_sine
{
  my $x = shift;  
  
  # sine is periodic, so map domain to [-x, x]
  $x -= int($x / $PI2) * $PI2;
  
  # recursion depth counter
  my $n = 0;
  
  # find our base case
  while (abs($x) > $ABSTOL and ++$n)
  {
    $x /= 3;    
  }
  
  # evaluate back up our expression tree
  while ($n--)
  {
    $x = 3 * $x - 4 * $x * $x * $x;
  }
  
  return $x;
}

die "\nusage:\n\t\tperl $0 <degrees>\n"
  unless scalar @ARGV and $ARGV[0] =~ /^[-0-9.]+$/;

printf("\ntaylor_sine(%f) = %4.12f\n", $ARGV[0], 
  taylor_sine($ARGV[0]));
  
printf("\ntriple_sine(%f) = %4.12f\n", $ARGV[0], 
  triple_sine($ARGV[0]));

ardnew said
October 7, 2011 at 4:44 PM
oops, usage says to provide input in degrees, but that should be radians

Programming Praxis