home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource3 / 141_01 / cmath.c < prev    next >
Text File  |  1985-03-11  |  15KB  |  471 lines
  1. /* cmath.c 1984 apr 30 pmk - 
  2.     combine fpsqrt + clogs + ctrig 
  3.     move consts to float.h & coef's to coef.h
  4.     speedup inner eval loops */
  5. #include "float.h"
  6. #include "coef.h"
  7.  
  8. char *fpsqrt (z, xin) char *z, *xin ; {    /* z = sqrt (x)    */
  9. #define ITS 6    /* iteration count */
  10.     float x, z2, x2 ;    unsigned n ;
  11.     fpasg (z, xin) ; fpasg (x, xin) ; 
  12.     z[4] = (z[4] >> 1) | (0x80 & z[4]) ;
  13.     for (n = ITS ; n ; --n)     {
  14.         /* Newton - Raphson, ITS iterations */
  15.         fpdiv (z2, x, z) ;    /* z2 = x / z ; */
  16.         fpadd (z, z, z2) ;    /* z += z2 ;    */
  17.         fpmult (z, HALF, z) ;    /* z *= half ;  */
  18.     } /* for ITS */
  19.     return z ;
  20. #undef ITS
  21. } /* fpsqrt */
  22.  
  23.  /* CTRIG.C *****
  24.     sine, cosine, tangent and arctangent
  25.     convert degrees to radians, convert radians to degrees
  26. These functions are discussed in detail in CTRIG.DOC
  27.  
  28. L. C. Calhoun
  29. 257 South Broadway
  30. Lebanon, Ohio 45036   513 932 4541/433 7510
  31.   */
  32.  
  33. /* simple ones first converting degrees - radians */
  34.  
  35. char *degtorad(rad,deg) /*obvious arguments in 5 char fp */
  36. char *rad, *deg;
  37. {
  38.     char *fpmult();
  39.     fpmult(rad,deg,RADINDEG);
  40.     return (rad);
  41. } /* degtorad */
  42.  
  43. char *radtodeg(deg,rad) /* 5 char fp arguments */
  44. char *deg, *rad;
  45. {
  46.     char *fpmult();
  47.     fpmult(deg,rad,DEGINRAD);
  48.     return (deg);
  49. } /* radtodeg */
  50.  
  51. /* service function sinev which evaluates when range of angle
  52. reduced to plus or minus pi/2 (90 deg) */
  53. char *sinev(result,x)    char *result, x[5];    {
  54.     char *fpmult(),xsq[5];    char *coef[5];
  55.     char *fpadd(),*fpasg();
  56.     int index;
  57.  
  58.     /*  use the exponent part of the floating point
  59.         to check for threat of underflow  use small
  60.         angle approximation if appropriate  */
  61.     if ( (x[4] > 128) && (x[4] < 226) )
  62.        {fpasg(result,x);       return (result);
  63.       }   /* solution to fpmult underflow problem */
  64.  
  65. /* series coef are 1., -.1666666, .008333026, -.0001980742,
  66.     .000002601887  determined from coefset program */
  67.     coef[0] = SINC0 ;    coef[1] = SINC1 ;
  68.     coef[2] = SINC2 ;    coef[3] = SINC3 ;
  69.     coef[4] = SINC4 ;    fpmult(xsq,x,x) ;
  70.  /* to this point the coef have been initialized, 
  71.     x squared computed. */
  72.  
  73. /* now to do the polynomial approximation */
  74.     fpasg (result, coef[4]) ;
  75.     for (index = 3 ; 0 <= index ; --index) {
  76.         fpadd (result, coef[index],
  77.             fpmult (result, result, xsq)) ;
  78.     } /* this loop approx 1/3 faster than original */
  79.     fpmult (result, result, x) ;
  80.     return (result);
  81. } /* sinev */
  82.  
  83.  /* here is sine(result,angle) with angle in radians */
  84. char *sine(result,angle)    char *result, *angle;    {
  85.     char *fpmult(),*twopi,*halfpi;
  86.     char *mtwopi,*mhalfpi,*fpasg(),*fpchs();
  87.     char *pi,*sinev(),*fpadd();    char y[5],*fpsub();
  88.     int fpcomp(), compar;    int signsine;
  89.  /* some initialization required here */
  90.     signsine = 1;
  91.     twopi    = TWOPI ;    halfpi    = HALFPI ;
  92.     pi    = PI ;        mtwopi    = MTWOPI ;
  93.     mhalfpi    = MHALFPI ;    fpasg(y,angle);
  94.     while (fpcomp(y,twopi) >= 0)
  95.        fpsub(y,y,twopi);
  96.     while (fpcomp(y,mtwopi)<= 0)
  97.        fpadd(y,y,twopi);
  98.     if(fpcomp(y,halfpi) > 0)
  99.        { signsine *=-1; fpsub (y,y,pi);       }
  100.     if(fpcomp(y,mhalfpi)<  0)
  101.        { signsine *=-1; fpadd (y,y,pi);       }
  102.     sinev(result,y);
  103.     if (signsine > 0) return (result);
  104.  /* minus so need to change sign */
  105.     else return ( fpchs(result,result) );
  106. } /* sine */
  107.  
  108.  /* cosine(result,angle) with angle in radians  - uses sine */
  109. char *cosine(result,angle)    char *result, *angle;    {
  110.     char *sine(),*fpsub(),y[5];
  111.     fpsub(y,HALFPI,angle);
  112.     sine(result,y);
  113.     return (result);
  114. } /* cosine */
  115.  
  116. /* tangent(result,angle) with angle in radians, returns very 
  117. large number for divide by zero condition */
  118. char *tangent(result,angle)    char *result, angle[5];    {
  119.     char *sine(), *cosine(), *fpdiv(), zero[5];
  120.     char sresult[5], cresult[5], intres[5], big[5];
  121.     char *fpmult(), *fpmag();
  122.  
  123.     sine(sresult,angle);
  124.     cosine(cresult,angle);
  125.     /* check magnitude of denominator :*/
  126.     /* check magnitude of denominator using exponent */
  127.     if ( (cresult[4] > 128) && (cresult[4] < 132) )
  128.        {initb(big,"30,207,228,127,128"); /*big number */
  129.         if ( sresult[3] > 127 ) /*use mantissa sign */
  130.            return ( fpchs(result,big) );
  131.         else return ( fpasg(result,big) );
  132.        }
  133.     /* check for small angle, use small angle approx to
  134.        avoid underflow */
  135.     if ( (angle[4] < 226) && (angle[4] > 128) )
  136.        return ( fpasg(result,angle) );
  137.     else
  138.        return ( fpdiv(result,sresult,cresult) );
  139. } /* tangent */
  140.  
  141. /* atanev(result,x) evaluates arctangent for 0 <= x < infinity,
  142.     result in radians */
  143. char *atanev(result,x)    char  result[5], x[5];    {
  144.     char *coef[8],y[5];    int index;
  145.     char *fpadd(),*fpmult(),*atof(),ysq[5];
  146.     char yterm[5],termreslt[5],*fpasg();
  147.  
  148.     /* small angle approximation */
  149.     if ( (x[4] > 128) && (x[4] < 226) )
  150.     /* use fp exponent to check size, use small angle */
  151.      return ( fpasg(result,x) );
  152.  
  153.     fpasg(result,PIOVER4);
  154.     /* check for argument near one */
  155.     fpsub(yterm,x,ONE);
  156.     if ( (yterm[4] > 128) && (yterm[4] < 243) )
  157.        return (result);
  158.  
  159.     coef[0] = ATNC0 ;    coef[1] = ATNC1 ;
  160.     coef[2] = ATNC2 ;    coef[3] = ATNC3 ;
  161.     coef[4] = ATNC4 ;    coef[5] = ATNC5 ;
  162.     coef[6] = ATNC6 ;    coef[7] = ATNC7 ;
  163.  
  164.     fpadd(termreslt,x,ONE);
  165.     fpdiv(y,yterm,termreslt);
  166.     fpmult(ysq,y,y);
  167.  /* do poly evaluation, fast loop */
  168.     for (index = 6, fpasg (result, coef[7]) ; 0 <= index ;
  169.         --index)
  170.         fpadd (result, coef[index], 
  171.             fpmult (result, result, ysq)) ;
  172.     fpadd (result, PIOVER4, fpmult (result, result, y) ) ;
  173.     return (result);
  174. } /* atanev */
  175.  
  176. /* arctan(result,angle) is floating point arctangent evaluation */
  177. char *arctan(result,x)    char result[5], x[5];    {
  178.     char *atanev(),*fpasg(),y[5];
  179.     char *fpchs();    int index;
  180.     /* check exponent for very large argument */
  181.     if ( (x[4] > 100) && (x[4] <= 128) )
  182.         fpasg(result,HALFPI);
  183.     else  /* go through evaluation */
  184.       { fpmag(y,x);   atanev(result,y);  }
  185.  
  186.     if ( x[3] > 127 )
  187.        return ( fpchs(result,result) ); 
  188.     else return (result);
  189. } /* arctan */
  190.  
  191. char *arctan2 (result, quadrant, opside, adjside)
  192. char result[5], opside[5], adjside[5];    int *quadrant;    {
  193.     char x[5], *fpmag(), *fpchs(); *arctan();
  194.     int opsign, adjsign;
  195.     char *zero ;    zero = ZERO ;
  196.  
  197.     opsign = fpcomp(opside,zero);
  198.     adjsign = fpcomp(adjside,zero);
  199.  
  200.     if((adjsign == 0) && (opsign == 0))
  201.        { *quadrant = 0;    fpasg(result,zero);
  202.         return(result);   }
  203.  
  204.     if ( ( (128 < adjside[4]) && (adjside[4] < 226) ) 
  205.     || (adjsign == 0) )
  206.        fpasg(result,HALFPI);
  207.     else
  208.        { fpdiv(x,opside,adjside);
  209.         fpmag(x,x);    arctan(result,x);   }
  210.  
  211.     if ( (adjsign == 0) && (opsign >  0) )
  212.        { *quadrant = 1;    return(result);   }
  213.     if ( (adjsign == 0) && (opsign <  0) )
  214.        { *quadrant = 4;    fpchs(result,result);
  215.         return(result);   }
  216.     if ( (adjsign >  0) && (opsign == 0) )
  217.        { *quadrant = 1;    return(result);   }
  218.     if ( (adjsign < 0 ) && (opsign == 0) )
  219.        { *quadrant = 2;    fpchs(result,result);
  220.         return(result);   }
  221.     if ( (adjsign > 0) && (opsign > 0) )
  222.        { *quadrant = 1;    return(result);   }
  223.     if ( (adjsign > 0) && (opsign < 0) )
  224.        { *quadrant = 4;    fpchs(result,result);
  225.         return(result);   }
  226.     if ( (adjsign < 0) && (opsign > 0) )
  227.        { *quadrant = 2;    fpchs(result,result);
  228.         return(result);   }
  229.     if ( (adjsign < 0) && (opsign < 0) )
  230.        { *quadrant = 3;    return (result);   }
  231. } /* arctan2 */
  232.  
  233. /*  CLOGS.C       */
  234.  
  235. char *pi(result)    char *result;    {
  236.     char *fpasg();
  237.     return (fpasg(result,PI) );
  238. } /* pi */
  239.  
  240. char *expe(result,xin)
  241.  /* computes the base of the natural log "e" raised to the "x'th"
  242.     power.  Checks are made for out of range values and result is
  243.     defaulted to 0, 1, or a large number as appropriate.  There are
  244.     no error flags.  The arguments are floating point character
  245.     arrays char result[5], x[5]; in calling program.  Return is
  246.     a pointer to result, and "e to the x" stored in result.
  247.  */
  248. char *result, *xin ;    {
  249.     char *one, *coef[7], x[5] ;
  250.     char  *fpmult(), *fpasg(), *fpdiv(), *fpchs();
  251.     int signx, index, bigexp, smallexp, zeroin;
  252.     int exprange();
  253.  
  254.     bigexp = smallexp = zeroin = 0;    one = ONE ;
  255.  
  256. /* preserve arg before transforms */
  257.     fpasg (x, xin) ;
  258.  
  259.     coef[0] = EXPC0 ;    coef[1] = EXPC1 ;
  260.     coef[2] = EXPC2 ;    coef[3] = EXPC3 ;
  261.     coef[4] = EXP