sys/i386/gnu/fpemul/poly_atan.c

   1 /*
   2  *  p_atan.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_atan.c,v 1.10 1999/08/28 00:42:53 peter Exp $
  60  * $DragonFly: src/sys/i386/gnu/fpemul/Attic/poly_atan.c,v 1.3 2003/08/07 21:17:20 dillon Exp $
  61  *
  62  */
  63
  64 #include "reg_constant.h"
  65 #include "control_w.h"
  66
  67
  68 #define HIPOWERon       6       /* odd poly, negative terms */
  69 static unsigned oddnegterms[HIPOWERon][2] =
  70 {
  71         {0x00000000, 0x00000000},       /* for + 1.0 */
  72         {0x763b6f3d, 0x1adc4428},
  73         {0x20f0630b, 0x0502909d},
  74         {0x4e825578, 0x0198ce38},
  75         {0x22b7cb87, 0x008da6e3},
  76         {0x9b30ca03, 0x00239c79}
  77 };
  78 #define HIPOWERop       6       /* odd poly, positive terms */
  79 static unsigned oddplterms[HIPOWERop][2] =
  80 {
  81         {0xa6f67cb8, 0x94d910bd},
  82         {0xa02ffab4, 0x0a43cb45},
  83         {0x04265e6b, 0x02bf5655},
  84         {0x0a728914, 0x00f280f7},
  85         {0x6d640e01, 0x004d6556},
  86         {0xf1dd2dbf, 0x000a530a}
  87 };
  88
  89
  90 static unsigned denomterm[2] =
  91 {0xfc4bd208, 0xea2e6612};
  92
  93
  94 static void poly_add_1(FPU_REG * src);
  95
  96 /*--- poly_atan() -----------------------------------------------------------+
  97  |                                                                           |
  98  +---------------------------------------------------------------------------*/
  99 void
 100 poly_atan(FPU_REG * arg)
 101 {
 102         char    recursions = 0;
 103         short   exponent;
 104         FPU_REG odd_poly, even_poly, pos_poly, neg_poly;
 105         FPU_REG argSq;
 106         long long arg_signif, argSqSq;
 107
 108
 109 #ifdef PARANOID
 110         if (arg->sign != 0) {   /* Can't hack a number < 0.0 */
 111                 arith_invalid(arg);
 112                 return;
 113         }                       /* Need a positive number */
 114 #endif                          /* PARANOID */
 115
 116         exponent = arg->exp - EXP_BIAS;
 117
 118         if (arg->tag == TW_Zero) {
 119                 /* Return 0.0 */
 120                 reg_move(&CONST_Z, arg);
 121                 return;
 122         }
 123         if (exponent >= -2) {
 124                 /* argument is in the range  [0.25 .. 1.0] */
 125                 if (exponent >= 0) {
 126 #ifdef PARANOID
 127                         if ((exponent == 0) &&
 128                             (arg->sigl == 0) && (arg->sigh == 0x80000000))
 129 #endif                          /* PARANOID */
 130                         {
 131                                 reg_move(&CONST_PI4, arg);
 132                                 return;
 133                         }
 134 #ifdef PARANOID
 135                         EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic
 136                                                          * error */
 137 #endif                          /* PARANOID */
 138                 }
 139                 /* If the argument is greater than sqrt(2)-1 (=0.414213562...) */
 140                 /* convert the argument by an identity for atan */
 141                 if ((exponent >= -1) || (arg->sigh > 0xd413ccd0)) {
 142                         FPU_REG numerator, denom;
 143
 144                         recursions++;
 145
 146                         arg_signif = *(long long *) &(arg->sigl);
 147                         if (exponent < -1) {
 148                                 if (shrx(&arg_signif, -1 - exponent) >= (unsigned)0x80000000)
 149                                         arg_signif++;   /* round up */
 150                         }
 151                         *(long long *) &(numerator.sigl) = -arg_signif;
 152                         numerator.exp = EXP_BIAS - 1;
 153                         normalize(&numerator);  /* 1 - arg */
 154
 155                         arg_signif = *(long long *) &(arg->sigl);
 156                         if (shrx(&arg_signif, -exponent) >= (unsigned)0x80000000)
 157                                 arg_signif++;   /* round up */
 158                         *(long long *) &(denom.sigl) = arg_signif;
 159                         denom.sigh |= 0x80000000;       /* 1 + arg */
 160
 161                         arg->exp = numerator.exp;
 162                         reg_u_div(&numerator, &denom, arg, FULL_PRECISION);
 163
 164                         exponent = arg->exp - EXP_BIAS;
 165                 }
 166         }
 167         *(long long *) &arg_signif = *(long long *) &(arg->sigl);
 168
 169 #ifdef PARANOID
 170         /* This must always be true */
 171         if (exponent >= -1) {
 172                 EXCEPTION(EX_INTERNAL | 0x120); /* There must be a logic error */
 173         }
 174 #endif                          /* PARANOID */
 175
 176         /* shift the argument right by the required places */
 177         if (shrx(&arg_signif, -1 - exponent) >= (unsigned)0x80000000)
 178                 arg_signif++;   /* round up */
 179
 180         /* Now have arg_signif with binary point at the left .1xxxxxxxx */
 181         mul64(&arg_signif, &arg_signif, (long long *) (&argSq.sigl));
 182         mul64((long long *) (&argSq.sigl), (long long *) (&argSq.sigl), &argSqSq);
 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,
 190             (unsigned short (*)[4]) oddplterms, HIPOWERop - 1);
 191         mul64((long long *) (&argSq.sigl), (long long *) (&pos_poly.sigl),
 192             (long long *) (&pos_poly.sigl));
 193
 194         /* will be a valid positive nr with expon = 0 */
 195         *(short *) &(neg_poly.sign) = 0;
 196         neg_poly.exp = EXP_BIAS;
 197
 198         /* Do the basic fixed point polynomial evaluation */
 199         polynomial((u_int *) &neg_poly.sigl, (unsigned *) &argSqSq,
 200             (unsigned short (*)[4]) oddnegterms, HIPOWERon - 1);
 201
 202         /* Subtract the mantissas */
 203         *((long long *) (&pos_poly.sigl)) -= *((long long *) (&neg_poly.sigl));
 204
 205         reg_move(&pos_poly, &odd_poly);
 206         poly_add_1(&odd_poly);
 207
 208         /* The complete odd polynomial */
 209         reg_u_mul(&odd_poly, arg, &odd_poly, FULL_PRECISION);
 210
 211         /* will be a valid positive nr with expon = 0 */
 212         *(short *) &(even_poly.sign) = 0;
 213
 214         mul64((long long *) (&argSq.sigl),
 215             (long long *) (&denomterm), (long long *) (&even_poly.sigl));
 216
 217         poly_add_1(&even_poly);
 218
 219         reg_div(&odd_poly, &even_poly, arg, FULL_PRECISION);
 220
 221         if (recursions)
 222                 reg_sub(&CONST_PI4, arg, arg, FULL_PRECISION);
 223 }
 224
 225
 226 /* The argument to this function must be polynomial() compatible,
 227    i.e. have an exponent (not checked) of EXP_BIAS-1 but need not
 228    be normalized.
 229    This function adds 1.0 to the (assumed positive) argument. */
 230 static void
 231 poly_add_1(FPU_REG * src)
 232 {
 233 /* Rounding in a consistent direction produces better results
 234    for the use of this function in poly_atan. Simple truncation
 235    is used here instead of round-to-nearest. */
 236
 237 #ifdef OBSOLETE
 238         char    round = (src->sigl & 3) == 3;
 239 #endif                          /* OBSOLETE */
 240
 241         shrx(&src->sigl, 1);
 242
 243 #ifdef OBSOLETE
 244         if (round)
 245                 (*(long long *) &src->sigl)++;  /* Round to even */
 246 #endif                          /* OBSOLETE */
 247
 248         src->sigh |= 0x80000000;
 249
 250         src->exp = EXP_BIAS;
 251
 252 }