Return to poly_sin.c CVS log

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

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

1.1     ! root        1: /*---------------------------------------------------------------------------+
        !             2:  |  poly_sin.c                                                               |
        !             3:  |                                                                           |
        !             4:  |  Computation of an approximation of the sin function by a polynomial      |
        !             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        HIPOWER 5
        !            19: static unsigned short  lterms[HIPOWER][4] =
        !            20:        {
        !            21:        { 0x846a, 0x42d1, 0xb544, 0x921f},
        !            22:        { 0xe110, 0x75aa, 0xbc67, 0x1466},
        !            23:        { 0x503d, 0xa43f, 0x83c1, 0x000a},
        !            24:        { 0x8f9d, 0x7a19, 0x00f4, 0x0000},
        !            25:        { 0xda03, 0x06aa, 0x0000, 0x0000},
        !            26:        };
        !            27: 
        !            28: static unsigned short  negterms[HIPOWER][4] =
        !            29:        {
        !            30:        { 0x95ed, 0x2df2, 0xe731, 0xa55d},
        !            31:        { 0xd159, 0xe62b, 0xd2cc, 0x0132},
        !            32:        { 0x6342, 0xe9fb, 0x3c60, 0x0000},
        !            33:        { 0x6256, 0xdf5a, 0x0002, 0x0000},
        !            34:        { 0xf279, 0x000b, 0x0000, 0x0000},
        !            35:        };
        !            36: 
        !            37: 
        !            38: /*--- poly_sine() -----------------------------------------------------------+
        !            39:  |                                                                           |
        !            40:  +---------------------------------------------------------------------------*/
        !            41: void   poly_sine(FPU_REG *arg, FPU_REG *result)
        !            42: {
        !            43:   short        exponent;
        !            44:   FPU_REG      Xx, Xx2, Xx4, accum, negaccum;
        !            45:   
        !            46:   
        !            47:   exponent = arg->exp - EXP_BIAS;
        !            48:   
        !            49:   if ( arg->tag == TW_Zero )
        !            50:     {
        !            51:       /* Return 0.0 */
        !            52:       reg_move(&CONST_Z, result);
        !            53:       return;
        !            54:     }
        !            55:   
        !            56: #ifdef PARANOID
        !            57:   if ( arg->sign != 0 )        /* Can't hack a number < 0.0 */
        !            58:     {
        !            59:       EXCEPTION(EX_Invalid);
        !            60:       reg_move(&CONST_QNaN, result);
        !            61:       return;
        !            62:     }
        !            63:   
        !            64:   if ( exponent >= 0 ) /* Can't hack a number > 1.0 */
        !            65:     {
        !            66:       if ( (exponent == 0) && (arg->sigl == 0) && (arg->sigh == 0x80000000) )
        !            67:        {
        !            68:          reg_move(&CONST_1, result);
        !            69:          return;
        !            70:        }
        !            71:       EXCEPTION(EX_Invalid);
        !            72:       reg_move(&CONST_QNaN, result);
        !            73:       return;
        !            74:     }
        !            75: #endif PARANOID
        !            76:   
        !            77:   Xx.sigl = arg->sigl;
        !            78:   Xx.sigh = arg->sigh;
        !            79:   if ( exponent < -1 )
        !            80:     {
        !            81:       /* shift the argument right by the required places */
        !            82:       if ( shrx(&(Xx.sigl), -1-exponent) >= 0x80000000U )
        !            83:        (*((long long *)(&(Xx.sigl))))++;       /* round up */
        !            84:     }
        !            85:   
        !            86:   mul64((long long *)&(Xx.sigl), (long long *)&(Xx.sigl),
        !            87:        (long long *)&(Xx2.sigl));
        !            88:   mul64((long long *)&(Xx2.sigl), (long long *)&(Xx2.sigl),
        !            89:        (long long *)&(Xx4.sigl));
        !            90:   
        !            91:   /* will be a valid positive nr with expon = 0 */
        !            92:   *(short *)&(accum.sign) = 0;
        !            93:   accum.exp = 0;
        !            94: 
        !            95:   /* Do the basic fixed point polynomial evaluation */
        !            96:   polynomial(&(accum.sigl), &(Xx4.sigl), lterms, HIPOWER-1);
        !            97:   
        !            98:   /* will be a valid positive nr with expon = 0 */
        !            99:   *(short *)&(negaccum.sign) = 0;
        !           100:   negaccum.exp = 0;
        !           101:   
        !           102:   /* Do the basic fixed point polynomial evaluation */
        !           103:   polynomial(&(negaccum.sigl), &(Xx4.sigl), negterms, HIPOWER-1);
        !           104:   mul64((long long *)&(Xx2.sigl), (long long *)&(negaccum.sigl),
        !           105:        (long long *)&(negaccum.sigl));
        !           106: 
        !           107:   /* Subtract the mantissas */
        !           108:   *((long long *)(&(accum.sigl))) -= *((long long *)(&(negaccum.sigl)));
        !           109:   
        !           110:   /* Convert to 64 bit signed-compatible */
        !           111:   accum.exp = EXP_BIAS - 1 + accum.exp;
        !           112:   
        !           113:   *(short *)&(result->sign) = *(short *)&(accum.sign);
        !           114:   result->exp = accum.exp;
        !           115:   result->sigl = accum.sigl;
        !           116:   result->sigh = accum.sigh;
        !           117:   
        !           118:   normalize(result);
        !           119: 
        !           120:   reg_mul(result, arg, result);
        !           121:   reg_u_add(result, arg, result);
        !           122:   
        !           123:   /* A small overflow may be possible... but an illegal result. */
        !           124:   if ( result->exp >= EXP_BIAS )
        !           125:     {
        !           126:       if (    (result->exp > EXP_BIAS) /* Larger or equal 2.0 */
        !           127:          || (result->sigl > 1)   /* Larger than 1.0+msb */
        !           128:          ||    (result->sigh != 0x80000000) /* Much > 1.0 */
        !           129:          )
        !           130:        {
        !           131: #ifdef DEBUGGING
        !           132:          RE_ENTRANT_CHECK_OFF
        !           133:          printf("\nEXP=%d, MS=%08x, LS=%08x\n", result->exp,
        !           134:                 result->sigh, result->sigl);
        !           135:          RE_ENTRANT_CHECK_ON
        !           136: #endif DEBUGGING
        !           137:          EXCEPTION(EX_INTERNAL|0x103);
        !           138:        }
        !           139:       
        !           140: #ifdef DEBUGGING
        !           141:       RE_ENTRANT_CHECK_OFF
        !           142:       printf("\n***CORRECTING ILLEGAL RESULT*** in poly_sin() computation\n");
        !           143:       printf("EXP=%d, MS=%08x, LS=%08x\n", result->exp,
        !           144:             result->sigh, result->sigl);
        !           145:       RE_ENTRANT_CHECK_ON
        !           146: #endif DEBUGGING
        !           147: 
        !           148:       result->sigl = 0;        /* Truncate the result to 1.00 */
        !           149:     }
        !           150: }

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