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_tan.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-19  |  6.7 KB  |  214 lines

  1. /*---------------------------------------------------------------------------+
  2.  |  poly_tan.c                                                               |
  3.  |                                                                           |
  4.  | Compute the tan 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 "control_w.h"
  17. #include "poly.h"
  18.  
  19.  
  20. #define    HiPOWERop    3    /* odd poly, positive terms */
  21. static const unsigned long long oddplterm[HiPOWERop] =
  22. {
  23.   0x0000000000000000LL,
  24.   0x0051a1cf08fca228LL,
  25.   0x0000000071284ff7LL
  26. };
  27.  
  28. #define    HiPOWERon    2    /* odd poly, negative terms */
  29. static const unsigned long long oddnegterm[HiPOWERon] =
  30. {
  31.    0x1291a9a184244e80LL,
  32.    0x0000583245819c21LL
  33. };
  34.  
  35. #define    HiPOWERep    2    /* even poly, positive terms */
  36. static const unsigned long long evenplterm[HiPOWERep] =
  37. {
  38.   0x0e848884b539e888LL,
  39.   0x00003c7f18b887daLL
  40. };
  41.  
  42. #define    HiPOWERen    2    /* even poly, negative terms */
  43. static const unsigned long long evennegterm[HiPOWERen] =
  44. {
  45.   0xf1f0200fd51569ccLL,
  46.   0x003afb46105c4432LL
  47. };
  48.  
  49. static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
  50.  
  51.  
  52. /*--- poly_tan() ------------------------------------------------------------+
  53.  |                                                                           |
  54.  +---------------------------------------------------------------------------*/
  55. void    poly_tan(FPU_REG const *arg, FPU_REG *result)
  56. {
  57.   long int            exponent;
  58.   int                   invert;
  59.   Xsig                  argSq, argSqSq, accumulatoro, accumulatore, accum,
  60.                         argSignif, fix_up;
  61.   unsigned long         adj;
  62.  
  63.   exponent = arg->exp - EXP_BIAS;
  64.  
  65. #ifdef PARANOID
  66.   if ( arg->sign != 0 )    /* Can't hack a number < 0.0 */
  67.     { arith_invalid(result); return; }  /* Need a positive number */
  68. #endif PARANOID
  69.  
  70.   /* Split the problem into two domains, smaller and larger than pi/4 */
  71.   if ( (exponent == 0) || ((exponent == -1) && (arg->sigh > 0xc90fdaa2)) )
  72.     {
  73.       /* The argument is greater than (approx) pi/4 */
  74.       invert = 1;
  75.       accum.lsw = 0;
  76.       XSIG_LL(accum) = significand(arg);
  77.  
  78.       if ( exponent == 0 )
  79.     {
  80.       /* The argument is >= 1.0 */
  81.       /* Put the binary point at the left. */
  82.       XSIG_LL(accum) <<= 1;
  83.     }
  84.       /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
  85.       XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
  86.  
  87.       argSignif.lsw = accum.lsw;
  88.       XSIG_LL(argSignif) = XSIG_LL(accum);
  89.       exponent = -1 + norm_Xsig(&argSignif);
  90.     }
  91.   else
  92.     {
  93.       invert = 0;
  94.       argSignif.lsw = 0;
  95.       XSIG_LL(accum) = XSIG_LL(argSignif) = significand(arg);
  96.  
  97.       if ( exponent < -1 )
  98.     {
  99.       /* shift the argument right by the required places */
  100.       if ( shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
  101.         XSIG_LL(accum) ++;    /* round up */
  102.     }
  103.     }
  104.  
  105.   XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
  106.   mul_Xsig_Xsig(&argSq, &argSq);
  107.   XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
  108.   mul_Xsig_Xsig(&argSqSq, &argSqSq);
  109.  
  110.   /* Compute the negative terms for the numerator polynomial */
  111.   accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
  112.   polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
  113.   mul_Xsig_Xsig(&accumulatoro, &argSq);
  114.   negate_Xsig(&accumulatoro);
  115.   /* Add the positive terms */
  116.   polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
  117.  
  118.   
  119.   /* Compute the positive terms for the denominator polynomial */
  120.   accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
  121.   polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
  122.   mul_Xsig_Xsig(&accumulatore, &argSq);
  123.   negate_Xsig(&accumulatore);
  124.   /* Add the negative terms */
  125.   polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
  126.   /* Multiply by arg^2 */
  127.   mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
  128.   mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
  129.   /* de-normalize and divide by 2 */
  130.   shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
  131.   negate_Xsig(&accumulatore);      /* This does 1 - accumulator */
  132.  
  133.   /* Now find the ratio. */
  134.   if ( accumulatore.msw == 0 )
  135.     {
  136.       /* accumulatoro must contain 1.0 here, (actually, 0) but it
  137.      really doesn't matter what value we use because it will
  138.      have negligible effect in later calculations
  139.      */
  140.       XSIG_LL(accum) = 0x8000000000000000LL;
  141.       accum.lsw = 0;
  142.     }
  143.   else
  144.     {
  145.       div_Xsig(&accumulatoro, &accumulatore, &accum);
  146.     }
  147.  
  148.   /* Multiply by 1/3 * arg^3 */
  149.   mul64_Xsig(&accum, &XSIG_LL(argSignif));
  150.   mul64_Xsig(&accum, &XSIG_LL(argSignif));
  151.   mul64_Xsig(&accum, &XSIG_LL(argSignif));
  152.   mul64_Xsig(&accum, &twothirds);
  153.   shr_Xsig(&accum, -2*(exponent+1));
  154.  
  155.   /* tan(arg) = arg + accum */
  156.   add_two_Xsig(&accum, &argSignif, &exponent);
  157.  
  158.   if ( invert )
  159.     {
  160.       /* We now have the value of tan(pi_2 - arg) where pi_2 is an
  161.      approximation for pi/2
  162.      */
  163.       /* The next step is to fix the answer to compensate for the
  164.      error due to the approximation used for pi/2
  165.      */
  166.  
  167.       /* This is (approx) delta, the error in our approx for pi/2
  168.      (see above). It has an exponent of -65
  169.      */
  170.       XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
  171.       fix_up.lsw = 0;
  172.  
  173.       if ( exponent == 0 )
  174.     adj = 0xffffffff;   /* We want approx 1.0 here, but
  175.                    this is close enough. */
  176.       else if ( exponent > -30 )
  177.     {
  178.       adj = accum.msw >> -(exponent+1);      /* tan */
  179.       mul_32_32(adj, adj, &adj);           /* tan^2 */
  180.     }
  181.       else
  182.     adj = 0;
  183.       mul_32_32(0x898cc517, adj, &adj);        /* delta * tan^2 */
  184.  
  185.       fix_up.msw += adj;
  186.       if ( !(fix_up.msw & 0x80000000) )   /* did fix_up overflow ? */
  187.     {
  188.       /* Yes, we need to add an msb */
  189.       shr_Xsig(&fix_up, 1);
  190.       fix_up.msw |= 0x80000000;
  191.       shr_Xsig(&fix_up, 64 + exponent);
  192.     }
  193.       else
  194.     shr_Xsig(&fix_up, 65 + exponent);
  195.  
  196.       add_two_Xsig(&accum, &fix_up, &exponent);
  197.  
  198.       /* accum now contains tan(pi/2 - arg).
  199.      Use tan(arg) = 1.0 / tan(pi/2 - arg)
  200.      */
  201.       accumulatoro.lsw = accumulatoro.midw = 0;
  202.       accumulatoro.msw = 0x80000000;
  203.       div_Xsig(&accumulatoro, &accum, &accum);
  204.       exponent = - exponent - 1;
  205.     }
  206.  
  207.   /* Transfer the result */
  208.   round_Xsig(&accum);
  209.   *(short *)&(result->sign) = 0;
  210.   significand(result) = XSIG_LL(accum);
  211.   result->exp = EXP_BIAS + exponent;
  212.  
  213. }
  214.