linux/arch/m68k/math-emu/fp_log.c
<<
>>
Prefs
   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
  20static const struct fp_ext fp_one =
  21{
  22        .exp = 0x3fff,
  23};
  24
  25extern struct fp_ext *fp_fadd(struct fp_ext *dest, const struct fp_ext *src);
  26extern struct fp_ext *fp_fdiv(struct fp_ext *dest, const struct fp_ext *src);
  27extern struct fp_ext *fp_fmul(struct fp_ext *dest, const struct fp_ext *src);
  28
  29struct fp_ext *
  30fp_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 around 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
 102struct fp_ext *
 103fp_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
 115struct fp_ext *
 116fp_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
 125struct fp_ext *
 126fp_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
 135struct fp_ext *
 136fp_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
 145struct fp_ext *
 146fp_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
 155struct fp_ext *
 156fp_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
 165struct fp_ext *
 166fp_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
 175struct fp_ext *
 176fp_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
 185struct fp_ext *
 186fp_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
 206struct fp_ext *
 207fp_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
 224