arch/i386/math-emu/poly_atan.c

   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,1997                                         |
   7  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
   8  |                  E-mail   billm@suburbia.net                              |
   9  |                                                                           |
  10  |                                                                           |
  11  +---------------------------------------------------------------------------*/
  12
  13 #include "exception.h"
  14 #include "reg_constant.h"
  15 #include "fpu_emu.h"
  16 #include "fpu_system.h"
  17 #include "status_w.h"
  18 #include "control_w.h"
  19 #include "poly.h"
  20
  21
  22 #define HIPOWERon       6       /* odd poly, negative terms */
  23 static const unsigned long long oddnegterms[HIPOWERon] =
  24 {
  25   0x0000000000000000LL, /* Dummy (not for - 1.0) */
  26   0x015328437f756467LL,
  27   0x0005dda27b73dec6LL,
  28   0x0000226bf2bfb91aLL,
  29   0x000000ccc439c5f7LL,
  30   0x0000000355438407LL
  31 } ;
  32
  33 #define HIPOWERop       6       /* odd poly, positive terms */
  34 static const unsigned long long oddplterms[HIPOWERop] =
  35 {
  36 /*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
  37   0x0db55a71875c9ac2LL,
  38   0x0029fce2d67880b0LL,
  39   0x0000dfd3908b4596LL,
  40   0x00000550fd61dab4LL,
  41   0x0000001c9422b3f9LL,
  42   0x000000003e3301e1LL
  43 };
  44
  45 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
  46
  47 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
  48
  49 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
  50
  51
  52 /*--- poly_atan() -----------------------------------------------------------+
  53  |                                                                           |
  54  +---------------------------------------------------------------------------*/
  55 void    poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
  56                   FPU_REG *st1_ptr, u_char st1_tag)
  57 {
  58   u_char        transformed, inverted,
  59                 sign1, sign2;
  60   int           exponent;
  61   long int      dummy_exp;
  62   Xsig          accumulator, Numer, Denom, accumulatore, argSignif,
  63                 argSq, argSqSq;
  64   u_char        tag;
  65
  66   sign1 = getsign(st0_ptr);
  67   sign2 = getsign(st1_ptr);
  68   if ( st0_tag == TAG_Valid )
  69     {
  70       exponent = exponent(st0_ptr);
  71     }
  72   else
  73     {
  74       /* This gives non-compatible stack contents... */
  75       FPU_to_exp16(st0_ptr, st0_ptr);
  76       exponent = exponent16(st0_ptr);
  77     }
  78   if ( st1_tag == TAG_Valid )
  79     {
  80       exponent -= exponent(st1_ptr);
  81     }
  82   else
  83     {
  84       /* This gives non-compatible stack contents... */
  85       FPU_to_exp16(st1_ptr, st1_ptr);
  86       exponent -= exponent16(st1_ptr);
  87     }
  88
  89   if ( (exponent < 0) || ((exponent == 0) &&
  90                           ((st0_ptr->sigh < st1_ptr->sigh) ||
  91                            ((st0_ptr->sigh == st1_ptr->sigh) &&
  92                             (st0_ptr->sigl < st1_ptr->sigl))) ) )
  93     {
  94       inverted = 1;
  95       Numer.lsw = Denom.lsw = 0;
  96       XSIG_LL(Numer) = significand(st0_ptr);
  97       XSIG_LL(Denom) = significand(st1_ptr);
  98     }
  99   else
 100     {
 101       inverted = 0;
 102       exponent = -exponent;
 103       Numer.lsw = Denom.lsw = 0;
 104       XSIG_LL(Numer) = significand(st1_ptr);
 105       XSIG_LL(Denom) = significand(st0_ptr);
 106      }
 107   div_Xsig(&Numer, &Denom, &argSignif);
 108   exponent += norm_Xsig(&argSignif);
 109
 110   if ( (exponent >= -1)
 111       || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
 112     {
 113       /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
 114       /* Convert the argument by an identity for atan */
 115       transformed = 1;
 116
 117       if ( exponent >= 0 )
 118         {
 119 #ifdef PARANOID
 120           if ( !( (exponent == 0) &&
 121                  (argSignif.lsw == 0) && (argSignif.midw == 0) &&
 122                  (argSignif.msw == 0x80000000) ) )
 123             {
 124               EXCEPTION(EX_INTERNAL|0x104);  /* There must be a logic error */
 125               return;
 126             }
 127 #endif /* PARANOID */
 128           argSignif.msw = 0;   /* Make the transformed arg -> 0.0 */
 129         }
 130       else
 131         {
 132           Numer.lsw = Denom.lsw = argSignif.lsw;
 133           XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
 134
 135           if ( exponent < -1 )
 136             shr_Xsig(&Numer, -1-exponent);
 137           negate_Xsig(&Numer);
 138
 139           shr_Xsig(&Denom, -exponent);
 140           Denom.msw |= 0x80000000;
 141
 142           div_Xsig(&Numer, &Denom, &argSignif);
 143
 144           exponent = -1 + norm_Xsig(&argSignif);
 145         }
 146     }
 147   else
 148     {
 149       transformed = 0;
 150     }
 151
 152   argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw;
 153   argSq.msw = argSignif.msw;
 154   mul_Xsig_Xsig(&argSq, &argSq);
 155
 156   argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw;
 157   mul_Xsig_Xsig(&argSqSq, &argSqSq);
 158
 159   accumulatore.lsw = argSq.lsw;
 160   XSIG_LL(accumulatore) = XSIG_LL(argSq);
 161
 162   shr_Xsig(&argSq, 2*(-1-exponent-1));
 163   shr_Xsig(&argSqSq, 4*(-1-exponent-1));
 164
 165   /* Now have argSq etc with binary point at the left
 166      .1xxxxxxxx */
 167
 168   /* Do the basic fixed point polynomial evaluation */
 169   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 170   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
 171                    oddplterms, HIPOWERop-1);
 172   mul64_Xsig(&accumulator, &XSIG_LL(argSq));
 173   negate_Xsig(&accumulator);
 174   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
 175   negate_Xsig(&accumulator);
 176   add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
 177
 178   mul64_Xsig(&accumulatore, &denomterm);
 179   shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
 180   accumulatore.msw |= 0x80000000;
 181
 182   div_Xsig(&accumulator, &accumulatore, &accumulator);
 183
 184   mul_Xsig_Xsig(&accumulator, &argSignif);
 185   mul_Xsig_Xsig(&accumulator, &argSq);
 186
 187   shr_Xsig(&accumulator, 3);
 188   negate_Xsig(&accumulator);
 189   add_Xsig_Xsig(&accumulator, &argSignif);
 190
 191   if ( transformed )
 192     {
 193       /* compute pi/4 - accumulator */
 194       shr_Xsig(&accumulator, -1-exponent);
 195       negate_Xsig(&accumulator);
 196       add_Xsig_Xsig(&accumulator, &pi_signif);
 197       exponent = -1;
 198     }
 199
 200   if ( inverted )
 201     {
 202       /* compute pi/2 - accumulator */
 203       shr_Xsig(&accumulator, -exponent);
 204       negate_Xsig(&accumulator);
 205       add_Xsig_Xsig(&accumulator, &pi_signif);
 206       exponent = 0;
 207     }
 208
 209   if ( sign1 )
 210     {
 211       /* compute pi - accumulator */
 212       shr_Xsig(&accumulator, 1 - exponent);
 213       negate_Xsig(&accumulator);
 214       add_Xsig_Xsig(&accumulator, &pi_signif);
 215       exponent = 1;
 216     }
 217
 218   exponent += round_Xsig(&accumulator);
 219
 220   significand(st1_ptr) = XSIG_LL(accumulator);
 221   setexponent16(st1_ptr, exponent);
 222
 223   tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
 224   FPU_settagi(1, tag);
 225
 226   set_precision_flag_up();  /* We do not really know if up or down,
 227                                use this as the default. */
 228
 229 }