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