Return to poly_atan.c CVS log

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

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

1.1     ! root        1: /*---------------------------------------------------------------------------+
        !             2:  |  p_atan.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        HIPOWERon       6       /* odd poly, negative terms */
        !            19: static unsigned oddnegterms[HIPOWERon][2] =
        !            20: {
        !            21:   { 0x00000000, 0x00000000 }, /* for + 1.0 */
        !            22:   { 0x763b6f3d, 0x1adc4428 },
        !            23:   { 0x20f0630b, 0x0502909d },
        !            24:   { 0x4e825578, 0x0198ce38 },
        !            25:   { 0x22b7cb87, 0x008da6e3 },
        !            26:   { 0x9b30ca03, 0x00239c79 }
        !            27: } ;
        !            28: 
        !            29: #define        HIPOWERop       6       /* odd poly, positive terms */
        !            30: static unsigned        oddplterms[HIPOWERop][2] =
        !            31: {
        !            32:   { 0xa6f67cb8, 0x94d910bd },
        !            33:   { 0xa02ffab4, 0x0a43cb45 },
        !            34:   { 0x04265e6b, 0x02bf5655 },
        !            35:   { 0x0a728914, 0x00f280f7 },
        !            36:   { 0x6d640e01, 0x004d6556 },
        !            37:   { 0xf1dd2dbf, 0x000a530a }
        !            38: };
        !            39: 
        !            40: 
        !            41: static unsigned denomterm[2] =
        !            42: { 0xfc4bd208, 0xea2e6612 };
        !            43: 
        !            44: 
        !            45: 
        !            46: /*--- poly_atan() -----------------------------------------------------------+
        !            47:  |                                                                           |
        !            48:  +---------------------------------------------------------------------------*/
        !            49: void   poly_atan(FPU_REG *arg)
        !            50: {
        !            51:   char         recursions = 0;
        !            52:   short                exponent;
        !            53:   FPU_REG       odd_poly, even_poly, pos_poly, neg_poly;
        !            54:   FPU_REG       argSq;
        !            55:   long long     arg_signif, argSqSq;
        !            56:   
        !            57: 
        !            58: #ifdef PARANOID
        !            59:   if ( arg->sign != 0 )        /* Can't hack a number < 0.0 */
        !            60:     { arith_invalid(arg); return; }
        !            61: #endif PARANOID
        !            62: 
        !            63:   exponent = arg->exp - EXP_BIAS;
        !            64:   
        !            65:   if ( arg->tag == TW_Zero )
        !            66:     {
        !            67:       /* Return 0.0 */
        !            68:       reg_move(&CONST_Z, arg);
        !            69:       return;
        !            70:     }
        !            71:   
        !            72:   if ( exponent >= -2 )
        !            73:     {
        !            74:       /* argument is in the range  [0.25 .. 1.0] */
        !            75:       if ( exponent >= 0 )
        !            76:        {
        !            77: #ifdef PARANOID
        !            78:          if ( (exponent == 0) && 
        !            79:              (arg->sigl == 0) && (arg->sigh == 0x80000000) )
        !            80: #endif PARANOID
        !            81:            {
        !            82:              reg_move(&CONST_PI4, arg);
        !            83:              return;
        !            84:            }
        !            85: #ifdef PARANOID
        !            86:          EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
        !            87: #endif PARANOID
        !            88:        }
        !            89: 
        !            90:       /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
        !            91:       /* convert the argument by an identity for atan */
        !            92:       if ( (exponent >= -1) || (arg->sigh > 0xd413ccd0) )
        !            93:        {
        !            94:          FPU_REG numerator, denom;
        !            95: 
        !            96:          recursions++;
        !            97: 
        !            98:          arg_signif = *(long long *)&(arg->sigl);
        !            99:          if ( exponent < -1 )
        !           100:            {
        !           101:              if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
        !           102:                arg_signif++;   /* round up */
        !           103:            }
        !           104:          *(long long *)&(numerator.sigl) = -arg_signif;
        !           105:          numerator.exp = EXP_BIAS - 1;
        !           106:          normalize(&numerator);                       /* 1 - arg */
        !           107: 
        !           108:          arg_signif = *(long long *)&(arg->sigl);
        !           109:          if ( shrx(&arg_signif, -exponent) >= 0x80000000U )
        !           110:            arg_signif++;       /* round up */
        !           111:          *(long long *)&(denom.sigl) = arg_signif;
        !           112:          denom.sigh |= 0x80000000;                    /* 1 + arg */
        !           113: 
        !           114:          arg->exp = numerator.exp;
        !           115:          reg_u_div((long long *)&(numerator.sigl),
        !           116:                    (long long *)&(denom.sigl), arg);
        !           117: 
        !           118:          exponent = arg->exp - EXP_BIAS;
        !           119:        }
        !           120:     }
        !           121: 
        !           122:   *(long long *)&arg_signif = *(long long *)&(arg->sigl);
        !           123: 
        !           124: #ifdef PARANOID
        !           125:   /* This must always be true */
        !           126:   if ( exponent >= -1 )
        !           127:     {
        !           128:       EXCEPTION(EX_INTERNAL|0x120);    /* There must be a logic error */
        !           129:     }
        !           130: #endif PARANOID
        !           131: 
        !           132:   /* shift the argument right by the required places */
        !           133:   if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
        !           134:     arg_signif++;      /* round up */
        !           135:   
        !           136:   /* Now have arg_signif with binary point at the left
        !           137:      .1xxxxxxxx */
        !           138:   mul64(&arg_signif, &arg_signif, (long long *)(&argSq.sigl));
        !           139:   mul64((long long *)(&argSq.sigl), (long long *)(&argSq.sigl), &argSqSq);
        !           140: 
        !           141:   /* will be a valid positive nr with expon = 0 */
        !           142:   *(short *)&(pos_poly.sign) = 0;
        !           143:   pos_poly.exp = EXP_BIAS;
        !           144: 
        !           145:   /* Do the basic fixed point polynomial evaluation */
        !           146:   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq,
        !           147:             (unsigned short (*)[4])oddplterms, HIPOWERop-1);
        !           148:   mul64((long long *)(&argSq.sigl), (long long *)(&pos_poly.sigl),
        !           149:        (long long *)(&pos_poly.sigl));
        !           150: 
        !           151:   /* will be a valid positive nr with expon = 0 */
        !           152:   *(short *)&(neg_poly.sign) = 0;
        !           153:   neg_poly.exp = EXP_BIAS;
        !           154: 
        !           155:   /* Do the basic fixed point polynomial evaluation */
        !           156:   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq,
        !           157:             (unsigned short (*)[4])oddnegterms, HIPOWERon-1);
        !           158: 
        !           159:   /* Subtract the mantissas */
        !           160:   *((long long *)(&pos_poly.sigl)) -= *((long long *)(&neg_poly.sigl));
        !           161: 
        !           162:   reg_move(&pos_poly, &odd_poly);
        !           163:   poly_add_1(&odd_poly);
        !           164: 
        !           165:   reg_u_mul(&odd_poly, arg, &odd_poly);        /* The complete odd polynomial */
        !           166:   odd_poly.exp -= EXP_BIAS - 1;
        !           167: 
        !           168:   /* will be a valid positive nr with expon = 0 */
        !           169:   *(short *)&(even_poly.sign) = 0;
        !           170: 
        !           171:   mul64((long long *)(&argSq.sigl),
        !           172:        (long long *)(&denomterm), (long long *)(&even_poly.sigl));
        !           173: 
        !           174:   poly_add_1(&even_poly);
        !           175: 
        !           176:   reg_div(&odd_poly, &even_poly, arg);
        !           177: 
        !           178:   if ( recursions )
        !           179:     reg_sub(&CONST_PI4, arg, arg);
        !           180: }
        !           181: 
        !           182: 
        !           183: /* The argument to this function must be polynomial() compatible,
        !           184:    i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
        !           185:    be normalized.
        !           186:    This function adds 1.0 to the (assumed positive) argument. */
        !           187: void poly_add_1(FPU_REG *src)
        !           188: {
        !           189: /* Rounding in a consistent direction produces better results
        !           190:    for the use of this function in poly_atan. Simple truncation
        !           191:    is used here instead of round-to-nearest. */
        !           192: 
        !           193: #ifdef OBSOLETE
        !           194: char round = (src->sigl & 3) == 3;
        !           195: #endif OBSOLETE
        !           196: 
        !           197: shrx(&src->sigl, 1);
        !           198: 
        !           199: #ifdef OBSOLETE
        !           200: if ( round ) (*(long long *)&src->sigl)++;   /* Round to even */
        !           201: #endif OBSOLETE
        !           202: 
        !           203: src->sigh |= 0x80000000;
        !           204: 
        !           205: src->exp = EXP_BIAS;
        !           206: 
        !           207: }

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