home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / KERNEL-S / V1.2 / LINUX-1.2 / LINUX-1 / linux / arch / i386 / math-emu / poly_sin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-01  |  11.2 KB  |  409 lines
  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. }
  409.