Appendix D. Computing Sine Using a Power Series

I l @ ve RuBoard

This program is designed to compute the sine function using a power series. A very limited floating-point format is used to demonstrate some of the problems that can occur when using floating point.

The program computes each term in the power series and displays the result. It continues computing terms until the last term is so small that it doesn't contribute to the final result. For comparison purposes, the result of the library function sin is displayed as well as the computed sine.

The program is invoked by:

 sine  value  

where value is an angle in radians. For example, to compute sin(0) we use the command:

 %  sine 0  x**1     0.000E+00 1!       1.000E+00 x**1/1! 0.000E+00 1 term computed sin(0.000E+00)=   0.000E+00 Actual sin(0)=0 

And to compute sin( figs/u03c0.gif ) we use the command:

 %  sine 3.141  x**1     3.141E+00 1!       1.000E+00 x**1/1! 3.141E+00   total   3.141E+00 x**3     3.099E+01 3!       6.000E+00 x**3/3! 5.165E+00   total   -2.024E+00 x**5     3.057E+02 5!       1.200E+02 x**5/5! 2.548E+00   total   5.239E-01 x**7     3.016E+03 7!       5.040E+03 x**7/7! 5.985E-01   total   -7.457E-02 x**9     2.976E+04 9!       3.629E+05 x**9/9! 8.201E-02   total   7.438E-03 x**11     2.936E+05 11!       3.992E+07 x**11/11! 7.355E-03   total   8.300E-05 x**13     2.897E+06 13!       6.227E+09 x**13/13! 4.652E-04   total   5.482E-04 x**15     2.858E+07 15!       1.308E+12 x**15/15! 2.185E-05   total   5.263E-04 x**17     2.819E+08 17!       3.557E+14 x**17/17! 7.927E-07   total   5.271E-04 x**19     2.782E+09 19!       1.217E+17 x**19/19! 2.287E-08   total   5.271E-04 x**21     2.744E+10 21!       5.109E+19 x**21/21! 5.371E-10 11 term computed sin(3.141E+00)=   5.271E-04 Actual sin(3.141)=0.000592654 

Example D-1 lists the Makefile for Unix.

Example D-1. sin/makefile.unx
 # # Makefile for many Unix compilers using the # "standard" command name CC # CC=CC CFLAGS=-g sine: sine.cpp         $(CC) $(CFLAGS) -o sine sine.cpp -lm clean:         rm sine 

Example D-2 lists the sine.cpp file.

Example D-2. sin/sine.cpp
 /********************************************************  * sine -- compute sine using very simple floating      *  *      arithmetic.                                     *  *                                                      *  * Usage:                                               *  *      sine <value>                                    *  *                                                      *  *      <value> is an angle in radians                  *  *                                                      *  * Format used in f.fffe+X                              *  *                                                      *  * f.fff is a 4 digit fraction                          *  *      + is a sign (+ or -)                            *  *      X is a single digit exponent                    *  *                                                      *  * sin(x) = x  - x**3 + x**5 - x**7                     *  *              -----   ----   ---- . . . .             *  *                3!     5!     7!                      *  *                                                      *  * Warning: This program is intended to show some of    *  *      problems with floating point.  It not intended  *  *      to be used to produce exact values for the      *  *      sin function.                                   *  *                                                      *  * Note: Even though we specify only one-digit for the  *  *       exponent, two are used for some calculations.  *  *       This is due to the fact that printf has no     *  *       format for a single digit exponent.            *  ********************************************************/ #include <iostream> #include <cstdlib>       #include <cmath> #include <cstdio> int main(int argc, char *argv[]) {     float   total;  // total of series so far      float   new_total;// newer version of total      float   term_top;// top part of term     float   term_bottom;// bottom of current term      float   term;   // current term      float   exp;    // exponent of current term      float   sign;   // +1 or -1 (changes on each term)      float   value;  // value of the argument to sin      int     index;  // index for counting terms      char    *float_2_ascii(float number);  // turn floating-point to ascii      float   fix_float(float number);       // round to correct digits      float   factorial(float number);       // compute n!      if (argc != 2) {         std::cerr << "Usage is:\n";         std::cerr << "  sine <value>\n";         exit (8);     }     value = fix_float(atof(&argv[1][0]));     total = 0.0;     exp = 1.0;     sign = 1.0;     for (index = 0; /* take care of below */ ; ++index) {         term_top = fix_float(pow(value, exp));         term_bottom = fix_float(factorial(exp));         term = fix_float(term_top / term_bottom);         std::cout << "x**" << static_cast<int>(exp) << "     " <<                  float_2_ascii(term_top) << '\n';         std::cout << exp << "!       " << float_2_ascii(term_bottom) << '\n';         std::cout << "x**" << static_cast<int>(exp) << "/" <<                  static_cast<int>(exp) << "! " <<                  float_2_ascii(term) << "\n";         new_total = fix_float(total + sign * term);         if (new_total == total)                 break;         total = new_total;         sign = -sign;         exp = exp + 2.0;         std::cout <<"  total   " << float_2_ascii(total) << '\n';         std::cout <<'\n';     }     std::cout << (index +1) << " term computed\n";     std::cout << "sin(" << float_2_ascii(value) << ")=\n";     std::cout << "  " << float_2_ascii(total) << '\n';     std::cout << "Actual sin(" << atof(&argv[1][0]) << ")=" <<              sin(atof(&argv[1][0])) << '\n';     return (0); } /********************************************************  * float_2_ascii -- turn a floating-point string        *  *      into ascii.                                     *  *                                                      *  * Parameters                                           *  *      number -- number to turn into ascii             *  *                                                      *  * Returns                                              *  *      Pointer to the string containing the number     *  *                                                      *  * Warning: Uses static storage, so later calls         *  *              overwrite earlier entries               *  ********************************************************/ char *float_2_ascii(float number) {     static char result[10];     //place to put the number      std::sprintf(result,"%8.3E", number);     return (result); } /********************************************************  * fix_float -- turn high precision numbers into        *  *              low precision numbers to simulate a     *  *              very dumb floating-point structure.     *  *                                                      *  * Parameters                                           *  *      number -- number to take care of                *  *                                                      *  * Returns                                              *  *      number accurate to 5 places only                *  *                                                      *  * Note: This works by changing a number into ascii and *  *       back.  Very slow, but it works.                *  ********************************************************/ float fix_float(float number) {     float   result;     // result of the conversion      char    ascii[10];  // ascii version of number      std::sprintf(ascii,"%8.4e", number);     std::sscanf(ascii, "%e", &result);     return (result); } /********************************************************  * factorial -- compute the factorial of a number.      *  *                                                      *  * Parameters                                           *  *      number -- number to use for factorial           *  *                                                      *  * Returns                                              *  *      factorial(number) or number!                    *  *                                                      *  * Note: Even though this is a floating-point routine,  *  *       using numbers that are not whole numbers       *  *       does not make sense.                           *  ********************************************************/ float factorial(float number) {     if (number <= 1.0)         return (number);     else         return (number *factorial(number - 1.0)); } 
I l @ ve RuBoard


Practical C++ Programming
Practical C Programming, 3rd Edition
ISBN: 1565923065
EAN: 2147483647
Year: 2003
Pages: 364

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net