Return to poly_l2.c CVS log

Up to [MW Coherent from dump] / coherent / b / kernel / emulator

Annotation of coherent/b/kernel/emulator/poly_l2.c, revision 1.1

1.1     ! root        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    W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
        !             7:  |                       Australia.  E-mail [email protected]    |
        !             8:  |                                                                           |
        !             9:  |                                                                           |
        !            10:  +---------------------------------------------------------------------------*/
        !            11: 
        !            12: #include "fpu_system.h"
        !            13: #include "exception.h"
        !            14: #include "reg_constant.h"
        !            15: #include "fpu_emu.h"
        !            16: 
        !            17: 
        !            18: 
        !            19: #define        HIPOWER 9
        !            20: static unsigned short  lterms[HIPOWER][4] =
        !            21:        {
        !            22:        /* Ideal computation with these coeffs gives about
        !            23:           64.6 bit rel accuracy. */
        !            24:        { 0xe177, 0xb82f, 0x7652, 0x7154 },
        !            25:        { 0xee0f, 0xe80f, 0x2770, 0x7b1c },
        !            26:        { 0x0fc0, 0xbe87, 0xb143, 0x49dd },
        !            27:        { 0x78b9, 0xdadd, 0xec54, 0x34c2 },
        !            28:        { 0x003a, 0x5de9, 0x628b, 0x2909 },
        !            29:        { 0x5588, 0xed16, 0x4abf, 0x2193 },
        !            30:        { 0xb461, 0x85f7, 0x347a, 0x1c6a },
        !            31:        { 0x0975, 0x87b3, 0xd5bf, 0x1876 },
        !            32:        { 0xe85c, 0xcec9, 0x84e7, 0x187d }
        !            33:        };
        !            34: 
        !            35: 
        !            36: 
        !            37: 
        !            38: /*--- poly_l2() -------------------------------------------------------------+
        !            39:  |   Base 2 logarithm by a polynomial approximation.                         |
        !            40:  +---------------------------------------------------------------------------*/
        !            41: void   poly_l2(FPU_REG *arg, FPU_REG *result)
        !            42: {
        !            43:   short                  exponent;
        !            44:   char           zero;         /* flag for an Xx == 0 */
        !            45:   unsigned short  bits, shift;
        !            46:   long long       Xsq;
        !            47:   FPU_REG        accum, denom, num, Xx;
        !            48: 
        !            49: 
        !            50:   exponent = arg->exp - EXP_BIAS;
        !            51: 
        !            52:   accum.tag = TW_Valid;        /* set the tags to Valid */
        !            53: 
        !            54:   if ( arg->sigh > (unsigned)0xb504f334 )
        !            55:     {
        !            56:       /* This is good enough for the computation of the polynomial
        !            57:         sum, but actually results in a loss of precision for
        !            58:         the computation of Xx. This will matter only if exponent
        !            59:         becomes zero. */
        !            60:       exponent++;
        !            61:       accum.sign = 1;  /* sign to negative */
        !            62:       num.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
        !            63:       reg_u_div((long long *)&(CONST_1.sigl), (long long *)&(arg->sigl), &num);
        !            64:     }
        !            65:   else
        !            66:     {
        !            67:       accum.sign = 0;  /* set the sign to positive */
        !            68:       num.sigl = arg->sigl;            /* copy the mantissa */
        !            69:       num.sigh = arg->sigh;
        !            70:     }
        !            71: 
        !            72:   /* shift num left, lose the ms bit */
        !            73:   num.sigh <<= 1;
        !            74:   if ( num.sigl & 0x80000000 ) num.sigh |= 1;
        !            75:   num.sigl <<= 1;
        !            76: 
        !            77:   denom.sigl = num.sigl;
        !            78:   denom.sigh = num.sigh;
        !            79:   poly_div4((long long *)&(denom.sigl));
        !            80:   denom.sigh += 0x80000000;                    /* set the msb */
        !            81:   Xx.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
        !            82:   reg_u_div((long long *)&num.sigl, (long long *)&(denom.sigl), &Xx);
        !            83: 
        !            84:   zero = !(Xx.sigh | Xx.sigl);
        !            85:   
        !            86:   mul64((long long *)&Xx.sigl, (long long *)&Xx.sigl, &Xsq);
        !            87:   poly_div16(&Xsq);
        !            88: 
        !            89:   accum.exp = -1;              /* exponent of accum */
        !            90: 
        !            91:   /* Do the basic fixed point polynomial evaluation */
        !            92:   polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);
        !            93: 
        !            94:   if ( !exponent )
        !            95:     {
        !            96:       /* If the exponent is zero, then we would lose precision by
        !            97:         sticking to fixed point computation here */
        !            98:       /* We need to re-compute Xx because of loss of precision. */
        !            99:       FPU_REG   lXx;
        !           100:       char     sign;
        !           101:       
        !           102:       sign = accum.sign;
        !           103:       accum.sign = 0;
        !           104: 
        !           105:       /* make accum compatible and normalize */
        !           106:       accum.exp = EXP_BIAS + accum.exp;
        !           107:       normalize(&accum);
        !           108: 
        !           109:       if ( zero )
        !           110:        {
        !           111:          reg_move(&CONST_Z, result);
        !           112:        }
        !           113:       else
        !           114:        {
        !           115:          /* we need to re-compute lXx to better accuracy */
        !           116:          num.tag = TW_Valid;           /* set the tags to Vaild */
        !           117:          num.sign = 0;         /* set the sign to positive */
        !           118:          num.exp = EXP_BIAS - 1;
        !           119:          if ( sign )
        !           120:            {
        !           121:              /* The argument is of the form 1-x */
        !           122:              /* Use  1-1/(1-x) = x/(1-x) */
        !           123:              *((long long *)&num.sigl) = - *((long long *)&(arg->sigl));
        !           124:              normalize(&num);
        !           125:              reg_div(&num, arg, &num);
        !           126:            }
        !           127:          else
        !           128:            {
        !           129:              normalize(&num);
        !           130:            }
        !           131:          
        !           132:          denom.tag = TW_Valid; /* set the tags to Valid */
        !           133:          denom.sign = SIGN_POS;        /* set the sign to positive */
        !           134:          denom.exp = EXP_BIAS;
        !           135:          
        !           136:          reg_div(&num, &denom, &lXx);
        !           137: 
        !           138:          reg_u_mul(&lXx, &accum, &accum);
        !           139:          accum.exp += - EXP_BIAS + 1;
        !           140: 
        !           141:          reg_u_add(&lXx, &accum, result);
        !           142:          
        !           143:          normalize(result);
        !           144:        }
        !           145:       result->sign = sign;
        !           146:       return;
        !           147:     }
        !           148: 
        !           149:   mul64((long long *)&accum.sigl,
        !           150:        (long long *)&Xx.sigl, (long long *)&accum.sigl);
        !           151: 
        !           152:   *((long long *)(&accum.sigl)) += *((long long *)(&Xx.sigl));
        !           153: 
        !           154:   if ( Xx.sigh > accum.sigh )
        !           155:     {
        !           156:       /* There was an overflow */
        !           157: 
        !           158:       poly_div2((long long *)&accum.sigl);
        !           159:       accum.sigh |= 0x80000000;
        !           160:       accum.exp++;
        !           161:     }
        !           162: 
        !           163:   /* When we add the exponent to the accum result later, we will
        !           164:      require that their signs are the same. Here we ensure that
        !           165:      this is so. */
        !           166:   if ( exponent && ((exponent < 0) ^ (accum.sign)) )
        !           167:     {
        !           168:       /* signs are different */
        !           169: 
        !           170:       accum.sign = !accum.sign;
        !           171: 
        !           172:       /* An exceptional case is when accum is zero */
        !           173:       if ( accum.sigl | accum.sigh )
        !           174:        {
        !           175:          /* find 1-accum */
        !           176:          /* Shift to get exponent == 0 */
        !           177:          if ( accum.exp < 0 )
        !           178:            {
        !           179:              poly_div2((long long *)&accum.sigl);
        !           180:              accum.exp++;
        !           181:            }
        !           182:          /* Just negate, but throw away the sign */
        !           183:          *((long long *)&(accum.sigl)) = - *((long long *)&(accum.sigl));
        !           184:          if ( exponent < 0 )
        !           185:            exponent++;
        !           186:          else
        !           187:            exponent--;
        !           188:        }
        !           189:     }
        !           190: 
        !           191:   shift = exponent >= 0 ? exponent : -exponent ;
        !           192:   bits = 0;
        !           193:   if ( shift )
        !           194:     {
        !           195:       if ( accum.exp )
        !           196:        {
        !           197:          accum.exp++;
        !           198:          poly_div2((long long *)&accum.sigl);
        !           199:        }
        !           200:       while ( shift )
        !           201:        {
        !           202:          poly_div2((long long *)&accum.sigl);
        !           203:          if ( shift & 1)
        !           204:            accum.sigh |= 0x80000000;
        !           205:          shift >>= 1;
        !           206:          bits++;
        !           207:        }
        !           208:     }
        !           209: 
        !           210:   /* Convert to 64 bit signed-compatible */
        !           211:   accum.exp += bits + EXP_BIAS - 1;
        !           212: 
        !           213:   reg_move(&accum, result);
        !           214:   normalize(result);
        !           215: 
        !           216:   return;
        !           217: }
        !           218: 
        !           219: 
        !           220: /*--- poly_l2p1() -----------------------------------------------------------+
        !           221:  |   Base 2 logarithm by a polynomial approximation.                         |
        !           222:  |   log2(x+1)                                                               |
        !           223:  +---------------------------------------------------------------------------*/
        !           224: int    poly_l2p1(FPU_REG *arg, FPU_REG *result)
        !           225: {
        !           226:   char         sign = 0;
        !           227:   long long     Xsq;
        !           228:   FPU_REG              arg_pl1, denom, accum, local_arg, poly_arg;
        !           229: 
        !           230: 
        !           231:   sign = arg->sign;
        !           232: 
        !           233:   reg_add(arg, &CONST_1, &arg_pl1);
        !           234: 
        !           235:   if ( (arg_pl1.sign) | (arg_pl1.tag) )
        !           236:     {                  /* We need a valid positive number! */
        !           237:       return 1;
        !           238:     }
        !           239: 
        !           240:   reg_add(&CONST_1, &arg_pl1, &denom);
        !           241:   reg_div(arg, &denom, &local_arg);
        !           242:   local_arg.sign = 0;  /* Make the sign positive */
        !           243: 
        !           244:   /* Now we need to check that  |local_arg| is less than
        !           245:      3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */
        !           246: 
        !           247:   if ( local_arg.exp >= EXP_BIAS - 3 )
        !           248:     {
        !           249:       if ( (local_arg.exp > EXP_BIAS - 3) ||
        !           250:          (local_arg.sigh > (unsigned)0xafb0ccc0) )
        !           251:        {
        !           252:          /* The argument is large */
        !           253:          poly_l2(&arg_pl1, result); return 0;
        !           254:        }
        !           255:     }
        !           256: 
        !           257:   /* Make a copy of local_arg */
        !           258:   reg_move(&local_arg, &poly_arg);
        !           259: 
        !           260:   /* Get poly_arg bits aligned as required */
        !           261:   shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
        !           262: 
        !           263:   mul64((long long *)&(poly_arg.sigl), (long long *)&(poly_arg.sigl), &Xsq);
        !           264:   poly_div16(&Xsq);
        !           265: 
        !           266:   /* Do the basic fixed point polynomial evaluation */
        !           267:   polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);
        !           268: 
        !           269:   accum.tag = TW_Valid;        /* set the tags to Valid */
        !           270:   accum.sign = SIGN_POS;       /* and make accum positive */
        !           271: 
        !           272:   /* make accum compatible and normalize */
        !           273:   accum.exp = EXP_BIAS - 1;
        !           274:   normalize(&accum);
        !           275: 
        !           276:   reg_u_mul(&local_arg, &accum, &accum);
        !           277:   accum.exp -= EXP_BIAS - 1;
        !           278: 
        !           279:   reg_u_add(&local_arg, &accum, result);
        !           280: 
        !           281:   /* Multiply the result by 2 */
        !           282:   result->exp++;
        !           283: 
        !           284:   result->sign = sign;
        !           285:   
        !           286:   return 0;
        !           287: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.