root/drivers/FPU-emu/poly_tan.c

/* */

DEFINITIONS

1. poly_tan

   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 }

/* */