PageRenderTime 15ms CodeModel.GetById 9ms app.highlight 4ms RepoModel.GetById 1ms app.codeStats 0ms

/arch/m68k/math-emu/fp_log.c

https://bitbucket.org/thekraven/iscream_thunderc-2.6.35
C | 222 lines | 123 code | 52 blank | 47 comment | 10 complexity | 63f7ec655c4d9abddc3fedab14f1c329 MD5 | raw file
Possible License(s): GPL-2.0, LGPL-2.0, AGPL-1.0
  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 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
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	if (IS_ZERO(dest))
109		return dest;
110
111	return dest;
112}
113
114struct fp_ext *
115fp_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
124struct fp_ext *
125fp_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
134struct fp_ext *
135fp_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
144struct fp_ext *
145fp_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
154struct fp_ext *
155fp_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
164struct fp_ext *
165fp_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
174struct fp_ext *
175fp_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
184struct fp_ext *
185fp_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
205struct fp_ext *
206fp_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