#include <iostream.h>
#include <math.h>

/*  A program to illustrate computer rounding.  One annoying thing about C
    relative to FORTRAN is the fact that C does not have h^n (power) as a 
    built in operation.  This program is "quick and dirty".  I don't think
    I will need to look at it later so I haven't written it carefully.      */

float summand(int i, int j );

int main() {

   double h, r_direct, r_asymp, r_rewrite, d2, d1;   
   double fact_2, fact_4, fact_6, fact_8, fact_10, fact_12, fact_14;
   fact_2  = 2;             
   fact_4  = fact_2*3*4;          //   fact_k is k!.  
   fact_6  = fact_4*5*6;          //   The computer automatically converts
   fact_8  = fact_6*7*8;          //   the integers 2, 4, ... into doubles
   fact_10 = fact_8*9*10;         //   so this gives the most accurate double
   fact_12 = fact_10*11*12;       //   precision values.
   fact_14 = fact_12*13*14;    



  while(1) {               // This is an infinite loop, type ^C to get out.
   cout << "Type h: ";
   cin >> h;
   r_direct = ( 1. - cos(h) ) / (h*h);                      
   r_rewrite = 2.*sin(h/2)*sin(h/2) / (h*h);
   r_asymp  = 1./fact_2 - h*h/fact_4 + h*h*h*h/fact_6        // This is the
             - h*h*h*h*h*h/fact_8 + h*h*h*h*h*h*h*h/fact_10  // Taylor series
             - h*h*h*h*h*h*h*h*h*h/fact_12                   // approximation.
             + h*h*h*h*h*h*h*h*h*h*h*h/fact_14;              

   d1 = r_direct - r_asymp;  
   d2 = r_rewrite - r_asymp;
   cout << "For h = " << h << " r_direct = " << r_direct << 
           ",  and r_direct - r_asymp = " << d1 << endl; 
   cout << "For h = " << h << " r_rewrite = " << r_rewrite << 
           ",  and r_rewrite - r_asymp = " << d2 << endl;}

}