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);
  27
  28struct fp_ext *
  29fp_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 whether 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
 101struct fp_ext *
 102fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
 103{
 104        uprint("fetoxm1\n");
 105
 106        fp_monadic_check(dest, src);
 107
 108        return dest;
 109}
 110
 111struct fp_ext *
 112fp_fetox(struct fp_ext *dest, struct fp_ext *src)
 113{
 114        uprint("fetox\n");
 115
 116        fp_monadic_check(dest, src);
 117
 118        return dest;
 119}
 120
 121struct fp_ext *
 122fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
 123{
 124        uprint("ftwotox\n");
 125
 126        fp_monadic_check(dest, src);
 127
 128        return dest;
 129}
 130
 131struct fp_ext *
 132fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
 133{
 134        uprint("ftentox\n");
 135
 136        fp_monadic_check(dest, src);
 137
 138        return dest;
 139}
 140
 141struct fp_ext *
 142fp_flogn(struct fp_ext *dest, struct fp_ext *src)
 143{
 144        uprint("flogn\n");
 145
 146        fp_monadic_check(dest, src);
 147
 148        return dest;
 149}
 150
 151struct fp_ext *
 152fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
 153{
 154        uprint("flognp1\n");
 155
 156        fp_monadic_check(dest, src);
 157
 158        return dest;
 159}
 160
 161struct fp_ext *
 162fp_flog10(struct fp_ext *dest, struct fp_ext *src)
 163{
 164        uprint("flog10\n");
 165
 166        fp_monadic_check(dest, src);
 167
 168        return dest;
 169}
 170
 171struct fp_ext *
 172fp_flog2(struct fp_ext *dest, struct fp_ext *src)
 173{
 174        uprint("flog2\n");
 175
 176        fp_monadic_check(dest, src);
 177
 178        return dest;
 179}
 180
 181struct fp_ext *
 182fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
 183{
 184        dprint(PINSTR, "fgetexp\n");
 185
 186        fp_monadic_check(dest, src);
 187
 188        if (IS_INF(dest)) {
 189                fp_set_nan(dest);
 190                return dest;
 191        }
 192        if (IS_ZERO(dest))
 193                return dest;
 194
 195        fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
 196
 197        fp_normalize_ext(dest);
 198
 199        return dest;
 200}
 201
 202struct fp_ext *
 203fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
 204{
 205        dprint(PINSTR, "fgetman\n");
 206
 207        fp_monadic_check(dest, src);
 208
 209        if (IS_ZERO(dest))
 210                return dest;
 211
 212        if (IS_INF(dest))
 213                return dest;
 214
 215        dest->exp = 0x3FFF;
 216
 217        return dest;
 218}
 219
 220