Return to poly_tan.c CVS log

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

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

1.1     ! root        1: /*---------------------------------------------------------------------------+
        !             2:  |  poly_tan.c                                                               |
        !             3:  |                                                                           |
        !             4:  | Compute the tan 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: #define        HIPOWERop       3       /* odd poly, positive terms */
        !            19: static unsigned short  oddplterms[HIPOWERop][4] =
        !            20:        {
        !            21:        { 0x846a, 0x42d1, 0xb544, 0x921f},
        !            22:        { 0x6fb2, 0x0215, 0x95c0, 0x099c},
        !            23:        { 0xfce6, 0x0cc8, 0x1c9a, 0x0000}
        !            24:        };
        !            25: 
        !            26: #define        HIPOWERon       2       /* odd poly, negative terms */
        !            27: static unsigned short  oddnegterms[HIPOWERon][4] =
        !            28:        {
        !            29:        { 0x6906, 0xe205, 0x25c8, 0x8838},
        !            30:        { 0x1dd7, 0x3fe3, 0x944e, 0x002c}
        !            31:        };
        !            32: 
        !            33: #define        HIPOWERep       2       /* even poly, positive terms */
        !            34: static unsigned short  evenplterms[HIPOWERep][4] =
        !            35:        {
        !            36:        { 0xdb8f, 0x3761, 0x1432, 0x2acf},
        !            37:        { 0x16eb, 0x13c1, 0x3099, 0x0003}
        !            38:        };
        !            39: 
        !            40: #define        HIPOWERen       2       /* even poly, negative terms */
        !            41: static unsigned short  evennegterms[HIPOWERen][4] =
        !            42:        {
        !            43:        { 0x3a7c, 0xe4c5, 0x7f87, 0x2945},
        !            44:        { 0x572b, 0x664c, 0xc543, 0x018c}
        !            45:        };
        !            46: 
        !            47: 
        !            48: /*--- poly_tan() ------------------------------------------------------------+
        !            49:  |                                                                           |
        !            50:  +---------------------------------------------------------------------------*/
        !            51: void   poly_tan(FPU_REG *arg, FPU_REG *y_reg)
        !            52: {
        !            53:   char         invert = 0;
        !            54:   short                exponent;
        !            55:   FPU_REG       odd_poly, even_poly, pos_poly, neg_poly;
        !            56:   FPU_REG       argSq;
        !            57:   long long     arg_signif, argSqSq;
        !            58:   
        !            59: 
        !            60:   exponent = arg->exp - EXP_BIAS;
        !            61:   
        !            62:   if ( arg->tag == TW_Zero )
        !            63:     {
        !            64:       /* Return 0.0 */
        !            65:       reg_move(&CONST_Z, y_reg);
        !            66:       return;
        !            67:     }
        !            68: 
        !            69:   if ( exponent >= -1 )
        !            70:     {
        !            71:       /* argument is in the range  [0.5 .. 1.0] */
        !            72:       if ( exponent >= 0 )
        !            73:        {
        !            74: #ifdef PARANOID
        !            75:          if ( (exponent == 0) && 
        !            76:              (arg->sigl == 0) && (arg->sigh == 0x80000000) )
        !            77: #endif PARANOID
        !            78:            {
        !            79:              arith_overflow(y_reg);
        !            80:              return;
        !            81:            }
        !            82: #ifdef PARANOID
        !            83:          EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
        !            84: #endif PARANOID
        !            85:        }
        !            86:       /* The argument is in the range  [0.5 .. 1.0) */
        !            87:       /* Convert the argument to a number in the range  (0.0 .. 0.5] */
        !            88:       *((long long *)(&arg->sigl)) = - *((long long *)(&arg->sigl));
        !            89:       normalize(arg);  /* Needed later */
        !            90:       exponent = arg->exp - EXP_BIAS;
        !            91:       invert = 1;
        !            92:     }
        !            93: 
        !            94: #ifdef PARANOID
        !            95:   if ( arg->sign != 0 )        /* Can't hack a number < 0.0 */
        !            96:     { arith_invalid(y_reg); return; }
        !            97: #endif PARANOID
        !            98: 
        !            99:   *(long long *)&arg_signif = *(long long *)&(arg->sigl);
        !           100:   if ( exponent < -1 )
        !           101:     {
        !           102:       /* shift the argument right by the required places */
        !           103:       if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
        !           104:        arg_signif++;   /* round up */
        !           105:     }
        !           106: 
        !           107:   mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
        !           108:   mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
        !           109: 
        !           110:   /* will be a valid positive nr with expon = 0 */
        !           111:   *(short *)&(pos_poly.sign) = 0;
        !           112:   pos_poly.exp = EXP_BIAS;
        !           113: 
        !           114:   /* Do the basic fixed point polynomial evaluation */
        !           115:   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq, oddplterms, HIPOWERop-1);
        !           116: 
        !           117:   /* will be a valid positive nr with expon = 0 */
        !           118:   *(short *)&(neg_poly.sign) = 0;
        !           119:   neg_poly.exp = EXP_BIAS;
        !           120: 
        !           121:   /* Do the basic fixed point polynomial evaluation */
        !           122:   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, oddnegterms, HIPOWERon-1);
        !           123:   mul64((long long *)(&argSq.sigl), (long long *)(&neg_poly.sigl),
        !           124:        (long long *)(&neg_poly.sigl));
        !           125: 
        !           126:   /* Subtract the mantissas */
        !           127:   *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
        !           128: 
        !           129:   /* Convert to 64 bit signed-compatible */
        !           130:   pos_poly.exp -= 1;
        !           131: 
        !           132:   reg_move(&pos_poly, &odd_poly);
        !           133:   normalize(&odd_poly);
        !           134:   
        !           135:   reg_mul(&odd_poly, arg, &odd_poly);
        !           136:   reg_u_add(&odd_poly, arg, &odd_poly);        /* This is just the odd polynomial */
        !           137: 
        !           138: 
        !           139:   /* will be a valid positive nr with expon = 0 */
        !           140:   *(short *)&(pos_poly.sign) = 0;
        !           141:   pos_poly.exp = EXP_BIAS;
        !           142:   
        !           143:   /* Do the basic fixed point polynomial evaluation */
        !           144:   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq, evenplterms, HIPOWERep-1);
        !           145:   mul64((long long *)(&argSq.sigl),
        !           146:        (long long *)(&pos_poly.sigl), (long long *)(&pos_poly.sigl));
        !           147:   
        !           148:   /* will be a valid positive nr with expon = 0 */
        !           149:   *(short *)&(neg_poly.sign) = 0;
        !           150:   neg_poly.exp = EXP_BIAS;
        !           151: 
        !           152:   /* Do the basic fixed point polynomial evaluation */
        !           153:   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, evennegterms, HIPOWERen-1);
        !           154: 
        !           155:   /* Subtract the mantissas */
        !           156:   *((long long *)(&neg_poly.sigl)) -= *((long long *)(&pos_poly.sigl));
        !           157:   /* and multiply by argSq */
        !           158: 
        !           159:   /* Convert argSq to a valid reg number */
        !           160:   *(short *)&(argSq.sign) = 0;
        !           161:   argSq.exp = EXP_BIAS - 1;
        !           162:   normalize(&argSq);
        !           163: 
        !           164:   /* Convert to 64 bit signed-compatible */
        !           165:   neg_poly.exp -= 1;
        !           166: 
        !           167:   reg_move(&neg_poly, &even_poly);
        !           168:   normalize(&even_poly);
        !           169: 
        !           170:   reg_mul(&even_poly, &argSq, &even_poly);
        !           171:   reg_add(&even_poly, &argSq, &even_poly);
        !           172:   reg_sub(&CONST_1, &even_poly, &even_poly);  /* This is just the even polynomial */
        !           173: 
        !           174:   /* Now ready to copy the results */
        !           175:   if ( invert )
        !           176:     { reg_div(&even_poly, &odd_poly, y_reg); }
        !           177:   else
        !           178:     { reg_div(&odd_poly, &even_poly, y_reg); }
        !           179: 
        !           180: }

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