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 *
fp_fsqrt(struct fp_ext * dest,struct fp_ext * src)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 *
fp_fetoxm1(struct fp_ext * dest,struct fp_ext * src)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 *
fp_fetox(struct fp_ext * dest,struct fp_ext * src)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 *
fp_ftwotox(struct fp_ext * dest,struct fp_ext * src)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 *
fp_ftentox(struct fp_ext * dest,struct fp_ext * src)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 *
fp_flogn(struct fp_ext * dest,struct fp_ext * src)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 *
fp_flognp1(struct fp_ext * dest,struct fp_ext * src)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 *
fp_flog10(struct fp_ext * dest,struct fp_ext * src)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 *
fp_flog2(struct fp_ext * dest,struct fp_ext * src)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 *
fp_fgetexp(struct fp_ext * dest,struct fp_ext * src)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 *
fp_fgetman(struct fp_ext * dest,struct fp_ext * src)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
223