root/arch/i386/math-emu/poly_atan.c

/* */

DEFINITIONS

1. poly_atan

   1 /*---------------------------------------------------------------------------+
   2  |  poly_atan.c                                                              |
   3  |                                                                           |
   4  | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
   5  |                                                                           |
   6  | Copyright (C) 1992,1993,1994                                              |
   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 "status_w.h"
  17 #include "control_w.h"
  18 #include "poly.h"
  19 
  20 
  21 #define HIPOWERon       6       /* odd poly, negative terms */
  22 static const unsigned long long oddnegterms[HIPOWERon] =
  23 {
  24   0x0000000000000000LL, /* Dummy (not for - 1.0) */
  25   0x015328437f756467LL,
  26   0x0005dda27b73dec6LL,
  27   0x0000226bf2bfb91aLL,
  28   0x000000ccc439c5f7LL,
  29   0x0000000355438407LL
  30 } ;
  31 
  32 #define HIPOWERop       6       /* odd poly, positive terms */
  33 static const unsigned long long oddplterms[HIPOWERop] =
  34 {
  35 /*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
  36   0x0db55a71875c9ac2LL,
  37   0x0029fce2d67880b0LL,
  38   0x0000dfd3908b4596LL,
  39   0x00000550fd61dab4LL,
  40   0x0000001c9422b3f9LL,
  41   0x000000003e3301e1LL
  42 };
  43 
  44 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
  45 
  46 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
  47 
  48 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
  49 
  50 
  51 /*--- poly_atan() -----------------------------------------------------------+
  52  |                                                                           |
  53  +---------------------------------------------------------------------------*/
  54 void    poly_atan(FPU_REG *arg1, FPU_REG *arg2, FPU_REG *result)
     /*  */
  55 {
  56   char                  transformed, inverted,
  57                         sign1 = arg1->sign, sign2 = arg2->sign;
  58   long int              exponent, dummy_exp;
  59   Xsig                  accumulator, Numer, Denom, accumulatore, argSignif,
  60                         argSq, argSqSq;
  61   
  62 
  63   arg1->sign = arg2->sign = SIGN_POS;
  64   if ( (compare(arg2) & ~COMP_Denormal) == COMP_A_lt_B )
  65     {
  66       inverted = 1;
  67       exponent = arg1->exp - arg2->exp;
  68       Numer.lsw = Denom.lsw = 0;
  69       XSIG_LL(Numer) = significand(arg1);
  70       XSIG_LL(Denom) = significand(arg2);
  71     }
  72   else
  73     {
  74       inverted = 0;
  75       exponent = arg2->exp - arg1->exp;
  76       Numer.lsw = Denom.lsw = 0;
  77       XSIG_LL(Numer) = significand(arg2);
  78       XSIG_LL(Denom) = significand(arg1);
  79      }
  80   div_Xsig(&Numer, &Denom, &argSignif);
  81   exponent += norm_Xsig(&argSignif);
  82 
  83   if ( (exponent >= -1)
  84       || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
  85     {
  86       /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
  87       /* Convert the argument by an identity for atan */
  88       transformed = 1;
  89 
  90       if ( exponent >= 0 )
  91         {
  92 #ifdef PARANOID
  93           if ( !( (exponent == 0) && 
  94                  (argSignif.lsw == 0) && (argSignif.midw == 0) &&
  95                  (argSignif.msw == 0x80000000) ) )
  96             {
  97               EXCEPTION(EX_INTERNAL|0x104);  /* There must be a logic error */
  98               return;
  99             }
 100 #endif PARANOID
 101           argSignif.msw = 0;   /* Make the transformed arg -> 0.0 */
 102         }
 103       else
 104         {
 105           Numer.lsw = Denom.lsw = argSignif.lsw;
 106           XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
 107 
 108           if ( exponent < -1 )
 109             shr_Xsig(&Numer, -1-exponent);
 110           negate_Xsig(&Numer);
 111       
 112           shr_Xsig(&Denom, -exponent);
 113           Denom.msw |= 0x80000000;
 114       
 115           div_Xsig(&Numer, &Denom, &argSignif);
 116 
 117           exponent = -1 + norm_Xsig(&argSignif);
 118         }
 119     }
 120   else
 121     {
 122       transformed = 0;
 123     }
 124 
 125   argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw;
 126   argSq.msw = argSignif.msw;
 127   mul_Xsig_Xsig(&argSq, &argSq);
 128   
 129   argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw;
 130   mul_Xsig_Xsig(&argSqSq, &argSqSq);
 131 
 132   accumulatore.lsw = argSq.lsw;
 133   XSIG_LL(accumulatore) = XSIG_LL(argSq);
 134 
 135   shr_Xsig(&argSq, 2*(-1-exponent-1));
 136   shr_Xsig(&argSqSq, 4*(-1-exponent-1));
 137 
 138   /* Now have argSq etc with binary point at the left
 139      .1xxxxxxxx */
 140 
 141   /* Do the basic fixed point polynomial evaluation */
 142   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 143   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
 144                    oddplterms, HIPOWERop-1);
 145   mul64_Xsig(&accumulator, &XSIG_LL(argSq));
 146   negate_Xsig(&accumulator);
 147   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
 148   negate_Xsig(&accumulator);
 149   add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
 150 
 151   mul64_Xsig(&accumulatore, &denomterm);
 152   shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
 153   accumulatore.msw |= 0x80000000;
 154 
 155   div_Xsig(&accumulator, &accumulatore, &accumulator);
 156 
 157   mul_Xsig_Xsig(&accumulator, &argSignif);
 158   mul_Xsig_Xsig(&accumulator, &argSq);
 159 
 160   shr_Xsig(&accumulator, 3);
 161   negate_Xsig(&accumulator);
 162   add_Xsig_Xsig(&accumulator, &argSignif);
 163 
 164   if ( transformed )
 165     {
 166       /* compute pi/4 - accumulator */
 167       shr_Xsig(&accumulator, -1-exponent);
 168       negate_Xsig(&accumulator);
 169       add_Xsig_Xsig(&accumulator, &pi_signif);
 170       exponent = -1;
 171     }
 172 
 173   if ( inverted )
 174     {
 175       /* compute pi/2 - accumulator */
 176       shr_Xsig(&accumulator, -exponent);
 177       negate_Xsig(&accumulator);
 178       add_Xsig_Xsig(&accumulator, &pi_signif);
 179       exponent = 0;
 180     }
 181 
 182   if ( sign1 )
 183     {
 184       /* compute pi - accumulator */
 185       shr_Xsig(&accumulator, 1 - exponent);
 186       negate_Xsig(&accumulator);
 187       add_Xsig_Xsig(&accumulator, &pi_signif);
 188       exponent = 1;
 189     }
 190 
 191   exponent += round_Xsig(&accumulator);
 192   significand(result) = XSIG_LL(accumulator);
 193   result->exp = exponent + EXP_BIAS;
 194   result->tag = TW_Valid;
 195   result->sign = sign2;
 196 
 197 }

/* */