home *** CD-ROM | disk | FTP | other *** search
/ PC Extra Super CD 1998 January / PCPLUS131.iso / DJGPP / V2MISC / WMEMU2S.ZIP / LINUX / DRIVERS / FPU-EMU / POLY_TAN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-01  |  4.9 KB  |  150 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  poly_tan.c                                                               |
  3.  |                                                                           |
  4.  | Compute the tan of a FPU_REG, using a polynomial approximation.           |
  5.  |                                                                           |
  6.  | Copyright (C) 1992,1993                                                   |
  7.  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
  8.  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12.  
  13. #include "exception.h"
  14. #include "reg_constant.h"
  15. #include "fpu_emu.h"
  16. #include "control_w.h"
  17.  
  18.  
  19. #define    HIPOWERop    3    /* odd poly, positive terms */
  20. static unsigned short const    oddplterms[HIPOWERop][4] =
  21.     {
  22.     { 0x846a, 0x42d1, 0xb544, 0x921f},
  23.     { 0x6fb2, 0x0215, 0x95c0, 0x099c},
  24.     { 0xfce6, 0x0cc8, 0x1c9a, 0x0000}
  25.     };
  26.  
  27. #define    HIPOWERon    2    /* odd poly, negative terms */
  28. static unsigned short const    oddnegterms[HIPOWERon][4] =
  29.     {
  30.     { 0x6906, 0xe205, 0x25c8, 0x8838},
  31.     { 0x1dd7, 0x3fe3, 0x944e, 0x002c}
  32.     };
  33.  
  34. #define    HIPOWERep    2    /* even poly, positive terms */
  35. static unsigned short const    evenplterms[HIPOWERep][4] =
  36.     {
  37.     { 0xdb8f, 0x3761, 0x1432, 0x2acf},
  38.     { 0x16eb, 0x13c1, 0x3099, 0x0003}
  39.     };
  40.  
  41. #define    HIPOWERen    2    /* even poly, negative terms */
  42. static unsigned short const    evennegterms[HIPOWERen][4] =
  43.     {
  44.     { 0x3a7c, 0xe4c5, 0x7f87, 0x2945},
  45.     { 0x572b, 0x664c, 0xc543, 0x018c}
  46.     };
  47.  
  48.  
  49. /*--- poly_tan() ------------------------------------------------------------+
  50.  |                                                                           |
  51.  +---------------------------------------------------------------------------*/
  52. void    poly_tan(FPU_REG const *arg, FPU_REG *result, int invert)
  53. {
  54.   short        exponent;
  55.   FPU_REG       odd_poly, even_poly, pos_poly, neg_poly;
  56.   FPU_REG       argSq;
  57.   unsigned long long     arg_signif, argSqSq;
  58.   
  59.  
  60.   exponent = arg->exp - EXP_BIAS;
  61.  
  62. #ifdef PARANOID
  63.   if ( arg->sign != 0 )    /* Can't hack a number < 0.0 */
  64.     { arith_invalid(result); return; }  /* Need a positive number */
  65. #endif PARANOID
  66.  
  67.   arg_signif = significand(arg);
  68.   if ( exponent < -1 )
  69.     {
  70.       /* shift the argument right by the required places */
  71.       if ( shrx(&arg_signif, -1-exponent) >= 0x80000000U )
  72.     arg_signif++;    /* round up */
  73.     }
  74.  
  75.   mul64(&arg_signif, &arg_signif, &significand(&argSq));
  76.   mul64(&significand(&argSq), &significand(&argSq), &argSqSq);
  77.  
  78.   /* will be a valid positive nr with expon = 0 */
  79.   *(short *)&(pos_poly.sign) = 0;
  80.   pos_poly.exp = EXP_BIAS;
  81.  
  82.   /* Do the basic fixed point polynomial evaluation */
  83.   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq, oddplterms, HIPOWERop-1);
  84.  
  85.   /* will be a valid positive nr with expon = 0 */
  86.   *(short *)&(neg_poly.sign) = 0;
  87.   neg_poly.exp = EXP_BIAS;
  88.  
  89.   /* Do the basic fixed point polynomial evaluation */
  90.   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, oddnegterms, HIPOWERon-1);
  91.   mul64(&significand(&argSq), &significand(&neg_poly),
  92.     &significand(&neg_poly));
  93.  
  94.   /* Subtract the mantissas */
  95.   significand(&pos_poly) -= significand(&neg_poly);
  96.  
  97.   /* Convert to 64 bit signed-compatible */
  98.   pos_poly.exp -= 1;
  99.  
  100.   reg_move(&pos_poly, &odd_poly);
  101.   normalize(&odd_poly);
  102.   
  103.   reg_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
  104.   /* Complete the odd polynomial. */
  105.   reg_u_add(&odd_poly, arg, &odd_poly, FULL_PRECISION);
  106.  
  107.   /* will be a valid positive nr with expon = 0 */
  108.   *(short *)&(pos_poly.sign) = 0;
  109.   pos_poly.exp = EXP_BIAS;
  110.   
  111.   /* Do the basic fixed point polynomial evaluation */
  112.   polynomial(&pos_poly.sigl, (unsigned *)&argSqSq, evenplterms, HIPOWERep-1);
  113.   mul64(&significand(&argSq),
  114.     &significand(&pos_poly), &significand(&pos_poly));
  115.   
  116.   /* will be a valid positive nr with expon = 0 */
  117.   *(short *)&(neg_poly.sign) = 0;
  118.   neg_poly.exp = EXP_BIAS;
  119.  
  120.   /* Do the basic fixed point polynomial evaluation */
  121.   polynomial(&neg_poly.sigl, (unsigned *)&argSqSq, evennegterms, HIPOWERen-1);
  122.  
  123.   /* Subtract the mantissas */
  124.   significand(&neg_poly) -= significand(&pos_poly);
  125.   /* and multiply by argSq */
  126.  
  127.   /* Convert argSq to a valid reg number */
  128.   *(short *)&(argSq.sign) = 0;
  129.   argSq.exp = EXP_BIAS - 1;
  130.   normalize(&argSq);
  131.  
  132.   /* Convert to 64 bit signed-compatible */
  133.   neg_poly.exp -= 1;
  134.  
  135.   reg_move(&neg_poly, &even_poly);
  136.   normalize(&even_poly);
  137.  
  138.   reg_mul(&even_poly, &argSq, &even_poly, FULL_PRECISION);
  139.   reg_add(&even_poly, &argSq, &even_poly, FULL_PRECISION);
  140.   /* Complete the even polynomial */
  141.   reg_sub(&CONST_1, &even_poly, &even_poly, FULL_PRECISION);
  142.  
  143.   /* Now ready to copy the results */
  144.   if ( invert )
  145.     { reg_div(&even_poly, &odd_poly, result, FULL_PRECISION); }
  146.   else
  147.     { reg_div(&odd_poly, &even_poly, result, FULL_PRECISION); }
  148.  
  149. }
  150.