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 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 
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 	return dest;
109 }
110 
111 struct fp_ext *
fp_fetox(struct fp_ext * dest,struct fp_ext * src)112 fp_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 
121 struct fp_ext *
fp_ftwotox(struct fp_ext * dest,struct fp_ext * src)122 fp_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 
131 struct fp_ext *
fp_ftentox(struct fp_ext * dest,struct fp_ext * src)132 fp_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 
141 struct fp_ext *
fp_flogn(struct fp_ext * dest,struct fp_ext * src)142 fp_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 
151 struct fp_ext *
fp_flognp1(struct fp_ext * dest,struct fp_ext * src)152 fp_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 
161 struct fp_ext *
fp_flog10(struct fp_ext * dest,struct fp_ext * src)162 fp_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 
171 struct fp_ext *
fp_flog2(struct fp_ext * dest,struct fp_ext * src)172 fp_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 
181 struct fp_ext *
fp_fgetexp(struct fp_ext * dest,struct fp_ext * src)182 fp_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 
202 struct fp_ext *
fp_fgetman(struct fp_ext * dest,struct fp_ext * src)203 fp_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