arch/m68k/math-emu/fp_log.c

   1 /*
   2
   3   fp_log.c: floating-point math routines for the Linux-m68k
   4   floating point emulator.
   5
   6   Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.
   7
   8   I hereby give permission, free of charge, to copy, modify, and
   9   redistribute this software, in source or binary form, provided that
  10   the above copyright notice and the following disclaimer are included
  11   in all such copies.
  12
  13   THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
  14   OR IMPLIED.
  15
  16 */
  17
  18 #include "fp_arith.h"
  19 #include "fp_emu.h"
  20 #include "fp_log.h"
  21
  22 static const struct fp_ext fp_one = {
  23         .exp = 0x3fff,
  24 };
  25
  26 struct fp_ext *fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
  27 {
  28         struct fp_ext tmp, src2;
  29         int i, exp;
  30
  31         dprint(PINSTR, "fsqrt\n");
  32
  33         fp_monadic_check(dest, src);
  34
  35         if (IS_ZERO(dest))
  36                 return dest;
  37
  38         if (dest->sign) {
  39                 fp_set_nan(dest);
  40                 return dest;
  41         }
  42         if (IS_INF(dest))
  43                 return dest;
  44
  45         /*
  46          *               sqrt(m) * 2^(p)        , if e = 2*p
  47          * sqrt(m*2^e) =
  48          *               sqrt(2*m) * 2^(p)      , if e = 2*p + 1
  49          *
  50          * So we use the last bit of the exponent to decide whether to
  51          * use the m or 2*m.
  52          *
  53          * Since only the fractional part of the mantissa is stored and
  54          * the integer part is assumed to be one, we place a 1 or 2 into
  55          * the fixed point representation.
  56          */
  57         exp = dest->exp;
  58         dest->exp = 0x3FFF;
  59         if (!(exp & 1))         /* lowest bit of exponent is set */
  60                 dest->exp++;
  61         fp_copy_ext(&src2, dest);
  62
  63         /*
  64          * The taylor row around a for sqrt(x) is:
  65          *      sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
  66          * With a=1 this gives:
  67          *      sqrt(x) = 1 + 1/2*(x-1)
  68          *              = 1/2*(1+x)
  69          */
  70         /* It is safe to cast away the constness, as fp_one is normalized */
  71         fp_fadd(dest, (struct fp_ext *)&fp_one);
  72         dest->exp--;            /* * 1/2 */
  73
  74         /*
  75          * We now apply the newton rule to the function
  76          *      f(x) := x^2 - r
  77          * which has a null point on x = sqrt(r).
  78          *
  79          * It gives:
  80          *      x' := x - f(x)/f'(x)
  81          *          = x - (x^2 -r)/(2*x)
  82          *          = x - (x - r/x)/2
  83          *          = (2*x - x + r/x)/2
  84          *          = (x + r/x)/2
  85          */
  86         for (i = 0; i < 9; i++) {
  87                 fp_copy_ext(&tmp, &src2);
  88
  89                 fp_fdiv(&tmp, dest);
  90                 fp_fadd(dest, &tmp);
  91                 dest->exp--;
  92         }
  93
  94         dest->exp += (exp - 0x3FFF) / 2;
  95
  96         return dest;
  97 }
  98
  99 struct fp_ext *fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
 100 {
 101         uprint("fetoxm1\n");
 102
 103         fp_monadic_check(dest, src);
 104
 105         return dest;
 106 }
 107
 108 struct fp_ext *fp_fetox(struct fp_ext *dest, struct fp_ext *src)
 109 {
 110         uprint("fetox\n");
 111
 112         fp_monadic_check(dest, src);
 113
 114         return dest;
 115 }
 116
 117 struct fp_ext *fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
 118 {
 119         uprint("ftwotox\n");
 120
 121         fp_monadic_check(dest, src);
 122
 123         return dest;
 124 }
 125
 126 struct fp_ext *fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
 127 {
 128         uprint("ftentox\n");
 129
 130         fp_monadic_check(dest, src);
 131
 132         return dest;
 133 }
 134
 135 struct fp_ext *fp_flogn(struct fp_ext *dest, struct fp_ext *src)
 136 {
 137         uprint("flogn\n");
 138
 139         fp_monadic_check(dest, src);
 140
 141         return dest;
 142 }
 143
 144 struct fp_ext *fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
 145 {
 146         uprint("flognp1\n");
 147
 148         fp_monadic_check(dest, src);
 149
 150         return dest;
 151 }
 152
 153 struct fp_ext *fp_flog10(struct fp_ext *dest, struct fp_ext *src)
 154 {
 155         uprint("flog10\n");
 156
 157         fp_monadic_check(dest, src);
 158
 159         return dest;
 160 }
 161
 162 struct fp_ext *fp_flog2(struct fp_ext *dest, struct fp_ext *src)
 163 {
 164         uprint("flog2\n");
 165
 166         fp_monadic_check(dest, src);
 167
 168         return dest;
 169 }
 170
 171 struct fp_ext *fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
 172 {
 173         dprint(PINSTR, "fgetexp\n");
 174
 175         fp_monadic_check(dest, src);
 176
 177         if (IS_INF(dest)) {
 178                 fp_set_nan(dest);
 179                 return dest;
 180         }
 181         if (IS_ZERO(dest))
 182                 return dest;
 183
 184         fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
 185
 186         fp_normalize_ext(dest);
 187
 188         return dest;
 189 }
 190
 191 struct fp_ext *fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
 192 {
 193         dprint(PINSTR, "fgetman\n");
 194
 195         fp_monadic_check(dest, src);
 196
 197         if (IS_ZERO(dest))
 198                 return dest;
 199
 200         if (IS_INF(dest))
 201                 return dest;
 202
 203         dest->exp = 0x3FFF;
 204
 205         return dest;
 206 }
 207