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_l2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-01  |  6.9 KB  |  256 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  poly_l2.c                                                                |
  3.  |                                                                           |
  4.  | Compute the base 2 log 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.  
  14. #include "exception.h"
  15. #include "reg_constant.h"
  16. #include "fpu_emu.h"
  17. #include "control_w.h"
  18. #include "poly.h"
  19.  
  20.  
  21.  
  22. static void log2_kernel(FPU_REG const *arg,
  23.             Xsig *accum_result, long int *expon);
  24.  
  25.  
  26. /*--- poly_l2() -------------------------------------------------------------+
  27.  |   Base 2 logarithm by a polynomial approximation.                         |
  28.  +---------------------------------------------------------------------------*/
  29. void    poly_l2(FPU_REG const *arg, FPU_REG const *y, FPU_REG *result)
  30. {
  31.   long int           exponent, expon, expon_expon;
  32.   Xsig                 accumulator, expon_accum, yaccum;
  33.   char               sign;
  34.   FPU_REG              x;
  35.  
  36.  
  37.   exponent = arg->exp - EXP_BIAS;
  38.  
  39.   /* From arg, make a number > sqrt(2)/2 and < sqrt(2) */
  40.   if ( arg->sigh > (unsigned)0xb504f334 )
  41.     {
  42.       /* Treat as  sqrt(2)/2 < arg < 1 */
  43.       significand(&x) = - significand(arg);
  44.       x.sign = SIGN_NEG;
  45.       x.tag = TW_Valid;
  46.       x.exp = EXP_BIAS-1;
  47.       exponent++;
  48.       normalize(&x);
  49.     }
  50.   else
  51.     {
  52.       /* Treat as  1 <= arg < sqrt(2) */
  53.       x.sigh = arg->sigh - 0x80000000;
  54.       x.sigl = arg->sigl;
  55.       x.sign = SIGN_POS;
  56.       x.tag = TW_Valid;
  57.       x.exp = EXP_BIAS;
  58.       normalize(&x);
  59.     }
  60.  
  61.   if ( x.tag == TW_Zero )
  62.     {
  63.       expon = 0;
  64.       accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  65.     }
  66.   else
  67.     {
  68.       log2_kernel(&x, &accumulator, &expon);
  69.     }
  70.  
  71.   sign = exponent < 0;
  72.   if ( sign ) exponent = -exponent;
  73.   expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
  74.   if ( exponent )
  75.     {
  76.       expon_expon = 31 + norm_Xsig(&expon_accum);
  77.       shr_Xsig(&accumulator, expon_expon - expon);
  78.  
  79.       if ( sign ^ (x.sign == SIGN_NEG) )
  80.     negate_Xsig(&accumulator);
  81.       add_Xsig_Xsig(&accumulator, &expon_accum);
  82.     }
  83.   else
  84.     {
  85.       expon_expon = expon;
  86.       sign = x.sign;
  87.     }
  88.  
  89.   yaccum.lsw = 0; XSIG_LL(yaccum) = significand(y);
  90.   mul_Xsig_Xsig(&accumulator, &yaccum);
  91.  
  92.   expon_expon += round_Xsig(&accumulator);
  93.  
  94.   if ( accumulator.msw == 0 )
  95.     {
  96.       reg_move(&CONST_Z, y);
  97.     }
  98.   else
  99.     {
  100.       result->exp = expon_expon + y->exp + 1;
  101.       significand(result) = XSIG_LL(accumulator);
  102.       result->tag = TW_Valid; /* set the tags to Valid */
  103.       result->sign = sign ^ y->sign;
  104.     }
  105.  
  106.   return;
  107. }
  108.  
  109.  
  110. /*--- poly_l2p1() -----------------------------------------------------------+
  111.  |   Base 2 logarithm by a polynomial approximation.                         |
  112.  |   log2(x+1)                                                               |
  113.  +---------------------------------------------------------------------------*/
  114. int    poly_l2p1(FPU_REG const *arg, FPU_REG const *y, FPU_REG *result)
  115. {
  116.   char                 sign;
  117.   long int             exponent;
  118.   Xsig                 accumulator, yaccum;
  119.  
  120.  
  121.   sign = arg->sign;
  122.  
  123.   if ( arg->exp < EXP_BIAS )
  124.     {
  125.       log2_kernel(arg, &accumulator, &exponent);
  126.  
  127.       yaccum.lsw = 0;
  128.       XSIG_LL(yaccum) = significand(y);
  129.       mul_Xsig_Xsig(&accumulator, &yaccum);
  130.  
  131.       exponent += round_Xsig(&accumulator);
  132.  
  133.       result->exp = exponent + y->exp + 1;
  134.       significand(result) = XSIG_LL(accumulator);
  135.       result->tag = TW_Valid; /* set the tags to Valid */
  136.       result->sign = sign ^ y->sign;
  137.  
  138.       return 0;
  139.     }
  140.   else
  141.     {
  142.       /* The magnitude of arg is far too large. */
  143.       reg_move(y, result);
  144.       if ( sign != SIGN_POS )
  145.     {
  146.       /* Trying to get the log of a negative number. */
  147.       return 1;
  148.     }
  149.       else
  150.     {
  151.       return 0;
  152.     }
  153.     }
  154.  
  155. }
  156.  
  157.  
  158.  
  159.  
  160. #undef HIPOWER
  161. #define    HIPOWER    10
  162. static const unsigned long long logterms[HIPOWER] =
  163. {
  164.   0x2a8eca5705fc2ef0LL,
  165.   0xf6384ee1d01febceLL,
  166.   0x093bb62877cdf642LL,
  167.   0x006985d8a9ec439bLL,
  168.   0x0005212c4f55a9c8LL,
  169.   0x00004326a16927f0LL,
  170.   0x0000038d1d80a0e7LL,
  171.   0x0000003141cc80c6LL,
  172.   0x00000002b1668c9fLL,
  173.   0x000000002c7a46aaLL
  174. };
  175.  
  176. static const unsigned long leadterm = 0xb8000000;
  177.  
  178.  
  179. /*--- log2_kernel() ---------------------------------------------------------+
  180.  |   Base 2 logarithm by a polynomial approximation.                         |
  181.  |   log2(x+1)                                                               |
  182.  +---------------------------------------------------------------------------*/
  183. static void log2_kernel(FPU_REG const *arg, Xsig *accum_result,
  184.             long int *expon)
  185. {
  186.   char                 sign;
  187.   long int             exponent, adj;
  188.   unsigned long long   Xsq;
  189.   Xsig                 accumulator, Numer, Denom, argSignif, arg_signif;
  190.  
  191.   sign = arg->sign;
  192.  
  193.   exponent = arg->exp - EXP_BIAS;
  194.   Numer.lsw = Denom.lsw = 0;
  195.   XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
  196.   if ( sign == SIGN_POS )
  197.     {
  198.       shr_Xsig(&Denom, 2 - (1 + exponent));
  199.       Denom.msw |= 0x80000000;
  200.       div_Xsig(&Numer, &Denom, &argSignif);
  201.     }
  202.   else
  203.     {
  204.       shr_Xsig(&Denom, 1 - (1 + exponent));
  205.       negate_Xsig(&Denom);
  206.       if ( Denom.msw & 0x80000000 )
  207.     {
  208.       div_Xsig(&Numer, &Denom, &argSignif);
  209.       exponent ++;
  210.     }
  211.       else
  212.     {
  213.       /* Denom must be 1.0 */
  214.       argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
  215.       argSignif.msw = Numer.msw;
  216.     }
  217.     }
  218.  
  219. #ifndef PECULIAR_486
  220.   /* Should check here that  |local_arg|  is within the valid range */
  221.   if ( exponent >= -2 )
  222.     {
  223.       if ( (exponent > -2) ||
  224.       (argSignif.msw > (unsigned)0xafb0ccc0) )
  225.     {
  226.       /* The argument is too large */
  227.     }
  228.     }
  229. #endif PECULIAR_486
  230.  
  231.   arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif);
  232.   adj = norm_Xsig(&argSignif);
  233.   accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif);
  234.   mul_Xsig_Xsig(&accumulator, &accumulator);
  235.   shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj)));
  236.   Xsq = XSIG_LL(accumulator);
  237.   if ( accumulator.lsw & 0x80000000 )
  238.     Xsq++;
  239.  
  240.   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  241.   /* Do the basic fixed point polynomial evaluation */
  242.   polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1);
  243.  
  244.   mul_Xsig_Xsig(&accumulator, &argSignif);
  245.   shr_Xsig(&accumulator, 6 - adj);
  246.  
  247.   mul32_Xsig(&arg_signif, leadterm);
  248.   add_two_Xsig(&accumulator, &arg_signif, &exponent);
  249.  
  250.   *expon = exponent + 1;
  251.   accum_result->lsw = accumulator.lsw;
  252.   accum_result->midw = accumulator.midw;
  253.   accum_result->msw = accumulator.msw;
  254.  
  255. }
  256.