sys/platform/pc32/gnu/fpemul/poly_tan.c

   1 /*
   2  *  poly_tan.c
   3  *
   4  * Compute the tan of a FPU_REG, using a polynomial approximation.
   5  *
   6  *
   7  * Copyright (C) 1992,1993,1994
   8  *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
   9  *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
  10  * All rights reserved.
  11  *
  12  * This copyright notice covers the redistribution and use of the
  13  * FPU emulator developed by W. Metzenthen. It covers only its use
  14  * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
  15  * use is not permitted under this copyright.
  16  *
  17  * Redistribution and use in source and binary forms, with or without
  18  * modification, are permitted provided that the following conditions
  19  * are met:
  20  * 1. Redistributions of source code must retain the above copyright
  21  *    notice, this list of conditions and the following disclaimer.
  22  * 2. Redistributions in binary form must include information specifying
  23  *    that source code for the emulator is freely available and include
  24  *    either:
  25  *      a) an offer to provide the source code for a nominal distribution
  26  *         fee, or
  27  *      b) list at least two alternative methods whereby the source
  28  *         can be obtained, e.g. a publically accessible bulletin board
  29  *         and an anonymous ftp site from which the software can be
  30  *         downloaded.
  31  * 3. All advertising materials specifically mentioning features or use of
  32  *    this emulator must acknowledge that it was developed by W. Metzenthen.
  33  * 4. The name of W. Metzenthen may not be used to endorse or promote
  34  *    products derived from this software without specific prior written
  35  *    permission.
  36  *
  37  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  38  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
  39  * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
  40  * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  41  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  42  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  43  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  44  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  45  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  46  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  47  *
  48  *
  49  * The purpose of this copyright, based upon the Berkeley copyright, is to
  50  * ensure that the covered software remains freely available to everyone.
  51  *
  52  * The software (with necessary differences) is also available, but under
  53  * the terms of the GNU copyleft, for the Linux operating system and for
  54  * the djgpp ms-dos extender.
  55  *
  56  * W. Metzenthen   June 1994.
  57  *
  58  *
  59  * $FreeBSD: src/sys/gnu/i386/fpemul/poly_tan.c,v 1.10 1999/08/28 00:42:54 peter Exp $
  60  *
  61  */
  62
  63 #include <gnu/i386/fpemul/reg_constant.h>
  64 #include <gnu/i386/fpemul/control_w.h>
  65
  66
  67 #define HIPOWERop       3       /* odd poly, positive terms */
  68 static unsigned short oddplterms[HIPOWERop][4] =
  69 {
  70         {0x846a, 0x42d1, 0xb544, 0x921f},
  71         {0x6fb2, 0x0215, 0x95c0, 0x099c},
  72         {0xfce6, 0x0cc8, 0x1c9a, 0x0000}
  73 };
  74 #define HIPOWERon       2       /* odd poly, negative terms */
  75 static unsigned short oddnegterms[HIPOWERon][4] =
  76 {
  77         {0x6906, 0xe205, 0x25c8, 0x8838},
  78         {0x1dd7, 0x3fe3, 0x944e, 0x002c}
  79 };
  80 #define HIPOWERep       2       /* even poly, positive terms */
  81 static unsigned short evenplterms[HIPOWERep][4] =
  82 {
  83         {0xdb8f, 0x3761, 0x1432, 0x2acf},
  84         {0x16eb, 0x13c1, 0x3099, 0x0003}
  85 };
  86 #define HIPOWERen       2       /* even poly, negative terms */
  87 static unsigned short evennegterms[HIPOWERen][4] =
  88 {
  89         {0x3a7c, 0xe4c5, 0x7f87, 0x2945},
  90         {0x572b, 0x664c, 0xc543, 0x018c}
  91 };
  92
  93
  94 /*--- poly_tan() ------------------------------------------------------------+
  95  |                                                                           |
  96  +---------------------------------------------------------------------------*/
  97 void
  98 poly_tan(FPU_REG * arg, FPU_REG * y_reg)
  99 {
 100         char    invert = 0;
 101         short   exponent;
 102         FPU_REG odd_poly, even_poly, pos_poly, neg_poly;
 103         FPU_REG argSq;
 104         long long arg_signif, argSqSq;
 105
 106
 107         exponent = arg->exp - EXP_BIAS;
 108
 109         if (arg->tag == TW_Zero) {
 110                 /* Return 0.0 */
 111                 reg_move(&CONST_Z, y_reg);
 112                 return;
 113         }
 114         if (exponent >= -1) {
 115                 /* argument is in the range  [0.5 .. 1.0] */
 116                 if (exponent >= 0) {
 117 #ifdef PARANOID
 118                         if ((exponent == 0) &&
 119                             (arg->sigl == 0) && (arg->sigh == 0x80000000))
 120 #endif                          /* PARANOID */
 121                         {
 122                                 arith_overflow(y_reg);
 123                                 return;
 124                         }
 125 #ifdef PARANOID
 126                         EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic
 127                                                          * error */
 128                         return;
 129 #endif                          /* PARANOID */
 130                 }
 131                 /* The argument is in the range  [0.5 .. 1.0) */
 132                 /* Convert the argument to a number in the range  (0.0 .. 0.5] */
 133                 *((long long *) (&arg->sigl)) = -*((long long *) (&arg->sigl));
 134                 normalize(arg); /* Needed later */
 135                 exponent = arg->exp - EXP_BIAS;
 136                 invert = 1;
 137         }
 138 #ifdef PARANOID
 139         if (arg->sign != 0) {   /* Can't hack a number < 0.0 */
 140                 arith_invalid(y_reg);
 141                 return;
 142         }                       /* Need a positive number */
 143 #endif                          /* PARANOID */
 144
 145         *(long long *) &arg_signif = *(long long *) &(arg->sigl);
 146         if (exponent < -1) {
 147                 /* shift the argument right by the required places */
 148                 if (shrx(&arg_signif, -1 - exponent) >= (unsigned)0x80000000)
 149                         arg_signif++;   /* round up */
 150         }
 151         mul64(&arg_signif, &arg_signif, (long long *) (&argSq.sigl));
 152         mul64((long long *) (&argSq.sigl), (long long *) (&argSq.sigl), &argSqSq);
 153
 154         /* will be a valid positive nr with expon = 0 */
 155         *(short *) &(pos_poly.sign) = 0;
 156         pos_poly.exp = EXP_BIAS;
 157
 158         /* Do the basic fixed point polynomial evaluation */
 159         polynomial((u_int *) &pos_poly.sigl, (unsigned *) &argSqSq, oddplterms, HIPOWERop - 1);
 160
 161         /* will be a valid positive nr with expon = 0 */
 162         *(short *) &(neg_poly.sign) = 0;
 163         neg_poly.exp = EXP_BIAS;
 164
 165         /* Do the basic fixed point polynomial evaluation */
 166         polynomial((u_int *) &neg_poly.sigl, (unsigned *) &argSqSq, oddnegterms, HIPOWERon - 1);
 167         mul64((long long *) (&argSq.sigl), (long long *) (&neg_poly.sigl),
 168             (long long *) (&neg_poly.sigl));
 169
 170         /* Subtract the mantissas */
 171         *((long long *) (&pos_poly.sigl)) -= *((long long *) (&neg_poly.sigl));
 172
 173         /* Convert to 64 bit signed-compatible */
 174         pos_poly.exp -= 1;
 175
 176         reg_move(&pos_poly, &odd_poly);
 177         normalize(&odd_poly);
 178
 179         reg_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
 180         reg_u_add(&odd_poly, arg, &odd_poly, FULL_PRECISION);   /* This is just the odd
 181                                                                  * polynomial */
 182
 183
 184         /* will be a valid positive nr with expon = 0 */
 185         *(short *) &(pos_poly.sign) = 0;
 186         pos_poly.exp = EXP_BIAS;
 187
 188         /* Do the basic fixed point polynomial evaluation */
 189         polynomial((u_int *) &pos_poly.sigl, (unsigned *) &argSqSq, evenplterms, HIPOWERep - 1);
 190         mul64((long long *) (&argSq.sigl),
 191             (long long *) (&pos_poly.sigl), (long long *) (&pos_poly.sigl));
 192
 193         /* will be a valid positive nr with expon = 0 */
 194         *(short *) &(neg_poly.sign) = 0;
 195         neg_poly.exp = EXP_BIAS;
 196
 197         /* Do the basic fixed point polynomial evaluation */
 198         polynomial((u_int *) &neg_poly.sigl, (unsigned *) &argSqSq, evennegterms, HIPOWERen - 1);
 199
 200         /* Subtract the mantissas */
 201         *((long long *) (&neg_poly.sigl)) -= *((long long *) (&pos_poly.sigl));
 202         /* and multiply by argSq */
 203
 204         /* Convert argSq to a valid reg number */
 205         *(short *) &(argSq.sign) = 0;
 206         argSq.exp = EXP_BIAS - 1;
 207         normalize(&argSq);
 208
 209         /* Convert to 64 bit signed-compatible */
 210         neg_poly.exp -= 1;
 211
 212         reg_move(&neg_poly, &even_poly);
 213         normalize(&even_poly);
 214
 215         reg_mul(&even_poly, &argSq, &even_poly, FULL_PRECISION);
 216         reg_add(&even_poly, &argSq, &even_poly, FULL_PRECISION);
 217         reg_sub(&CONST_1, &even_poly, &even_poly, FULL_PRECISION);      /* This is just the even
 218                                                                          * polynomial */
 219
 220         /* Now ready to copy the results */
 221         if (invert) {
 222                 reg_div(&even_poly, &odd_poly, y_reg, FULL_PRECISION);
 223         } else {
 224                 reg_div(&odd_poly, &even_poly, y_reg, FULL_PRECISION);
 225         }
 226
 227 }