b1033ae0d6f084a357ddf318a8e67b10d789098a
[linux-3.10.git] / 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 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 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
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