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