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