poly_l2.c
上传用户:lgb322
上传日期:2013-02-24
资源大小:30529k
文件大小:8k
源码类别:

嵌入式Linux

开发平台:

Unix_Linux

  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,1997                                         |
  7.  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
  8.  |                  E-mail   billm@suburbia.net                              |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12. #include "exception.h"
  13. #include "reg_constant.h"
  14. #include "fpu_emu.h"
  15. #include "fpu_system.h"
  16. #include "control_w.h"
  17. #include "poly.h"
  18. static void log2_kernel(FPU_REG const *arg, u_char argsign,
  19. Xsig *accum_result, long int *expon);
  20. /*--- poly_l2() -------------------------------------------------------------+
  21.  |   Base 2 logarithm by a polynomial approximation.                         |
  22.  +---------------------------------------------------------------------------*/
  23. void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
  24. {
  25.   long int        exponent, expon, expon_expon;
  26.   Xsig                 accumulator, expon_accum, yaccum;
  27.   u_char        sign, argsign;
  28.   FPU_REG              x;
  29.   int                  tag;
  30.   exponent = exponent16(st0_ptr);
  31.   /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
  32.   if ( st0_ptr->sigh > (unsigned)0xb504f334 )
  33.     {
  34.       /* Treat as  sqrt(2)/2 < st0_ptr < 1 */
  35.       significand(&x) = - significand(st0_ptr);
  36.       setexponent16(&x, -1);
  37.       exponent++;
  38.       argsign = SIGN_NEG;
  39.     }
  40.   else
  41.     {
  42.       /* Treat as  1 <= st0_ptr < sqrt(2) */
  43.       x.sigh = st0_ptr->sigh - 0x80000000;
  44.       x.sigl = st0_ptr->sigl;
  45.       setexponent16(&x, 0);
  46.       argsign = SIGN_POS;
  47.     }
  48.   tag = FPU_normalize_nuo(&x);
  49.   if ( tag == TAG_Zero )
  50.     {
  51.       expon = 0;
  52.       accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  53.     }
  54.   else
  55.     {
  56.       log2_kernel(&x, argsign, &accumulator, &expon);
  57.     }
  58.   if ( exponent < 0 )
  59.     {
  60.       sign = SIGN_NEG;
  61.       exponent = -exponent;
  62.     }
  63.   else
  64.     sign = SIGN_POS;
  65.   expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
  66.   if ( exponent )
  67.     {
  68.       expon_expon = 31 + norm_Xsig(&expon_accum);
  69.       shr_Xsig(&accumulator, expon_expon - expon);
  70.       if ( sign ^ argsign )
  71. negate_Xsig(&accumulator);
  72.       add_Xsig_Xsig(&accumulator, &expon_accum);
  73.     }
  74.   else
  75.     {
  76.       expon_expon = expon;
  77.       sign = argsign;
  78.     }
  79.   yaccum.lsw = 0; XSIG_LL(yaccum) = significand(st1_ptr);
  80.   mul_Xsig_Xsig(&accumulator, &yaccum);
  81.   expon_expon += round_Xsig(&accumulator);
  82.   if ( accumulator.msw == 0 )
  83.     {
  84.       FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
  85.       return;
  86.     }
  87.   significand(st1_ptr) = XSIG_LL(accumulator);
  88.   setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
  89.   tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
  90.   FPU_settagi(1, tag);
  91.   set_precision_flag_up();  /* 80486 appears to always do this */
  92.   return;
  93. }
  94. /*--- poly_l2p1() -----------------------------------------------------------+
  95.  |   Base 2 logarithm by a polynomial approximation.                         |
  96.  |   log2(x+1)                                                               |
  97.  +---------------------------------------------------------------------------*/
  98. int poly_l2p1(u_char sign0, u_char sign1,
  99.   FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest)
  100. {
  101.   u_char              tag;
  102.   long int         exponent;
  103.   Xsig               accumulator, yaccum;
  104.   if ( exponent16(st0_ptr) < 0 )
  105.     {
  106.       log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
  107.       yaccum.lsw = 0;
  108.       XSIG_LL(yaccum) = significand(st1_ptr);
  109.       mul_Xsig_Xsig(&accumulator, &yaccum);
  110.       exponent += round_Xsig(&accumulator);
  111.       exponent += exponent16(st1_ptr) + 1;
  112.       if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER;
  113.       significand(dest) = XSIG_LL(accumulator);
  114.       setexponent16(dest, exponent);
  115.       tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
  116.       FPU_settagi(1, tag);
  117.       if ( tag == TAG_Valid )
  118. set_precision_flag_up();   /* 80486 appears to always do this */
  119.     }
  120.   else
  121.     {
  122.       /* The magnitude of st0_ptr is far too large. */
  123.       if ( sign0 != SIGN_POS )
  124. {
  125.   /* Trying to get the log of a negative number. */
  126. #ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
  127.   changesign(st1_ptr);
  128. #else
  129.   if ( arith_invalid(1) < 0 )
  130.     return 1;
  131. #endif /* PECULIAR_486 */
  132. }
  133.       /* 80486 appears to do this */
  134.       if ( sign0 == SIGN_NEG )
  135. set_precision_flag_down();
  136.       else
  137. set_precision_flag_up();
  138.     }
  139.   if ( exponent(dest) <= EXP_UNDER )
  140.     EXCEPTION(EX_Underflow);
  141.   return 0;
  142. }
  143. #undef HIPOWER
  144. #define HIPOWER 10
  145. static const unsigned long long logterms[HIPOWER] =
  146. {
  147.   0x2a8eca5705fc2ef0LL,
  148.   0xf6384ee1d01febceLL,
  149.   0x093bb62877cdf642LL,
  150.   0x006985d8a9ec439bLL,
  151.   0x0005212c4f55a9c8LL,
  152.   0x00004326a16927f0LL,
  153.   0x0000038d1d80a0e7LL,
  154.   0x0000003141cc80c6LL,
  155.   0x00000002b1668c9fLL,
  156.   0x000000002c7a46aaLL
  157. };
  158. static const unsigned long leadterm = 0xb8000000;
  159. /*--- log2_kernel() ---------------------------------------------------------+
  160.  |   Base 2 logarithm by a polynomial approximation.                         |
  161.  |   log2(x+1)                                                               |
  162.  +---------------------------------------------------------------------------*/
  163. static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
  164. long int *expon)
  165. {
  166.   long int             exponent, adj;
  167.   unsigned long long   Xsq;
  168.   Xsig                 accumulator, Numer, Denom, argSignif, arg_signif;
  169.   exponent = exponent16(arg);
  170.   Numer.lsw = Denom.lsw = 0;
  171.   XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
  172.   if ( argsign == SIGN_POS )
  173.     {
  174.       shr_Xsig(&Denom, 2 - (1 + exponent));
  175.       Denom.msw |= 0x80000000;
  176.       div_Xsig(&Numer, &Denom, &argSignif);
  177.     }
  178.   else
  179.     {
  180.       shr_Xsig(&Denom, 1 - (1 + exponent));
  181.       negate_Xsig(&Denom);
  182.       if ( Denom.msw & 0x80000000 )
  183. {
  184.   div_Xsig(&Numer, &Denom, &argSignif);
  185.   exponent ++;
  186. }
  187.       else
  188. {
  189.   /* Denom must be 1.0 */
  190.   argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
  191.   argSignif.msw = Numer.msw;
  192. }
  193.     }
  194. #ifndef PECULIAR_486
  195.   /* Should check here that  |local_arg|  is within the valid range */
  196.   if ( exponent >= -2 )
  197.     {
  198.       if ( (exponent > -2) ||
  199.   (argSignif.msw > (unsigned)0xafb0ccc0) )
  200. {
  201.   /* The argument is too large */
  202. }
  203.     }
  204. #endif /* PECULIAR_486 */
  205.   arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif);
  206.   adj = norm_Xsig(&argSignif);
  207.   accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif);
  208.   mul_Xsig_Xsig(&accumulator, &accumulator);
  209.   shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj)));
  210.   Xsq = XSIG_LL(accumulator);
  211.   if ( accumulator.lsw & 0x80000000 )
  212.     Xsq++;
  213.   accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  214.   /* Do the basic fixed point polynomial evaluation */
  215.   polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1);
  216.   mul_Xsig_Xsig(&accumulator, &argSignif);
  217.   shr_Xsig(&accumulator, 6 - adj);
  218.   mul32_Xsig(&arg_signif, leadterm);
  219.   add_two_Xsig(&accumulator, &arg_signif, &exponent);
  220.   *expon = exponent + 1;
  221.   accum_result->lsw = accumulator.lsw;
  222.   accum_result->midw = accumulator.midw;
  223.   accum_result->msw = accumulator.msw;
  224. }