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_atan.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-01  |  5.9 KB  |  198 lines
  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                                              |
  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 "status_w.h"
  17. #include "control_w.h"
  18. #include "poly.h"
  19.  
  20.  
  21. #define    HIPOWERon    6    /* odd poly, negative terms */
  22. static const unsigned long long oddnegterms[HIPOWERon] =
  23. {
  24.   0x0000000000000000LL, /* Dummy (not for - 1.0) */
  25.   0x015328437f756467LL,
  26.   0x0005dda27b73dec6LL,
  27.   0x0000226bf2bfb91aLL,
  28.   0x000000ccc439c5f7LL,
  29.   0x0000000355438407LL
  30. } ;
  31.  
  32. #define    HIPOWERop    6    /* odd poly, positive terms */
  33. static const unsigned long long oddplterms[HIPOWERop] =
  34. {
  35. /*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
  36.   0x0db55a71875c9ac2LL,
  37.   0x0029fce2d67880b0LL,
  38.   0x0000dfd3908b4596LL,
  39.   0x00000550fd61dab4LL,
  40.   0x0000001c9422b3f9LL,
  41.   0x000000003e3301e1LL
  42. };
  43.  
  44. static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
  45.  
  46. static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
  47.  
  48. static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
  49.  
  50.  
  51. /*--- poly_atan() -----------------------------------------------------------+
  52.  |                                                                           |
  53.  +---------------------------------------------------------------------------*/
  54. void    poly_atan(FPU_REG *arg1, FPU_REG *arg2, FPU_REG *result)
  55. {
  56.   char                transformed, inverted,
  57.                         sign1 = arg1->sign, sign2 = arg2->sign;
  58.   long int           exponent, dummy_exp;
  59.   Xsig                  accumulator, Numer, Denom, accumulatore, argSignif,
  60.                         argSq, argSqSq;
  61.   
  62.  
  63.   arg1->sign = arg2->sign = SIGN_POS;
  64.   if ( (compare(arg2) & ~COMP_Denormal) == COMP_A_lt_B )
  65.     {
  66.       inverted = 1;
  67.       exponent = arg1->exp - arg2->exp;
  68.       Numer.lsw = Denom.lsw = 0;
  69.       XSIG_LL(Numer) = significand(arg1);
  70.       XSIG_LL(Denom) = significand(arg2);
  71.     }
  72.   else
  73.     {
  74.       inverted = 0;
  75.       exponent = arg2->exp - arg1->exp;
  76.       Numer.lsw = Denom.lsw = 0;
  77.       XSIG_LL(Numer) = significand(arg2);
  78.       XSIG_LL(Denom) = significand(arg1);
  79.      }
  80.   div_Xsig(&Numer, &Denom, &argSignif);
  81.   exponent += norm_Xsig(&argSignif);
  82.  
  83.   if ( (exponent >= -1)
  84.       || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
  85.     {
  86.       /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
  87.       /* Convert the argument by an identity for atan */
  88.       transformed = 1;
  89.  
  90.       if ( exponent >= 0 )
  91.     {
  92. #ifdef PARANOID
  93.       if ( !( (exponent == 0) && 
  94.          (argSignif.lsw == 0) && (argSignif.midw == 0) &&
  95.          (argSignif.msw == 0x80000000) ) )
  96.         {
  97.           EXCEPTION(EX_INTERNAL|0x104);  /* There must be a logic error */
  98.           return;
  99.         }
  100. #endif PARANOID
  101.       argSignif.msw = 0;   /* Make the transformed arg -> 0.0 */
  102.     }
  103.       else
  104.     {
  105.       Numer.lsw = Denom.lsw = argSignif.lsw;
  106.       XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
  107.  
  108.       if ( exponent < -1 )
  109.         shr_Xsig(&Numer, -1-exponent);
  110.       negate_Xsig(&Numer);
  111.       
  112.       shr_Xsig(&Denom, -exponent);
  113.       Denom.msw |= 0x80000000;
  114.       
  115.       div_Xsig(&Numer, &Denom, &argSignif);
  116.  
  117.       exponent = -1 + norm_Xsig(&argSignif);
  118.     }
  119.     }
  120.   else
  121.     {
  122.       transformed = 0;
  123.     }
  124.  
  125.   argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw;
  126.   argSq.msw = argSignif.msw;
  127.   mul_Xsig_Xsig(&argSq, &argSq);
  128.   
  129.   argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw;
  130.   mul_Xsig_Xsig(&argSqSq, &argSqSq);
  131.  
  132.   accumulatore.lsw = argSq.lsw;
  133.   XSIG_LL(accumulatore) = XSIG_LL(argSq);
  134.  
  135.   shr_Xsig(&argSq, 2*(-1-exponent-1));
  136.   shr_Xsig(&argSqSq, 4*(-1-exponent-1));
  137.  
  138.   /* Now have argSq etc with binary point at the left
  139.      .1xxxxxxxx */
  140.  
  141.   /* Do the basic fixed point polynomial evaluation */
  142.   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  143.   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
  144.            oddplterms, HIPOWERop-1);
  145.   mul64_Xsig(&accumulator, &XSIG_LL(argSq));
  146.   negate_Xsig(&accumulator);
  147.   polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
  148.   negate_Xsig(&accumulator);
  149.   add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
  150.  
  151.   mul64_Xsig(&accumulatore, &denomterm);
  152.   shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
  153.   accumulatore.msw |= 0x80000000;
  154.  
  155.   div_Xsig(&accumulator, &accumulatore, &accumulator);
  156.  
  157.   mul_Xsig_Xsig(&accumulator, &argSignif);
  158.   mul_Xsig_Xsig(&accumulator, &argSq);
  159.  
  160.   shr_Xsig(&accumulator, 3);
  161.   negate_Xsig(&accumulator);
  162.   add_Xsig_Xsig(&accumulator, &argSignif);
  163.  
  164.   if ( transformed )
  165.     {
  166.       /* compute pi/4 - accumulator */
  167.       shr_Xsig(&accumulator, -1-exponent);
  168.       negate_Xsig(&accumulator);
  169.       add_Xsig_Xsig(&accumulator, &pi_signif);
  170.       exponent = -1;
  171.     }
  172.  
  173.   if ( inverted )
  174.     {
  175.       /* compute pi/2 - accumulator */
  176.       shr_Xsig(&accumulator, -exponent);
  177.       negate_Xsig(&accumulator);
  178.       add_Xsig_Xsig(&accumulator, &pi_signif);
  179.       exponent = 0;
  180.     }
  181.  
  182.   if ( sign1 )
  183.     {
  184.       /* compute pi - accumulator */
  185.       shr_Xsig(&accumulator, 1 - exponent);
  186.       negate_Xsig(&accumulator);
  187.       add_Xsig_Xsig(&accumulator, &pi_signif);
  188.       exponent = 1;
  189.     }
  190.  
  191.   exponent += round_Xsig(&accumulator);
  192.   significand(result) = XSIG_LL(accumulator);
  193.   result->exp = exponent + EXP_BIAS;
  194.   result->tag = TW_Valid;
  195.   result->sign = sign2;
  196.  
  197. }
  198.