arch/i386/math-emu/poly

/* */
This source file includes following definitions.
poly_sine
poly_cos
   1 /*---------------------------------------------------------------------------+
   2  |  poly_sin.c                                                               |
   3  |                                                                           |
   4  |  Computation of an approximation of the sin function and the cosine       |
   5  |  function by a polynomial.                                                |
   6  |                                                                           |
   7  | Copyright (C) 1992,1993,1994                                              |
   8  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
   9  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  10  |                                                                           |
  11  |                                                                           |
  12  +---------------------------------------------------------------------------*/
  13 
  14 
  15 #include "exception.h"
  16 #include "reg_constant.h"
  17 #include "fpu_emu.h"
  18 #include "control_w.h"
  19 #include "poly.h"
  20 
  21 
  22 #define N_COEFF_P       4
  23 #define N_COEFF_N       4
  24 
  25 static const unsigned long long pos_terms_l[N_COEFF_P] =
  26 {
  27   0xaaaaaaaaaaaaaaabLL,
  28   0x00d00d00d00cf906LL,
  29   0x000006b99159a8bbLL,
  30   0x000000000d7392e6LL
  31 };
  32 
  33 static const unsigned long long neg_terms_l[N_COEFF_N] =
  34 {
  35   0x2222222222222167LL,
  36   0x0002e3bc74aab624LL,
  37   0x0000000b09229062LL,
  38   0x00000000000c7973LL
  39 };
  40 
  41 
  42 
  43 #define N_COEFF_PH      4
  44 #define N_COEFF_NH      4
  45 static const unsigned long long pos_terms_h[N_COEFF_PH] =
  46 {
  47   0x0000000000000000LL,
  48   0x05b05b05b05b0406LL,
  49   0x000049f93edd91a9LL,
  50   0x00000000c9c9ed62LL
  51 };
  52 
  53 static const unsigned long long neg_terms_h[N_COEFF_NH] =
  54 {
  55   0xaaaaaaaaaaaaaa98LL,
  56   0x001a01a01a019064LL,
  57   0x0000008f76c68a77LL,
  58   0x0000000000d58f5eLL
  59 };
  60 
  61 
  62 /*--- poly_sine() -----------------------------------------------------------+
  63  |                                                                           |
  64  +---------------------------------------------------------------------------*/
  65 void    poly_sine(FPU_REG const *arg, FPU_REG *result)
     /*  */
  66 {
  67   int                 exponent, echange;
  68   Xsig                accumulator, argSqrd, argTo4;
  69   unsigned long       fix_up, adj;
  70   unsigned long long  fixed_arg;
  71 
  72 
  73 #ifdef PARANOID
  74   if ( arg->tag == TW_Zero )
  75     {
  76       /* Return 0.0 */
  77       reg_move(&CONST_Z, result);
  78       return;
  79     }
  80 #endif PARANOID
  81 
  82   exponent = arg->exp - EXP_BIAS;
  83 
  84   accumulator.lsw = accumulator.midw = accumulator.msw = 0;
  85 
  86   /* Split into two ranges, for arguments below and above 1.0 */
  87   /* The boundary between upper and lower is approx 0.88309101259 */
  88   if ( (exponent < -1) || ((exponent == -1) && (arg->sigh <= 0xe21240aa)) )
  89     {
  90       /* The argument is <= 0.88309101259 */
  91 
  92       argSqrd.msw = arg->sigh; argSqrd.midw = arg->sigl; argSqrd.lsw = 0;
  93       mul64_Xsig(&argSqrd, &significand(arg));
  94       shr_Xsig(&argSqrd, 2*(-1-exponent));
  95       argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
  96       argTo4.lsw = argSqrd.lsw;
  97       mul_Xsig_Xsig(&argTo4, &argTo4);
  98 
  99       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
 100                       N_COEFF_N-1);
 101       mul_Xsig_Xsig(&accumulator, &argSqrd);
 102       negate_Xsig(&accumulator);
 103 
 104       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
 105                       N_COEFF_P-1);
 106 
 107       shr_Xsig(&accumulator, 2);    /* Divide by four */
 108       accumulator.msw |= 0x80000000;  /* Add 1.0 */
 109 
 110       mul64_Xsig(&accumulator, &significand(arg));
 111       mul64_Xsig(&accumulator, &significand(arg));
 112       mul64_Xsig(&accumulator, &significand(arg));
 113 
 114       /* Divide by four, FPU_REG compatible, etc */
 115       exponent = 3*exponent + EXP_BIAS;
 116 
 117       /* The minimum exponent difference is 3 */
 118       shr_Xsig(&accumulator, arg->exp - exponent);
 119 
 120       negate_Xsig(&accumulator);
 121       XSIG_LL(accumulator) += significand(arg);
 122 
 123       echange = round_Xsig(&accumulator);
 124 
 125       result->exp = arg->exp + echange;
 126     }
 127   else
 128     {
 129       /* The argument is > 0.88309101259 */
 130       /* We use sin(arg) = cos(pi/2-arg) */
 131 
 132       fixed_arg = significand(arg);
 133 
 134       if ( exponent == 0 )
 135         {
 136           /* The argument is >= 1.0 */
 137 
 138           /* Put the binary point at the left. */
 139           fixed_arg <<= 1;
 140         }
 141       /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
 142       fixed_arg = 0x921fb54442d18469LL - fixed_arg;
 143 
 144       XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
 145       mul64_Xsig(&argSqrd, &fixed_arg);
 146 
 147       XSIG_LL(argTo4) = XSIG_LL(argSqrd); argTo4.lsw = argSqrd.lsw;
 148       mul_Xsig_Xsig(&argTo4, &argTo4);
 149 
 150       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
 151                       N_COEFF_NH-1);
 152       mul_Xsig_Xsig(&accumulator, &argSqrd);
 153       negate_Xsig(&accumulator);
 154 
 155       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
 156                       N_COEFF_PH-1);
 157       negate_Xsig(&accumulator);
 158 
 159       mul64_Xsig(&accumulator, &fixed_arg);
 160       mul64_Xsig(&accumulator, &fixed_arg);
 161 
 162       shr_Xsig(&accumulator, 3);
 163       negate_Xsig(&accumulator);
 164 
 165       add_Xsig_Xsig(&accumulator, &argSqrd);
 166 
 167       shr_Xsig(&accumulator, 1);
 168 
 169       accumulator.lsw |= 1;  /* A zero accumulator here would cause problems */
 170       negate_Xsig(&accumulator);
 171 
 172       /* The basic computation is complete. Now fix the answer to
 173          compensate for the error due to the approximation used for
 174          pi/2
 175          */
 176 
 177       /* This has an exponent of -65 */
 178       fix_up = 0x898cc517;
 179       /* The fix-up needs to be improved for larger args */
 180       if ( argSqrd.msw & 0xffc00000 )
 181         {
 182           /* Get about 32 bit precision in these: */
 183           mul_32_32(0x898cc517, argSqrd.msw, &adj);
 184           fix_up -= adj/6;
 185         }
 186       mul_32_32(fix_up, LL_MSW(fixed_arg), &fix_up);
 187 
 188       adj = accumulator.lsw;    /* temp save */
 189       accumulator.lsw -= fix_up;
 190       if ( accumulator.lsw > adj )
 191         XSIG_LL(accumulator) --;
 192 
 193       echange = round_Xsig(&accumulator);
 194 
 195       result->exp = EXP_BIAS - 1 + echange;
 196     }
 197 
 198   significand(result) = XSIG_LL(accumulator);
 199   result->tag = TW_Valid;
 200   result->sign = arg->sign;
 201 
 202 #ifdef PARANOID
 203   if ( (result->exp >= EXP_BIAS)
 204       && (significand(result) > 0x8000000000000000LL) )
 205     {
 206       EXCEPTION(EX_INTERNAL|0x150);
 207     }
 208 #endif PARANOID
 209 
 210 }
 211 
 212 
 213 
 214 /*--- poly_cos() ------------------------------------------------------------+
 215  |                                                                           |
 216  +---------------------------------------------------------------------------*/
 217 void    poly_cos(FPU_REG const *arg, FPU_REG *result)
     /*  */
 218 {
 219   long int            exponent, exp2, echange;
 220   Xsig                accumulator, argSqrd, fix_up, argTo4;
 221   unsigned long       adj;
 222   unsigned long long  fixed_arg;
 223 
 224 
 225 #ifdef PARANOID
 226   if ( arg->tag == TW_Zero )
 227     {
 228       /* Return 1.0 */
 229       reg_move(&CONST_1, result);
 230       return;
 231     }
 232 
 233   if ( (arg->exp > EXP_BIAS)
 234       || ((arg->exp == EXP_BIAS)
 235           && (significand(arg) > 0xc90fdaa22168c234LL)) )
 236     {
 237       EXCEPTION(EX_Invalid);
 238       reg_move(&CONST_QNaN, result);
 239       return;
 240     }
 241 #endif PARANOID
 242 
 243   exponent = arg->exp - EXP_BIAS;
 244 
 245   accumulator.lsw = accumulator.midw = accumulator.msw = 0;
 246 
 247   if ( (exponent < -1) || ((exponent == -1) && (arg->sigh <= 0xb00d6f54)) )
 248     {
 249       /* arg is < 0.687705 */
 250 
 251       argSqrd.msw = arg->sigh; argSqrd.midw = arg->sigl; argSqrd.lsw = 0;
 252       mul64_Xsig(&argSqrd, &significand(arg));
 253 
 254       if ( exponent < -1 )
 255         {
 256           /* shift the argument right by the required places */
 257           shr_Xsig(&argSqrd, 2*(-1-exponent));
 258         }
 259 
 260       argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
 261       argTo4.lsw = argSqrd.lsw;
 262       mul_Xsig_Xsig(&argTo4, &argTo4);
 263 
 264       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
 265                       N_COEFF_NH-1);
 266       mul_Xsig_Xsig(&accumulator, &argSqrd);
 267       negate_Xsig(&accumulator);
 268 
 269       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
 270                       N_COEFF_PH-1);
 271       negate_Xsig(&accumulator);
 272 
 273       mul64_Xsig(&accumulator, &significand(arg));
 274       mul64_Xsig(&accumulator, &significand(arg));
 275       shr_Xsig(&accumulator, -2*(1+exponent));
 276 
 277       shr_Xsig(&accumulator, 3);
 278       negate_Xsig(&accumulator);
 279 
 280       add_Xsig_Xsig(&accumulator, &argSqrd);
 281 
 282       shr_Xsig(&accumulator, 1);
 283 
 284       /* It doesn't matter if accumulator is all zero here, the
 285          following code will work ok */
 286       negate_Xsig(&accumulator);
 287 
 288       if ( accumulator.lsw & 0x80000000 )
 289         XSIG_LL(accumulator) ++;
 290       if ( accumulator.msw == 0 )
 291         {
 292           /* The result is 1.0 */
 293           reg_move(&CONST_1, result);
 294         }
 295       else
 296         {
 297           significand(result) = XSIG_LL(accumulator);
 298       
 299           /* will be a valid positive nr with expon = -1 */
 300           *(short *)&(result->sign) = 0;
 301           result->exp = EXP_BIAS - 1;
 302         }
 303     }
 304   else
 305     {
 306       fixed_arg = significand(arg);
 307 
 308       if ( exponent == 0 )
 309         {
 310           /* The argument is >= 1.0 */
 311 
 312           /* Put the binary point at the left. */
 313           fixed_arg <<= 1;
 314         }
 315       /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
 316       fixed_arg = 0x921fb54442d18469LL - fixed_arg;
 317 
 318       exponent = -1;
 319       exp2 = -1;
 320 
 321       /* A shift is needed here only for a narrow range of arguments,
 322          i.e. for fixed_arg approx 2^-32, but we pick up more... */
 323       if ( !(LL_MSW(fixed_arg) & 0xffff0000) )
 324         {
 325           fixed_arg <<= 16;
 326           exponent -= 16;
 327           exp2 -= 16;
 328         }
 329 
 330       XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
 331       mul64_Xsig(&argSqrd, &fixed_arg);
 332 
 333       if ( exponent < -1 )
 334         {
 335           /* shift the argument right by the required places */
 336           shr_Xsig(&argSqrd, 2*(-1-exponent));
 337         }
 338 
 339       argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
 340       argTo4.lsw = argSqrd.lsw;
 341       mul_Xsig_Xsig(&argTo4, &argTo4);
 342 
 343       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
 344                       N_COEFF_N-1);
 345       mul_Xsig_Xsig(&accumulator, &argSqrd);
 346       negate_Xsig(&accumulator);
 347 
 348       polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
 349                       N_COEFF_P-1);
 350 
 351       shr_Xsig(&accumulator, 2);    /* Divide by four */
 352       accumulator.msw |= 0x80000000;  /* Add 1.0 */
 353 
 354       mul64_Xsig(&accumulator, &fixed_arg);
 355       mul64_Xsig(&accumulator, &fixed_arg);
 356       mul64_Xsig(&accumulator, &fixed_arg);
 357 
 358       /* Divide by four, FPU_REG compatible, etc */
 359       exponent = 3*exponent;
 360 
 361       /* The minimum exponent difference is 3 */
 362       shr_Xsig(&accumulator, exp2 - exponent);
 363 
 364       negate_Xsig(&accumulator);
 365       XSIG_LL(accumulator) += fixed_arg;
 366 
 367       /* The basic computation is complete. Now fix the answer to
 368          compensate for the error due to the approximation used for
 369          pi/2
 370          */
 371 
 372       /* This has an exponent of -65 */
 373       XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
 374       fix_up.lsw = 0;
 375 
 376       /* The fix-up needs to be improved for larger args */
 377       if ( argSqrd.msw & 0xffc00000 )
 378         {
 379           /* Get about 32 bit precision in these: */
 380           mul_32_32(0x898cc517, argSqrd.msw, &adj);
 381           fix_up.msw -= adj/2;
 382           mul_32_32(0x898cc517, argTo4.msw, &adj);
 383           fix_up.msw += adj/24;
 384         }
 385 
 386       exp2 += norm_Xsig(&accumulator);
 387       shr_Xsig(&accumulator, 1); /* Prevent overflow */
 388       exp2++;
 389       shr_Xsig(&fix_up, 65 + exp2);
 390 
 391       add_Xsig_Xsig(&accumulator, &fix_up);
 392 
 393       echange = round_Xsig(&accumulator);
 394 
 395       result->exp = exp2 + EXP_BIAS + echange;
 396       *(short *)&(result->sign) = 0;      /* Is a valid positive nr */
 397       significand(result) = XSIG_LL(accumulator);
 398     }
 399 
 400 #ifdef PARANOID
 401   if ( (result->exp >= EXP_BIAS)
 402       && (significand(result) > 0x8000000000000000LL) )
 403     {
 404       EXCEPTION(EX_INTERNAL|0x151);
 405     }
 406 #endif PARANOID
 407 
 408 }
/* */
root/arch/i386/math-emu/poly_sin.c

DEFINITIONS
