home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / polyn / euclid.c next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.8 KB  |  254 lines
  1. /*                            euclid.c
  2.  *
  3.  *    Rational arithmetic routines
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * 
  10.  * typedef struct
  11.  *      {
  12.  *      double n;  numerator
  13.  *      double d;  denominator
  14.  *      }fract;
  15.  *
  16.  * radd( a, b, c )      c = b + a
  17.  * rsub( a, b, c )      c = b - a
  18.  * rmul( a, b, c )      c = b * a
  19.  * rdiv( a, b, c )      c = b / a
  20.  * euclid( &n, &d )     Reduce n/d to lowest terms,
  21.  *                      return greatest common divisor.
  22.  *
  23.  * Arguments of the routines are pointers to the structures.
  24.  * The double precision numbers are assumed, without checking,
  25.  * to be integer valued.  Overflow conditions are reported.
  26.  */
  27.  
  28.  
  29. #include "mconf.h"
  30.  
  31. extern double MACHEP;
  32. #define BIG (1.0/MACHEP)
  33.  
  34. typedef struct
  35.     {
  36.     double n; /* numerator */
  37.     double d; /* denominator */
  38.     }fract;
  39.  
  40.  
  41. double euclid();
  42.  
  43.  
  44.  
  45.  
  46. /* Add fractions. */
  47.  
  48. radd( f1, f2, f3 )
  49. fract *f1, *f2, *f3;
  50. {
  51. double gcd, d1, d2, gcn, n1, n2;
  52.  
  53. n1 = f1->n;
  54. d1 = f1->d;
  55. n2 = f2->n;
  56. d2 = f2->d;
  57. if( n1 == 0.0 )
  58.     {
  59.     f3->n = n2;
  60.     f3->d = d2;
  61.     return;
  62.     }
  63. if( n2 == 0.0 )
  64.     {
  65.     f3->n = n1;
  66.     f3->d = d1;
  67.     return;
  68.     }
  69.  
  70. gcd = euclid( &d1, &d2 ); /* common divisors of denominators */
  71. gcn = euclid( &n1, &n2 ); /* common divisors of numerators */
  72. /* Note, factoring the numerators
  73.  * makes overflow slightly less likely.
  74.  */
  75. f3->n = ( n1 * d2 + n2 * d1) * gcn;
  76. f3->d = d1 * d2 * gcd;
  77. euclid( &f3->n, &f3->d );
  78. }
  79.  
  80.  
  81. /* Subtract fractions. */
  82.  
  83. rsub( f1, f2, f3 )
  84. fract *f1, *f2, *f3;
  85. {
  86. double gcd, d1, d2, gcn, n1, n2;
  87.  
  88. n1 = f1->n;
  89. d1 = f1->d;
  90. n2 = f2->n;
  91. d2 = f2->d;
  92. if( n1 == 0.0 )
  93.     {
  94.     f3->n = n2;
  95.     f3->d = d2;
  96.     return;
  97.     }
  98. if( n2 == 0.0 )
  99.     {
  100.     f3->n = -n1;
  101.     f3->d = d1;
  102.     return;
  103.     }
  104.  
  105. gcd = euclid( &d1, &d2 );
  106. gcn = euclid( &n1, &n2 );
  107. f3->n = (n2 * d1 - n1 * d2) * gcn;
  108. f3->d = d1 * d2 * gcd;
  109. euclid( &f3->n, &f3->d );
  110. }
  111.  
  112.  
  113.  
  114.  
  115. /* Multiply fractions. */
  116.  
  117. rmul( ff1, ff2, ff3 )
  118. fract *ff1, *ff2, *ff3;
  119. {
  120. double d1, d2, n1, n2;
  121. double fabs();
  122.  
  123. n1 = ff1->n;
  124. d1 = ff1->d;
  125. n2 = ff2->n;
  126. d2 = ff2->d;
  127.  
  128. if( (n1 == 0.0) || (n2 == 0.0) )
  129.     {
  130.     ff3->n = 0.0;
  131.     ff3->d = 1.0;
  132.     return;
  133.     }
  134. euclid( &n1, &d2 ); /* cross cancel common divisors */
  135. euclid( &n2, &d1 );
  136. ff3->n = n1 * n2;
  137. ff3->d = d1 * d2;
  138. /* Report overflow. */
  139. if( (fabs(ff3->n) >= BIG) || (fabs(ff3->d) >= BIG) )
  140.     {
  141.     mtherr( "rmul", OVERFLOW );
  142.     return(1.0);
  143.     }
  144. /* euclid( &ff3->n, &ff3->d );*/
  145. }
  146.  
  147.  
  148.  
  149. /* Divide fractions. */
  150.  
  151. rdiv( ff1, ff2, ff3 )
  152. fract *ff1, *ff2, *ff3;
  153. {
  154. double d1, d2, n1, n2;
  155. double fabs();
  156.  
  157. n1 = ff1->d;    /* Invert ff1, then multiply */
  158. d1 = ff1->n;
  159. if( d1 < 0.0 )
  160.     { /* keep denominator positive */
  161.     n1 = -n1;
  162.     d1 = -d1;
  163.     }
  164. n2 = ff2->n;
  165. d2 = ff2->d;
  166. if( (n1 == 0.0) || (n2 == 0.0) )
  167.     {
  168.     ff3->n = 0.0;
  169.     ff3->d = 1.0;
  170.     return;
  171.     }
  172.  
  173. euclid( &n1, &d2 ); /* cross cancel any common divisors */
  174. euclid( &n2, &d1 );
  175. ff3->n = n1 * n2;
  176. ff3->d = d1 * d2;
  177. /* Report overflow. */
  178. if( (fabs(ff3->n) >= BIG) || (fabs(ff3->d) >= BIG) )
  179.     {
  180.     mtherr( "rdiv", OVERFLOW );
  181.     return(1.0);
  182.     }
  183. /* euclid( &ff3->n, &ff3->d );*/
  184. }
  185.  
  186.  
  187.  
  188.  
  189.  
  190. /* Euclidean algorithm
  191.  *   reduces fraction to lowest terms,
  192.  *   returns greatest common divisor.
  193.  */
  194.  
  195.  
  196. double euclid( num, den )
  197. double *num, *den;
  198. {
  199. double n, d, q, r;
  200. double floor();
  201.  
  202. n = *num; /* Numerator. */
  203. d = *den; /* Denominator. */
  204.  
  205. /* Make numbers positive, locally. */
  206. if( n < 0.0 )
  207.     n = -n;
  208. if( d < 0.0 )
  209.     d = -d;
  210.  
  211. /* Abort if numbers are too big for integer arithmetic. */
  212. if( (n >= BIG) || (d >= BIG) )
  213.     {
  214.     mtherr( "euclid", OVERFLOW );
  215.     return(1.0);
  216.     }
  217.  
  218. /* Divide by zero, gcd = 1. */
  219. if(d == 0.0)
  220.     return( 1.0 );
  221.  
  222. /* Zero. Return 0/1, gcd = denominator. */
  223. if(n == 0.0)
  224.     {
  225. /*
  226.     if( *den < 0.0 )
  227.         *den = -1.0;
  228.     else
  229.         *den = 1.0;
  230. */
  231.     *den = 1.0;
  232.     return( d );
  233.     }
  234.  
  235. while( d > 0.5 )
  236.     {
  237. /* Find integer part of n divided by d. */
  238.     q = floor( n/d );
  239. /* Find remainder after dividing n by d. */
  240.     r = n - d * q;
  241. /* The next fraction is d/r. */
  242.     n = d;
  243.     d = r;
  244.     }
  245.  
  246. if( n < 0.0 )
  247.     mtherr( "euclid", UNDERFLOW );
  248.  
  249. *num /= n;
  250. *den /= n;
  251. return( n );
  252. }
  253.  
  254.