arch/x86/math-emu/poly_l2.c

   1 // SPDX-License-Identifier: GPL-2.0
   2 /*---------------------------------------------------------------------------+
   3  |  poly_l2.c                                                                |
   4  |                                                                           |
   5  | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
   6  |                                                                           |
   7  | Copyright (C) 1992,1993,1994,1997                                         |
   8  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
   9  |                  E-mail   billm@suburbia.net                              |
  10  |                                                                           |
  11  |                                                                           |
  12  +---------------------------------------------------------------------------*/
  13
  14 #include "exception.h"
  15 #include "reg_constant.h"
  16 #include "fpu_emu.h"
  17 #include "fpu_system.h"
  18 #include "control_w.h"
  19 #include "poly.h"
  20
  21 static void log2_kernel(FPU_REG const *arg, u_char argsign,
  22                         Xsig * accum_result, long int *expon);
  23
  24 /*--- poly_l2() -------------------------------------------------------------+
  25  |   Base 2 logarithm by a polynomial approximation.                         |
  26  +---------------------------------------------------------------------------*/
  27 void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
  28 {
  29         long int exponent, expon, expon_expon;
  30         Xsig accumulator, expon_accum, yaccum;
  31         u_char sign, argsign;
  32         FPU_REG x;
  33         int tag;
  34
  35         exponent = exponent16(st0_ptr);
  36
  37         /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
  38         if (st0_ptr->sigh > (unsigned)0xb504f334) {
  39                 /* Treat as  sqrt(2)/2 < st0_ptr < 1 */
  40                 significand(&x) = -significand(st0_ptr);
  41                 setexponent16(&x, -1);
  42                 exponent++;
  43                 argsign = SIGN_NEG;
  44         } else {
  45                 /* Treat as  1 <= st0_ptr < sqrt(2) */
  46                 x.sigh = st0_ptr->sigh - 0x80000000;
  47                 x.sigl = st0_ptr->sigl;
  48                 setexponent16(&x, 0);
  49                 argsign = SIGN_POS;
  50         }
  51         tag = FPU_normalize_nuo(&x);
  52
  53         if (tag == TAG_Zero) {
  54                 expon = 0;
  55                 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
  56         } else {
  57                 log2_kernel(&x, argsign, &accumulator, &expon);
  58         }
  59
  60         if (exponent < 0) {
  61                 sign = SIGN_NEG;
  62                 exponent = -exponent;
  63         } else
  64                 sign = SIGN_POS;
  65         expon_accum.msw = exponent;
  66         expon_accum.midw = expon_accum.lsw = 0;
  67         if (exponent) {
  68                 expon_expon = 31 + norm_Xsig(&expon_accum);
  69                 shr_Xsig(&accumulator, expon_expon - expon);
  70
  71                 if (sign ^ argsign)
  72                         negate_Xsig(&accumulator);
  73                 add_Xsig_Xsig(&accumulator, &expon_accum);
  74         } else {
  75                 expon_expon = expon;
  76                 sign = argsign;
  77         }
  78
  79         yaccum.lsw = 0;
  80         XSIG_LL(yaccum) = significand(st1_ptr);
  81         mul_Xsig_Xsig(&accumulator, &yaccum);
  82
  83         expon_expon += round_Xsig(&accumulator);
  84
  85         if (accumulator.msw == 0) {
  86                 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
  87                 return;
  88         }
  89
  90         significand(st1_ptr) = XSIG_LL(accumulator);
  91         setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
  92
  93         tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
  94         FPU_settagi(1, tag);
  95
  96         set_precision_flag_up();        /* 80486 appears to always do this */
  97
  98         return;
  99
 100 }
 101
 102 /*--- poly_l2p1() -----------------------------------------------------------+
 103  |   Base 2 logarithm by a polynomial approximation.                         |
 104  |   log2(x+1)                                                               |
 105  +---------------------------------------------------------------------------*/
 106 int poly_l2p1(u_char sign0, u_char sign1,
 107               FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest)
 108 {
 109         u_char tag;
 110         long int exponent;
 111         Xsig accumulator, yaccum;
 112
 113         if (exponent16(st0_ptr) < 0) {
 114                 log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
 115
 116                 yaccum.lsw = 0;
 117                 XSIG_LL(yaccum) = significand(st1_ptr);
 118                 mul_Xsig_Xsig(&accumulator, &yaccum);
 119
 120                 exponent += round_Xsig(&accumulator);
 121
 122                 exponent += exponent16(st1_ptr) + 1;
 123                 if (exponent < EXP_WAY_UNDER)
 124                         exponent = EXP_WAY_UNDER;
 125
 126                 significand(dest) = XSIG_LL(accumulator);
 127                 setexponent16(dest, exponent);
 128
 129                 tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
 130                 FPU_settagi(1, tag);
 131
 132                 if (tag == TAG_Valid)
 133                         set_precision_flag_up();        /* 80486 appears to always do this */
 134         } else {
 135                 /* The magnitude of st0_ptr is far too large. */
 136
 137                 if (sign0 != SIGN_POS) {
 138                         /* Trying to get the log of a negative number. */
 139 #ifdef PECULIAR_486             /* Stupid 80486 doesn't worry about log(negative). */
 140                         changesign(st1_ptr);
 141 #else
 142                         if (arith_invalid(1) < 0)
 143                                 return 1;
 144 #endif /* PECULIAR_486 */
 145                 }
 146
 147                 /* 80486 appears to do this */
 148                 if (sign0 == SIGN_NEG)
 149                         set_precision_flag_down();
 150                 else
 151                         set_precision_flag_up();
 152         }
 153
 154         if (exponent(dest) <= EXP_UNDER)
 155                 EXCEPTION(EX_Underflow);
 156
 157         return 0;
 158
 159 }
 160
 161 #undef HIPOWER
 162 #define HIPOWER 10
 163 static const unsigned long long logterms[HIPOWER] = {
 164         0x2a8eca5705fc2ef0LL,
 165         0xf6384ee1d01febceLL,
 166         0x093bb62877cdf642LL,
 167         0x006985d8a9ec439bLL,
 168         0x0005212c4f55a9c8LL,
 169         0x00004326a16927f0LL,
 170         0x0000038d1d80a0e7LL,
 171         0x0000003141cc80c6LL,
 172         0x00000002b1668c9fLL,
 173         0x000000002c7a46aaLL
 174 };
 175
 176 static const unsigned long leadterm = 0xb8000000;
 177
 178 /*--- log2_kernel() ---------------------------------------------------------+
 179  |   Base 2 logarithm by a polynomial approximation.                         |
 180  |   log2(x+1)                                                               |
 181  +---------------------------------------------------------------------------*/
 182 static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
 183                         long int *expon)
 184 {
 185         long int exponent, adj;
 186         unsigned long long Xsq;
 187         Xsig accumulator, Numer, Denom, argSignif, arg_signif;
 188
 189         exponent = exponent16(arg);
 190         Numer.lsw = Denom.lsw = 0;
 191         XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
 192         if (argsign == SIGN_POS) {
 193                 shr_Xsig(&Denom, 2 - (1 + exponent));
 194                 Denom.msw |= 0x80000000;
 195                 div_Xsig(&Numer, &Denom, &argSignif);
 196         } else {
 197                 shr_Xsig(&Denom, 1 - (1 + exponent));
 198                 negate_Xsig(&Denom);
 199                 if (Denom.msw & 0x80000000) {
 200                         div_Xsig(&Numer, &Denom, &argSignif);
 201                         exponent++;
 202                 } else {
 203                         /* Denom must be 1.0 */
 204                         argSignif.lsw = Numer.lsw;
 205                         argSignif.midw = Numer.midw;
 206                         argSignif.msw = Numer.msw;
 207                 }
 208         }
 209
 210 #ifndef PECULIAR_486
 211         /* Should check here that  |local_arg|  is within the valid range */
 212         if (exponent >= -2) {
 213                 if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) {
 214                         /* The argument is too large */
 215                 }
 216         }
 217 #endif /* PECULIAR_486 */
 218
 219         arg_signif.lsw = argSignif.lsw;
 220         XSIG_LL(arg_signif) = XSIG_LL(argSignif);
 221         adj = norm_Xsig(&argSignif);
 222         accumulator.lsw = argSignif.lsw;
 223         XSIG_LL(accumulator) = XSIG_LL(argSignif);
 224         mul_Xsig_Xsig(&accumulator, &accumulator);
 225         shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj)));
 226         Xsq = XSIG_LL(accumulator);
 227         if (accumulator.lsw & 0x80000000)
 228                 Xsq++;
 229
 230         accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 231         /* Do the basic fixed point polynomial evaluation */
 232         polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1);
 233
 234         mul_Xsig_Xsig(&accumulator, &argSignif);
 235         shr_Xsig(&accumulator, 6 - adj);
 236
 237         mul32_Xsig(&arg_signif, leadterm);
 238         add_two_Xsig(&accumulator, &arg_signif, &exponent);
 239
 240         *expon = exponent + 1;
 241         accum_result->lsw = accumulator.lsw;
 242         accum_result->midw = accumulator.midw;
 243         accum_result->msw = accumulator.msw;
 244
 245 }