home *** CD-ROM | disk | FTP | other *** search
/ PC Extra Super CD 1998 January / PCPLUS131.iso / DJGPP / V2MISC / WMEMU2S.ZIP / LINUX / DRIVERS / FPU-EMU / POLY_ATA.C < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-01  |  5.9 KB  |  205 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  p_atan.c                                                                 |
  3.  |                                                                           |
  4.  | Compute the tan of a FPU_REG, using a polynomial approximation.           |
  5.  |                                                                           |
  6.  | Copyright (C) 1992,1993                                                   |
  7.  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
  8.  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12.  
  13. #include "exception.h"
  14. #include "reg_constant.h"
  15. #include "fpu_emu.h"
  16. #include "control_w.h"
  17.  
  18.  
  19. #define    HIPOWERon    6    /* odd poly, negative terms */
  20. static unsigned const oddnegterms[HIPOWERon][2] =
  21. {
  22.   { 0x00000000, 0x00000000 }, /* for + 1.0 */
  23.   { 0x763b6f3d, 0x1adc4428 },
  24.   { 0x20f0630b, 0x0502909d },
  25.   { 0x4e825578, 0x0198ce38 },
  26.   { 0x22b7cb87, 0x008da6e3 },
  27.   { 0x9b30ca03, 0x00239c79 }
  28. } ;
  29.  
  30. #define    HIPOWERop    6    /* odd poly, positive terms */
  31. static unsigned const    oddplterms[HIPOWERop][2] =
  32. {
  33.   { 0xa6f67cb8, 0x94d910bd },
  34.   { 0xa02ffab4, 0x0a43cb45 },
  35.   { 0x04265e6b, 0x02bf5655 },
  36.   { 0x0a728914, 0x00f280f7 },
  37.   { 0x6d640e01, 0x004d6556 },
  38.   { 0xf1dd2dbf, 0x000a530a }
  39. };
  40.  
  41. static unsigned long long const denomterm = 0xea2e6612fc4bd208LL;
  42.  
  43.  
  44. /*--- poly_atan() -----------------------------------------------------------+
  45.  |                                                                           |
  46.  +---------------------------------------------------------------------------*/
  47. void    poly_atan(FPU_REG *arg)
  48. {
  49.   char        recursions = 0;
  50.   short        exponent;
  51.   FPU_REG       odd_poly, even_poly, pos_poly, neg_poly, ratio;
  52.   FPU_REG       argSq;
  53.   unsigned long long     arg_signif, argSqSq;
  54.   
  55.  
  56. #ifdef PARANOID
  57.   if ( arg->sign != 0 )    /* Can't hack a number < 0.0 */
  58.     { arith_invalid(arg); return; }  /* Need a positive number */
  59. #endif PARANOID
  60.  
  61.   exponent = arg->exp - EXP_BIAS;
  62.   
  63.   if ( arg->tag == TW_Zero )
  64.     {
  65.       /* Return 0.0 */
  66.       reg_move(&CONST_Z, arg);
  67.       return;
  68.     }
  69.   
  70.   if ( exponent >= -2 )
  71.     {
  72.       /* argument is in the range  [0.25 .. 1.0] */
  73.       if ( exponent >= 0 )
  74.     {
  75. #ifdef PARANOID
  76.       if ( (exponent == 0) && 
  77.           (arg->sigl == 0) && (arg->sigh == 0x80000000) )
  78. #endif PARANOID
  79.         {
  80.           reg_move(&CONST_PI4, arg);
  81.           return;
  82.         }
  83. #ifdef PARANOID
  84.       EXCEPTION(EX_INTERNAL|0x104);    /* There must be a logic error */
  85.       return;
  86. #endif PARANOID
  87.     }
  88.  
  89.       /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
  90.       /* convert the argument by an identity for atan */
  91.       if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
  92.     {
  93.       FPU_REG numerator, denom;
  94.  
  95.       recursions++;
  96.  
  97.       arg_signif = significand(arg);
  98.       if ( exponent < -1 )
  99.         {
  100.           if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
  101.         arg_signif++;    /* round up */
  102.         }
  103.       significand(&numerator) = -arg_signif;
  104.       numerator.exp = EXP_BIAS - 1;
  105.       normalize(&numerator);                       /* 1 - arg */
  106.  
  107.       arg_signif = significand(arg);
  108.       if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
  109.         arg_signif++;    /* round up */
  110.       significand(&denom) = arg_signif;
  111.       denom.sigh |= 0x80000000;                    /* 1 + arg */
  112.  
  113.       arg->exp = numerator.exp;
  114.       reg_u_div(&numerator, &denom, arg, FULL_PRECISION);
  115.  
  116.       exponent = arg->exp - EXP_BIAS;
  117.     }
  118.     }
  119.  
  120.   arg_signif = significand(arg);
  121.  
  122. #ifdef PARANOID
  123.   /* This must always be true */
  124.   if ( exponent >= -1 )
  125.     {
  126.       EXCEPTION(EX_INTERNAL|0x120);    /* There must be a logic error */
  127.     }
  128. #endif PARANOID
  129.  
  130.   /* shift the argument right by the required places */
  131.   if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
  132.     arg_signif++;    /* round up */
  133.   
  134.   /* Now have arg_signif with binary point at the left
  135.      .1xxxxxxxx */
  136.   mul64(&arg_signif, &arg_signif, &significand(&argSq));
  137.   mul64(&significand(&argSq), &significand(&argSq), &argSqSq);
  138.  
  139.   /* will be a valid positive nr with expon = 0 */
  140.   *(short *)&(pos_poly.sign) = 0;
  141.   pos_poly.exp = EXP_BIAS;
  142.  
  143.   /* Do the basic fixed point polynomial evaluation */
  144.   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
  145.          (unsigned short (*)[4])oddplterms, HIPOWERop-1);
  146.   mul64(&significand(&argSq), &significand(&pos_poly),
  147.     &significand(&pos_poly));
  148.  
  149.   /* will be a valid positive nr with expon = 0 */
  150.   *(short *)&(neg_poly.sign) = 0;
  151.   neg_poly.exp = EXP_BIAS;
  152.  
  153.   /* Do the basic fixed point polynomial evaluation */
  154.   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
  155.          (unsigned short (*)[4])oddnegterms, HIPOWERon-1);
  156.  
  157.   /* Subtract the mantissas */
  158.   significand(&pos_poly) -= significand(&neg_poly);
  159.  
  160.   reg_move(&pos_poly, &odd_poly);
  161.   poly_add_1(&odd_poly);
  162.  
  163.   /* will be a valid positive nr with expon = 0 */
  164.   *(short *)&(even_poly.sign) = 0;
  165.  
  166.   mul64(&significand(&argSq), &denomterm, &significand(&even_poly));
  167.  
  168.   poly_add_1(&even_poly);
  169.  
  170.   reg_div(&odd_poly, &even_poly, &ratio, FULL_PRECISION);
  171.  
  172.   reg_u_mul(&ratio, arg, arg, FULL_PRECISION);
  173.  
  174.   if ( recursions )
  175.     reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION);
  176.  
  177. }
  178.  
  179.  
  180. /* The argument to this function must be polynomial() compatible,
  181.    i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
  182.    be normalized.
  183.    This function adds 1.0 to the (assumed positive) argument. */
  184. void poly_add_1(FPU_REG *src)
  185. {
  186. /* Rounding in a consistent direction produces better results
  187.    for the use of this function in poly_atan. Simple truncation
  188.    is used here instead of round-to-nearest. */
  189.  
  190. #ifdef OBSOLETE
  191. char round = (src->sigl & 3) == 3;
  192. #endif OBSOLETE
  193.  
  194. shrx(&src->sigl, 1);
  195.  
  196. #ifdef OBSOLETE
  197. if ( round ) significand(src)++;   /* Round to even */
  198. #endif OBSOLETE
  199.  
  200. src->sigh |= 0x80000000;
  201.  
  202. src->exp = EXP_BIAS;
  203.  
  204. }
  205.