home *** CD-ROM | disk | FTP | other *** search
/ PC Extra Super CD 1998 January / PCPLUS131.iso / DJGPP / V2MISC / WMEMU2S.ZIP / LINUX / DRIVERS / FPU-EMU / POLY_SIN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1994-01-25  |  4.5 KB  |  151 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  poly_sin.c                                                               |
  3.  |                                                                           |
  4.  |  Computation of an approximation of the sin function by a polynomial      |
  5.  |                                                                           |
  6.  | Copyright (C) 1992,1993,1994                                              |
  7.  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
  8.  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12.  
  13.  
  14. #include "exception.h"
  15. #include "reg_constant.h"
  16. #include "fpu_emu.h"
  17. #include "control_w.h"
  18.  
  19.  
  20. #define    HIPOWER    5
  21. static unsigned short const    lterms[HIPOWER][4] =
  22.     {
  23.     { 0x846a, 0x42d1, 0xb544, 0x921f},
  24.     { 0xe110, 0x75aa, 0xbc67, 0x1466},
  25.     { 0x503d, 0xa43f, 0x83c1, 0x000a},
  26.     { 0x8f9d, 0x7a19, 0x00f4, 0x0000},
  27.     { 0xda03, 0x06aa, 0x0000, 0x0000},
  28.     };
  29.  
  30. static unsigned short const    negterms[HIPOWER][4] =
  31.     {
  32.     { 0x95ed, 0x2df2, 0xe731, 0xa55d},
  33.     { 0xd159, 0xe62b, 0xd2cc, 0x0132},
  34.     { 0x6342, 0xe9fb, 0x3c60, 0x0000},
  35.     { 0x6256, 0xdf5a, 0x0002, 0x0000},
  36.     { 0xf279, 0x000b, 0x0000, 0x0000},
  37.     };
  38.  
  39.  
  40. /*--- poly_sine() -----------------------------------------------------------+
  41.  |                                                                           |
  42.  +---------------------------------------------------------------------------*/
  43. void    poly_sine(FPU_REG const *arg, FPU_REG *result)
  44. {
  45.   short    exponent;
  46.   FPU_REG    fixed_arg, arg_sqrd, arg_to_4, accum, negaccum;
  47.   
  48.   
  49.   exponent = arg->exp - EXP_BIAS;
  50.   
  51.   if ( arg->tag == TW_Zero )
  52.     {
  53.       /* Return 0.0 */
  54.       reg_move(&CONST_Z, result);
  55.       return;
  56.     }
  57.   
  58. #ifdef PARANOID
  59.   if ( arg->sign != 0 )    /* Can't hack a number < 0.0 */
  60.     {
  61.       EXCEPTION(EX_Invalid);
  62.       reg_move(&CONST_QNaN, result);
  63.       return;
  64.     }
  65.   
  66.   if ( exponent >= 0 )    /* Can't hack a number > 1.0 */
  67.     {
  68.       if ( (exponent == 0) && (arg->sigl == 0) && (arg->sigh == 0x80000000) )
  69.     {
  70.       reg_move(&CONST_1, result);
  71.       return;
  72.     }
  73.       EXCEPTION(EX_Invalid);
  74.       reg_move(&CONST_QNaN, result);
  75.       return;
  76.     }
  77. #endif PARANOID
  78.   
  79.   fixed_arg.sigl = arg->sigl;
  80.   fixed_arg.sigh = arg->sigh;
  81.   if ( exponent < -1 )
  82.     {
  83.       /* shift the argument right by the required places */
  84.       if ( shrx(&(fixed_arg.sigl), -1-exponent) >= 0x80000000U )
  85.     significand(&fixed_arg)++;    /* round up */
  86.     }
  87.   
  88.   mul64(&significand(&fixed_arg), &significand(&fixed_arg),
  89.     &significand(&arg_sqrd));
  90.   mul64(&significand(&arg_sqrd), &significand(&arg_sqrd),
  91.     &significand(&arg_to_4));
  92.   
  93.   /* will be a valid positive nr with expon = 0 */
  94.   *(short *)&(accum.sign) = 0;
  95.   accum.exp = 0;
  96.  
  97.   /* Do the basic fixed point polynomial evaluation */
  98.   polynomial(&(accum.sigl), &(arg_to_4.sigl), lterms, HIPOWER-1);
  99.   
  100.   /* will be a valid positive nr with expon = 0 */
  101.   *(short *)&(negaccum.sign) = 0;
  102.   negaccum.exp = 0;
  103.   
  104.   /* Do the basic fixed point polynomial evaluation */
  105.   polynomial(&(negaccum.sigl), &(arg_to_4.sigl), negterms, HIPOWER-1);
  106.   mul64(&significand(&arg_sqrd), &significand(&negaccum),
  107.     &significand(&negaccum));
  108.  
  109.   /* Subtract the mantissas */
  110.   significand(&accum) -= significand(&negaccum);
  111.   
  112.   /* Convert to 64 bit signed-compatible */
  113.   accum.exp = EXP_BIAS - 1 + accum.exp;
  114.  
  115.   reg_move(&accum, result);
  116.  
  117.   normalize(result);
  118.  
  119.   reg_mul(result, arg, result, FULL_PRECISION);
  120.   reg_u_add(result, arg, result, FULL_PRECISION);
  121.   
  122.   if ( result->exp >= EXP_BIAS )
  123.     {
  124.       /* A small overflow may be possible... but an illegal result. */
  125.       if (    (result->exp > EXP_BIAS) /* Larger or equal 2.0 */
  126.       || (result->sigl > 1)      /* Larger than 1.0+msb */
  127.       ||    (result->sigh != 0x80000000) /* Much > 1.0 */
  128.       )
  129.     {
  130. #ifdef DEBUGGING
  131.       RE_ENTRANT_CHECK_OFF;
  132.       printk("\nEXP=%d, MS=%08x, LS=%08x\n", result->exp,
  133.          result->sigh, result->sigl);
  134.       RE_ENTRANT_CHECK_ON;
  135. #endif DEBUGGING
  136.       EXCEPTION(EX_INTERNAL|0x103);
  137.     }
  138.       
  139. #ifdef DEBUGGING
  140.       RE_ENTRANT_CHECK_OFF;
  141.       printk("\n***CORRECTING ILLEGAL RESULT*** in poly_sin() computation\n");
  142.       printk("EXP=%d, MS=%08x, LS=%08x\n", result->exp,
  143.          result->sigh, result->sigl);
  144.       RE_ENTRANT_CHECK_ON;
  145. #endif DEBUGGING
  146.  
  147.       result->sigl = 0;    /* Truncate the result to 1.00 */
  148.     }
  149.  
  150. }
  151.