drivers/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   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 }
/* */
root/drivers/FPU-emu/poly_atan.c

DEFINITIONS
