FreeBSD/Linux Kernel Cross Reference
sys/arch/i386/math-emu/poly_sin.c

     /*---------------------------------------------------------------------------+
      |  poly_sin.c                                                               |
      |                                                                           |
      |  Computation of an approximation of the sin function and the cosine       |
      |  function by a polynomial.                                                |
      |                                                                           |
      | Copyright (C) 1992,1993,1994,1997,1999                                    |
      |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
      |                  E-mail   billm@melbpc.org.au                             |
     |                                                                           |
     |                                                                           |
     +---------------------------------------------------------------------------*/
    
    
    #include "exception.h"
    #include "reg_constant.h"
    #include "fpu_emu.h"
    #include "fpu_system.h"
    #include "control_w.h"
    #include "poly.h"
    
    
    #define N_COEFF_P       4
    #define N_COEFF_N       4
    
    static const unsigned long long pos_terms_l[N_COEFF_P] =
    {
      0xaaaaaaaaaaaaaaabLL,
      0x00d00d00d00cf906LL,
      0x000006b99159a8bbLL,
      0x000000000d7392e6LL
    };
    
    static const unsigned long long neg_terms_l[N_COEFF_N] =
    {
      0x2222222222222167LL,
      0x0002e3bc74aab624LL,
      0x0000000b09229062LL,
      0x00000000000c7973LL
    };
    
    
    
    #define N_COEFF_PH      4
    #define N_COEFF_NH      4
    static const unsigned long long pos_terms_h[N_COEFF_PH] =
    {
      0x0000000000000000LL,
      0x05b05b05b05b0406LL,
      0x000049f93edd91a9LL,
      0x00000000c9c9ed62LL
    };
    
    static const unsigned long long neg_terms_h[N_COEFF_NH] =
    {
      0xaaaaaaaaaaaaaa98LL,
      0x001a01a01a019064LL,
      0x0000008f76c68a77LL,
      0x0000000000d58f5eLL
    };
    
    
    /*--- poly_sine() -----------------------------------------------------------+
     |                                                                           |
     +---------------------------------------------------------------------------*/
    void    poly_sine(FPU_REG *st0_ptr)
    {
      int                 exponent, echange;
      Xsig                accumulator, argSqrd, argTo4;
      unsigned long       fix_up, adj;
      unsigned long long  fixed_arg;
      FPU_REG             result;
    
      exponent = exponent(st0_ptr);
    
      accumulator.lsw = accumulator.midw = accumulator.msw = 0;
    
      /* Split into two ranges, for arguments below and above 1.0 */
      /* The boundary between upper and lower is approx 0.88309101259 */
      if ( (exponent < -1) || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa)) )
        {
          /* The argument is <= 0.88309101259 */
    
          argSqrd.msw = st0_ptr->sigh; argSqrd.midw = st0_ptr->sigl; argSqrd.lsw = 0;
          mul64_Xsig(&argSqrd, &significand(st0_ptr));
          shr_Xsig(&argSqrd, 2*(-1-exponent));
          argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
          argTo4.lsw = argSqrd.lsw;
          mul_Xsig_Xsig(&argTo4, &argTo4);
    
          polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
                          N_COEFF_N-1);
          mul_Xsig_Xsig(&accumulator, &argSqrd);
          negate_Xsig(&accumulator);
    
          polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
                          N_COEFF_P-1);
    
          shr_Xsig(&accumulator, 2);    /* Divide by four */
         accumulator.msw |= 0x80000000;  /* Add 1.0 */
   
         mul64_Xsig(&accumulator, &significand(st0_ptr));
         mul64_Xsig(&accumulator, &significand(st0_ptr));
         mul64_Xsig(&accumulator, &significand(st0_ptr));
   
         /* Divide by four, FPU_REG compatible, etc */
         exponent = 3*exponent;
   
         /* The minimum exponent difference is 3 */
         shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
   
         negate_Xsig(&accumulator);
         XSIG_LL(accumulator) += significand(st0_ptr);
   
         echange = round_Xsig(&accumulator);
   
         setexponentpos(&result, exponent(st0_ptr) + echange);
       }
     else
       {
         /* The argument is > 0.88309101259 */
         /* We use sin(st(0)) = cos(pi/2-st(0)) */
   
         fixed_arg = significand(st0_ptr);
   
         if ( exponent == 0 )
           {
             /* The argument is >= 1.0 */
   
             /* Put the binary point at the left. */
             fixed_arg <<= 1;
           }
         /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
         fixed_arg = 0x921fb54442d18469LL - fixed_arg;
         /* There is a special case which arises due to rounding, to fix here. */
         if ( fixed_arg == 0xffffffffffffffffLL )
           fixed_arg = 0;
   
         XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
         mul64_Xsig(&argSqrd, &fixed_arg);
   
         XSIG_LL(argTo4) = XSIG_LL(argSqrd); argTo4.lsw = argSqrd.lsw;
         mul_Xsig_Xsig(&argTo4, &argTo4);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
                         N_COEFF_NH-1);
         mul_Xsig_Xsig(&accumulator, &argSqrd);
         negate_Xsig(&accumulator);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
                         N_COEFF_PH-1);
         negate_Xsig(&accumulator);
   
         mul64_Xsig(&accumulator, &fixed_arg);
         mul64_Xsig(&accumulator, &fixed_arg);
   
         shr_Xsig(&accumulator, 3);
         negate_Xsig(&accumulator);
   
         add_Xsig_Xsig(&accumulator, &argSqrd);
   
         shr_Xsig(&accumulator, 1);
   
         accumulator.lsw |= 1;  /* A zero accumulator here would cause problems */
         negate_Xsig(&accumulator);
   
         /* The basic computation is complete. Now fix the answer to
            compensate for the error due to the approximation used for
            pi/2
            */
   
         /* This has an exponent of -65 */
         fix_up = 0x898cc517;
         /* The fix-up needs to be improved for larger args */
         if ( argSqrd.msw & 0xffc00000 )
           {
             /* Get about 32 bit precision in these: */
             fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
           }
         fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
   
         adj = accumulator.lsw;    /* temp save */
         accumulator.lsw -= fix_up;
         if ( accumulator.lsw > adj )
           XSIG_LL(accumulator) --;
   
         echange = round_Xsig(&accumulator);
   
         setexponentpos(&result, echange - 1);
       }
   
     significand(&result) = XSIG_LL(accumulator);
     setsign(&result, getsign(st0_ptr));
     FPU_copy_to_reg0(&result, TAG_Valid);
   
   #ifdef PARANOID
     if ( (exponent(&result) >= 0)
         && (significand(&result) > 0x8000000000000000LL) )
       {
         EXCEPTION(EX_INTERNAL|0x150);
       }
   #endif /* PARANOID */
   
   }
   
   
   
   /*--- poly_cos() ------------------------------------------------------------+
    |                                                                           |
    +---------------------------------------------------------------------------*/
   void    poly_cos(FPU_REG *st0_ptr)
   {
     FPU_REG             result;
     long int            exponent, exp2, echange;
     Xsig                accumulator, argSqrd, fix_up, argTo4;
     unsigned long long  fixed_arg;
   
   #ifdef PARANOID
     if ( (exponent(st0_ptr) > 0)
         || ((exponent(st0_ptr) == 0)
             && (significand(st0_ptr) > 0xc90fdaa22168c234LL)) )
       {
         EXCEPTION(EX_Invalid);
         FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
         return;
       }
   #endif /* PARANOID */
   
     exponent = exponent(st0_ptr);
   
     accumulator.lsw = accumulator.midw = accumulator.msw = 0;
   
     if ( (exponent < -1) || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54)) )
       {
         /* arg is < 0.687705 */
   
         argSqrd.msw = st0_ptr->sigh; argSqrd.midw = st0_ptr->sigl;
         argSqrd.lsw = 0;
         mul64_Xsig(&argSqrd, &significand(st0_ptr));
   
         if ( exponent < -1 )
           {
             /* shift the argument right by the required places */
             shr_Xsig(&argSqrd, 2*(-1-exponent));
           }
   
         argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
         argTo4.lsw = argSqrd.lsw;
         mul_Xsig_Xsig(&argTo4, &argTo4);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
                         N_COEFF_NH-1);
         mul_Xsig_Xsig(&accumulator, &argSqrd);
         negate_Xsig(&accumulator);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
                         N_COEFF_PH-1);
         negate_Xsig(&accumulator);
   
         mul64_Xsig(&accumulator, &significand(st0_ptr));
         mul64_Xsig(&accumulator, &significand(st0_ptr));
         shr_Xsig(&accumulator, -2*(1+exponent));
   
         shr_Xsig(&accumulator, 3);
         negate_Xsig(&accumulator);
   
         add_Xsig_Xsig(&accumulator, &argSqrd);
   
         shr_Xsig(&accumulator, 1);
   
         /* It doesn't matter if accumulator is all zero here, the
            following code will work ok */
         negate_Xsig(&accumulator);
   
         if ( accumulator.lsw & 0x80000000 )
           XSIG_LL(accumulator) ++;
         if ( accumulator.msw == 0 )
           {
             /* The result is 1.0 */
             FPU_copy_to_reg0(&CONST_1, TAG_Valid);
             return;
           }
         else
           {
             significand(&result) = XSIG_LL(accumulator);
         
             /* will be a valid positive nr with expon = -1 */
             setexponentpos(&result, -1);
           }
       }
     else
       {
         fixed_arg = significand(st0_ptr);
   
         if ( exponent == 0 )
           {
             /* The argument is >= 1.0 */
   
             /* Put the binary point at the left. */
             fixed_arg <<= 1;
           }
         /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
         fixed_arg = 0x921fb54442d18469LL - fixed_arg;
         /* There is a special case which arises due to rounding, to fix here. */
         if ( fixed_arg == 0xffffffffffffffffLL )
           fixed_arg = 0;
   
         exponent = -1;
         exp2 = -1;
   
         /* A shift is needed here only for a narrow range of arguments,
            i.e. for fixed_arg approx 2^-32, but we pick up more... */
         if ( !(LL_MSW(fixed_arg) & 0xffff0000) )
           {
             fixed_arg <<= 16;
             exponent -= 16;
             exp2 -= 16;
           }
   
         XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
         mul64_Xsig(&argSqrd, &fixed_arg);
   
         if ( exponent < -1 )
           {
             /* shift the argument right by the required places */
             shr_Xsig(&argSqrd, 2*(-1-exponent));
           }
   
         argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
         argTo4.lsw = argSqrd.lsw;
         mul_Xsig_Xsig(&argTo4, &argTo4);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
                         N_COEFF_N-1);
         mul_Xsig_Xsig(&accumulator, &argSqrd);
         negate_Xsig(&accumulator);
   
         polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
                         N_COEFF_P-1);
   
         shr_Xsig(&accumulator, 2);    /* Divide by four */
         accumulator.msw |= 0x80000000;  /* Add 1.0 */
   
         mul64_Xsig(&accumulator, &fixed_arg);
         mul64_Xsig(&accumulator, &fixed_arg);
         mul64_Xsig(&accumulator, &fixed_arg);
   
         /* Divide by four, FPU_REG compatible, etc */
         exponent = 3*exponent;
   
         /* The minimum exponent difference is 3 */
         shr_Xsig(&accumulator, exp2 - exponent);
   
         negate_Xsig(&accumulator);
         XSIG_LL(accumulator) += fixed_arg;
   
         /* The basic computation is complete. Now fix the answer to
            compensate for the error due to the approximation used for
            pi/2
            */
   
         /* This has an exponent of -65 */
         XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
         fix_up.lsw = 0;
   
         /* The fix-up needs to be improved for larger args */
         if ( argSqrd.msw & 0xffc00000 )
           {
             /* Get about 32 bit precision in these: */
             fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
             fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
           }
   
         exp2 += norm_Xsig(&accumulator);
         shr_Xsig(&accumulator, 1); /* Prevent overflow */
         exp2++;
         shr_Xsig(&fix_up, 65 + exp2);
   
         add_Xsig_Xsig(&accumulator, &fix_up);
   
         echange = round_Xsig(&accumulator);
   
         setexponentpos(&result, exp2 + echange);
         significand(&result) = XSIG_LL(accumulator);
       }
   
     FPU_copy_to_reg0(&result, TAG_Valid);
   
   #ifdef PARANOID
     if ( (exponent(&result) >= 0)
         && (significand(&result) > 0x8000000000000000LL) )
       {
         EXCEPTION(EX_INTERNAL|0x151);
       }
   #endif /* PARANOID */
   
   }

Cache object: bcb24bb77e3efefeaeabf4ccf6212492