/*

sin(x)  =  + x/1 - x*x*x/1*2*3 + x*x*x*x*x/1*2*3*4*5 - x*x*x*x*x*x*x/1*2*3*4*5*6*7 + etc etc

We want a function that can take the approximation to as many terms as we want, so as well
as the value of X, sine will also need a parameter to induicate the maximum term to go to.


writing it sideways:

   each term has three components:
       pom (plus or minus) +1, -1, +1, -1, ...
       xpn (x to the power of n) x, x*x*x, x*x*x*x*x, ...
       fac (N factorial) 1, 1*2*3, 1*2*3*4*5, ...

      term-value        N   pom         xpn           fac
    -------------      ---  ---         ---           ---
          x             1    +1          x            1
       -xxx/123         3    -1        x*x*x          1*2*3
     +xxxxx/12345       5    +1      x*x*x*x*x        1*2*3*4*5
   -xxxxxxx/1234567     7    -1    x*x*x*x*x*x*x      1*2*3*4*5*6*7
 +xxxxxxxxx/123456789   9    +1  x*x*x*x*x*x*x*x*x    1*2*3*4*5*6*7*8*9

To progress from term N to the next term:
   N   is increased by 2
   pom is negated
   xpn is multiplied by x squared
   fac is multiplied by (N+1)*(N+2)
all of these must be parameters too.

sine(x, N, MaxN, pom, xpn, fac)
  if given the correct parameter values,
  will add up all the terms between N and MaxN, and return that is its result

When N is 1, it will compute the entire approsimation to sin(x), and will
need these initial values
  N    1
  pom +1
  xpn  x
  fac  1

*/

#include <library.h>



double sine(const double x,
            const    int N,
            const    int MaxN,
            const    int pom,
            const double xpn,
            const double fac)
{
  if (N > MaxN)
    return 0.0;

  const double thisterm   = pom * xpn / fac;
  const double otherterms = sine(x, 
                                 N+2, 
                                 MaxN,
                                 - pom, 
                                 xpn * x * x, 
                                 fac * (N+1) * (N+2));

  return thisterm + otherterms;
}



double sine(const double x, const int MaxN)
{ return sine(x, 1, MaxN, 1, x, 1); }


const double pi = acos(-1.0);

void main()
{ while (true)
  { print("degrees? ");
    const double angle = read_double() * pi / 180.0;

    print("sine up to term 3  = ");
    precise_print(sine(angle, 3), 17);
    new_line();

    print("sine up to term 7  = ");
    precise_print(sine(angle, 7), 17);
    new_line();

    print("sine up to term 11 = ");
    precise_print(sine(angle, 11), 17);
    new_line();

    print("sine up to term 15 = ");
    precise_print(sine(angle, 15), 17);
    new_line();

    print("sine up to term 33 = ");
    precise_print(sine(angle, 33), 17);
    new_line();

    print("sine up to term 99 = ");
    precise_print(sine(angle, 99), 17);
    new_line();

    print("official sine      = ");
    precise_print(sin(angle), 17);
    new_line();

    new_line(); } }