kernel/FPU-emu/poly

/* */
This source file includes following definitions.
poly_atan
poly_add_1
   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 apm233m@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 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 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 
  42 static unsigned denomterm[2] =
  43 { 0xfc4bd208, 0xea2e6612 };
  44 
  45 
  46 
  47 /*--- poly_atan() -----------------------------------------------------------+
  48  |                                                                           |
  49  +---------------------------------------------------------------------------*/
  50 void    poly_atan(FPU_REG *arg)
     /*  */
  51 {
  52   char          recursions = 0;
  53   short         exponent;
  54   FPU_REG       odd_poly, even_poly, pos_poly, neg_poly;
  55   FPU_REG       argSq;
  56   long long     arg_signif, argSqSq;
  57   
  58 
  59 #ifdef PARANOID
  60   if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
  61     { arith_invalid(arg); return; }
  62 #endif PARANOID
  63 
  64   exponent = arg->exp - EXP_BIAS;
  65   
  66   if ( arg->tag == TW_Zero )
  67     {
  68       /* Return 0.0 */
  69       reg_move(&CONST_Z, arg);
  70       return;
  71     }
  72   
  73   if ( exponent >= -2 )
  74     {
  75       /* argument is in the range  [0.25 .. 1.0] */
  76       if ( exponent >= 0 )
  77         {
  78 #ifdef PARANOID
  79           if ( (exponent == 0) && 
  80               (arg->sigl == 0) && (arg->sigh == 0x80000000) )
  81 #endif PARANOID
  82             {
  83               reg_move(&CONST_PI4, arg);
  84               return;
  85             }
  86 #ifdef PARANOID
  87           EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
  88 #endif PARANOID
  89         }
  90 
  91       /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
  92       /* convert the argument by an identity for atan */
  93       if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
  94         {
  95           FPU_REG numerator, denom;
  96 
  97           recursions++;
  98 
  99           arg_signif = *(long long *)&(arg->sigl);
 100           if ( exponent < -1 )
 101             {
 102               if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
 103                 arg_signif++;   /* round up */
 104             }
 105           *(long long *)&(numerator.sigl) = -arg_signif;
 106           numerator.exp = EXP_BIAS - 1;
 107           normalize(&numerator);                       /* 1 - arg */
 108 
 109           arg_signif = *(long long *)&(arg->sigl);
 110           if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
 111             arg_signif++;       /* round up */
 112           *(long long *)&(denom.sigl) = arg_signif;
 113           denom.sigh |= 0x80000000;                    /* 1 + arg */
 114 
 115           arg->exp = numerator.exp;
 116           reg_u_div(&numerator, &denom, arg, FULL_PRECISION);
 117 
 118           exponent = arg->exp - EXP_BIAS;
 119         }
 120     }
 121 
 122   *(long long *)&arg_signif = *(long long *)&(arg->sigl);
 123 
 124 #ifdef PARANOID
 125   /* This must always be true */
 126   if ( exponent >= -1 )
 127     {
 128       EXCEPTION(EX_INTERNAL|0x120);     /* There must be a logic error */
 129     }
 130 #endif PARANOID
 131 
 132   /* shift the argument right by the required places */
 133   if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
 134     arg_signif++;       /* round up */
 135   
 136   /* Now have arg_signif with binary point at the left
 137      .1xxxxxxxx */
 138   mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
 139   mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
 140 
 141   /* will be a valid positive nr with expon = 0 */
 142   *(short *)&(pos_poly.sign) = 0;
 143   pos_poly.exp = EXP_BIAS;
 144 
 145   /* Do the basic fixed point polynomial evaluation */
 146   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
 147              (unsigned short (*)[4])oddplterms, HIPOWERop-1);
 148   mul64((long long *)(&argSq.sigl), (long long *)(&pos_poly.sigl),
 149         (long long *)(&pos_poly.sigl));
 150 
 151   /* will be a valid positive nr with expon = 0 */
 152   *(short *)&(neg_poly.sign) = 0;
 153   neg_poly.exp = EXP_BIAS;
 154 
 155   /* Do the basic fixed point polynomial evaluation */
 156   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
 157              (unsigned short (*)[4])oddnegterms, HIPOWERon-1);
 158 
 159   /* Subtract the mantissas */
 160   *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
 161 
 162   reg_move(&pos_poly, &odd_poly);
 163   poly_add_1(&odd_poly);
 164 
 165   /* The complete odd polynomial */
 166   reg_u_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
 167   odd_poly.exp -= EXP_BIAS - 1;
 168 
 169   /* will be a valid positive nr with expon = 0 */
 170   *(short *)&(even_poly.sign) = 0;
 171 
 172   mul64((long long *)(&argSq.sigl),
 173         (long long *)(&denomterm), (long long *)(&even_poly.sigl));
 174 
 175   poly_add_1(&even_poly);
 176 
 177   reg_div(&odd_poly, &even_poly, arg, FULL_PRECISION);
 178 
 179   if ( recursions )
 180     reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION);
 181 }
 182 
 183 
 184 /* The argument to this function must be polynomial() compatible,
 185    i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
 186    be normalized.
 187    This function adds 1.0 to the (assumed positive) argument. */
 188 void poly_add_1(FPU_REG *src)
     /*  */
 189 {
 190 /* Rounding in a consistent direction produces better results
 191    for the use of this function in poly_atan. Simple truncation
 192    is used here instead of round-to-nearest. */
 193 
 194 #ifdef OBSOLETE
 195 char round = (src->sigl & 3) == 3;
 196 #endif OBSOLETE
 197 
 198 shrx(&src->sigl, 1);
 199 
 200 #ifdef OBSOLETE
 201 if ( round ) (*(long long *)&src->sigl)++;   /* Round to even */
 202 #endif OBSOLETE
 203 
 204 src->sigh |= 0x80000000;
 205 
 206 src->exp = EXP_BIAS;
 207 
 208 }
/* */
root/kernel/FPU-emu/poly_atan.c

DEFINITIONS
