arch/m68k/math-emu/fp_log.c

   1 /*
   2
   3   fp_trig.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_emu.h"
  19
  20 static const struct fp_ext fp_one =
  21 {
  22         0, 0, 0x3fff, { 0 }
  23 };
  24
  25 extern struct fp_ext *fp_fadd(struct fp_ext *dest, const struct fp_ext *src);
  26 extern struct fp_ext *fp_fdiv(struct fp_ext *dest, const struct fp_ext *src);
  27 extern struct fp_ext *fp_fmul(struct fp_ext *dest, const struct fp_ext *src);
  28
  29 struct fp_ext *
  30 fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
  31 {
  32         struct fp_ext tmp, src2;
  33         int i, exp;
  34
  35         dprint(PINSTR, "fsqrt\n");
  36
  37         fp_monadic_check(dest, src);
  38
  39         if (IS_ZERO(dest))
  40                 return dest;
  41
  42         if (dest->sign) {
  43                 fp_set_nan(dest);
  44                 return dest;
  45         }
  46         if (IS_INF(dest))
  47                 return dest;
  48
  49         /*
  50          *               sqrt(m) * 2^(p)        , if e = 2*p
  51          * sqrt(m*2^e) =
  52          *               sqrt(2*m) * 2^(p)      , if e = 2*p + 1
  53          *
  54          * So we use the last bit of the exponent to decide wether to
  55          * use the m or 2*m.
  56          *
  57          * Since only the fractional part of the mantissa is stored and
  58          * the integer part is assumed to be one, we place a 1 or 2 into
  59          * the fixed point representation.
  60          */
  61         exp = dest->exp;
  62         dest->exp = 0x3FFF;
  63         if (!(exp & 1))         /* lowest bit of exponent is set */
  64                 dest->exp++;
  65         fp_copy_ext(&src2, dest);
  66
  67         /*
  68          * The taylor row arround a for sqrt(x) is:
  69          *      sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
  70          * With a=1 this gives:
  71          *      sqrt(x) = 1 + 1/2*(x-1)
  72          *              = 1/2*(1+x)
  73          */
  74         fp_fadd(dest, &fp_one);
  75         dest->exp--;            /* * 1/2 */
  76
  77         /*
  78          * We now apply the newton rule to the function
  79          *      f(x) := x^2 - r
  80          * which has a null point on x = sqrt(r).
  81          *
  82          * It gives:
  83          *      x' := x - f(x)/f'(x)
  84          *          = x - (x^2 -r)/(2*x)
  85          *          = x - (x - r/x)/2
  86          *          = (2*x - x + r/x)/2
  87          *          = (x + r/x)/2
  88          */
  89         for (i = 0; i < 9; i++) {
  90                 fp_copy_ext(&tmp, &src2);
  91
  92                 fp_fdiv(&tmp, dest);
  93                 fp_fadd(dest, &tmp);
  94                 dest->exp--;
  95         }
  96
  97         dest->exp += (exp - 0x3FFF) / 2;
  98
  99         return dest;
 100 }
 101
 102 struct fp_ext *
 103 fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
 104 {
 105         uprint("fetoxm1\n");
 106
 107         fp_monadic_check(dest, src);
 108
 109         if (IS_ZERO(dest))
 110                 return dest;
 111
 112         return dest;
 113 }
 114
 115 struct fp_ext *
 116 fp_fetox(struct fp_ext *dest, struct fp_ext *src)
 117 {
 118         uprint("fetox\n");
 119
 120         fp_monadic_check(dest, src);
 121
 122         return dest;
 123 }
 124
 125 struct fp_ext *
 126 fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
 127 {
 128         uprint("ftwotox\n");
 129
 130         fp_monadic_check(dest, src);
 131
 132         return dest;
 133 }
 134
 135 struct fp_ext *
 136 fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
 137 {
 138         uprint("ftentox\n");
 139
 140         fp_monadic_check(dest, src);
 141
 142         return dest;
 143 }
 144
 145 struct fp_ext *
 146 fp_flogn(struct fp_ext *dest, struct fp_ext *src)
 147 {
 148         uprint("flogn\n");
 149
 150         fp_monadic_check(dest, src);
 151
 152         return dest;
 153 }
 154
 155 struct fp_ext *
 156 fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
 157 {
 158         uprint("flognp1\n");
 159
 160         fp_monadic_check(dest, src);
 161
 162         return dest;
 163 }
 164
 165 struct fp_ext *
 166 fp_flog10(struct fp_ext *dest, struct fp_ext *src)
 167 {
 168         uprint("flog10\n");
 169
 170         fp_monadic_check(dest, src);
 171
 172         return dest;
 173 }
 174
 175 struct fp_ext *
 176 fp_flog2(struct fp_ext *dest, struct fp_ext *src)
 177 {
 178         uprint("flog2\n");
 179
 180         fp_monadic_check(dest, src);
 181
 182         return dest;
 183 }
 184
 185 struct fp_ext *
 186 fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
 187 {
 188         dprint(PINSTR, "fgetexp\n");
 189
 190         fp_monadic_check(dest, src);
 191
 192         if (IS_INF(dest)) {
 193                 fp_set_nan(dest);
 194                 return dest;
 195         }
 196         if (IS_ZERO(dest))
 197                 return dest;
 198
 199         fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
 200
 201         fp_normalize_ext(dest);
 202
 203         return dest;
 204 }
 205
 206 struct fp_ext *
 207 fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
 208 {
 209         dprint(PINSTR, "fgetman\n");
 210
 211         fp_monadic_check(dest, src);
 212
 213         if (IS_ZERO(dest))
 214                 return dest;
 215
 216         if (IS_INF(dest))
 217                 return dest;
 218
 219         dest->exp = 0x3FFF;
 220
 221         return dest;
 222 }
 223