root/kernel/FPU-emu/poly_atan.c

/* [previous][next][first][last][top][bottom][index][help] */

DEFINITIONS

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

/* [previous][next][first][last][top][bottom][index][help] */