home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / nihcl-3.0 / ex / Partial.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-15  |  4.5 KB  |  185 lines
  1. //  Partial.c --  vector of partial derivatives
  2.  
  3. //$Header: /afs/alw.nih.gov/unix/sun4_40c/usr/local/src/nihcl-3.0/share/ex/RCS/Partial.c,v 3.0 90/05/15 22:43:54 kgorlen Rel $
  4.  
  5. /*
  6. Author:
  7.     S.M.Orlow
  8.     Systex, Inc.
  9.     Beltsville, MD 20705
  10.     301-474-0111
  11.     sandy@alw.nih.gov
  12.  
  13. Function:
  14.  
  15. Class Partial implements a method of computing an arbitrary
  16. number of partial derivatives from a known function. The
  17. function must be realizable as a syntactically correct
  18. expression for the C-language.  Analytically, the function must
  19. be a member of the set of functions generated by a fixed number
  20. of independent variables, and exp, log, sin, and cos and closed
  21. under rational operations and function composition.
  22.  
  23. The number (order) of independent variables (partial
  24. derivatives) is set by constructors. This number is restricted
  25. to be <= the maximum order, MAX_ORD, which may be changed by
  26. recompiling Partial.c.
  27.  
  28. The member variable, 'du', will be an array of MAX_ORD+1 double
  29. values representing one function value and a list of partial
  30. derivative values.  Class Partial has been implemented in this
  31. way, rather than using dynamic allocation to allow for an
  32. indeterminate number of variables, to improve speed of execution
  33. of member functions by eliminating the need for repeated call to
  34. allocate/free memory each time a Partial object is constructed.
  35.  
  36. Class Partial is not intended to be used directly in an
  37. application program. Instead class ArrayPartial should be used
  38. to construct multiple Partial objects.
  39.  
  40. */
  41.  
  42. #include <libc.h>
  43. #include <iostream.h>
  44. #include "Partial.h"
  45.  
  46. void err_exit(char* msg) 
  47. {
  48.     cerr << " Partial: " << msg << endl;
  49.     abort();
  50. }
  51.  
  52. int max_order(const Partial& x, const Partial& y)
  53. {
  54.     int xord = x.order(), yord=y.order();
  55.     return (xord>yord)? xord:yord;
  56. }
  57. Partial::Partial()
  58. {
  59.     ord = -1;
  60.     for(int i=0; i<=MAX_ORD; i++)
  61.       du[i] = 0;
  62. }
  63. int Partial::order(int neword)
  64. {
  65.     if ( ord>=0 ) err_exit("can't reassign order");
  66.     return (ord=neword);
  67. }
  68. Partial::Partial(double val,int neword)
  69. {
  70.     if ( (neword<0)||(neword>MAX_ORD) )
  71.       err_exit("order out of bounds");
  72.     for(int i=0; i<=MAX_ORD; i++)
  73.       du[i] = 0;
  74.     ord = neword;
  75.     du[0] = val;    
  76. }
  77.     
  78. Partial::Partial(const Partial& x) 
  79. {
  80.     *this = x;
  81. }
  82.  
  83. void Partial::operator=(const Partial& x) 
  84. {
  85.     ord = x.ord;
  86.         for(int i=0; i<=MAX_ORD; i++)
  87.        if ( i<=ord ) du[i] = x[i];
  88.         else du[i] = 0;
  89. }
  90.     
  91. Partial operator+(const Partial& x,const Partial& y) 
  92. {
  93.     Partial v(0,max_order(x,y));
  94.     for(int i=0; i<=v.order(); i++) 
  95.        v[i] = x[i]+y[i];
  96.     return v;
  97. }
  98. Partial operator-(const Partial&x, const Partial& y) 
  99. {
  100.     return (x + (-y));
  101. }
  102. Partial operator*(const Partial& x,const Partial& y) 
  103. {
  104.     Partial v(x[0]*y[0],max_order(x,y));
  105.     for(int i=1; i<=v.order(); i++)
  106.        v[i] = x[0]*y[i] + x[i]*y[0];
  107.     return v;
  108. }
  109. Partial operator/(const Partial& x,const Partial& y) 
  110. {
  111.     if ( y[0] == 0 ) err_exit("division by zero");
  112.     Partial v(x[0]/y[0],max_order(x,y));
  113.     for(int i=1; i<=v.order(); i++) 
  114.        v[i] = (y[0]*x[i] - x[0]*y[i])/(y[0]*y[0]);
  115.     return v;
  116. }
  117.  
  118. Partial Partial::operator-() const 
  119. {
  120.     Partial v(0,ord);
  121.     for(int i=0; i<=ord; i++)
  122.       v[i] = - du[i];
  123.     return v;
  124. }
  125.  
  126. Partial Partial::pow(int k) const 
  127. {
  128.     if ( k==0 ) return Partial(1.0);
  129.     if ( k==1 ) return *this;
  130.  
  131.     Partial v(*this);
  132.     double du0k_1 = 1;        // to hold (k-1)th power
  133.     for(int kk=1; kk<=k-1; kk++)
  134.       du0k_1 *= du[0];
  135.     v[0] = du0k_1*du[0];        // kth power
  136.     for(int i=1; i<=ord; i++)
  137.         v[i] = k*du0k_1*du[i];
  138.     return v;
  139. }
  140.  
  141. void Partial::printOn(ostream& strm) const
  142. {
  143.     strm << "[";
  144.     for(int i=0; i<=ord; i++) {
  145.       strm << du[i];
  146.       if ( i<ord ) strm << " ";
  147.       }
  148.     strm << "]";
  149. }
  150.  
  151. ostream& operator<<(ostream& strm, const Partial& p)
  152. {
  153.     p.printOn(strm);
  154.     return strm;
  155. }
  156. Partial exp(const Partial& x) 
  157. {
  158.     Partial v(::exp(x[0]),x.order());
  159.     for(int i=1; i<=x.order(); i++)
  160.       v[i] = x[i]*v[0];    // d(exp(F)) = dF*exp(F)
  161.     return v;
  162. }
  163. Partial log(const Partial& x) 
  164. {
  165.     if ( x[0] == 0 ) err_exit("0 argument for log");
  166.     Partial v(::log(x[0]),x.order());
  167.     for(int i=1; i<=x.order(); i++)
  168.       v[i] = x[i]/x[0];
  169.     return v;
  170. }
  171. Partial sin(const Partial& x) {
  172.     Partial v(::sin(x[0]),x.order());
  173.     double cosx = ::cos(x[0]);
  174.     for(int i=1; i<=x.order(); i++)
  175.       v[i] = x[i]*cosx;
  176.     return v;
  177. }
  178. Partial cos(const Partial& x) {
  179.     Partial v(::cos(x[0]),x.order());
  180.     double sinx = ::sin(x[0]);
  181.     for(int i=1; i<=x.order(); i++)
  182.       v[i] = -x[i]*sinx;
  183.     return v;
  184. }
  185.