fp_log.c

  1/*
  2
  3  fp_log.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_arith.h"
 19#include "fp_emu.h"
 20#include "fp_log.h"
 21
 22static const struct fp_ext fp_one = {
 23	.exp = 0x3fff,
 24};
 25
 26struct fp_ext *fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
 27{
 28	struct fp_ext tmp, src2;
 29	int i, exp;
 30
 31	dprint(PINSTR, "fsqrt\n");
 32
 33	fp_monadic_check(dest, src);
 34
 35	if (IS_ZERO(dest))
 36		return dest;
 37
 38	if (dest->sign) {
 39		fp_set_nan(dest);
 40		return dest;
 41	}
 42	if (IS_INF(dest))
 43		return dest;
 44
 45	/*
 46	 *		 sqrt(m) * 2^(p)	, if e = 2*p
 47	 * sqrt(m*2^e) =
 48	 *		 sqrt(2*m) * 2^(p)	, if e = 2*p + 1
 49	 *
 50	 * So we use the last bit of the exponent to decide whether to
 51	 * use the m or 2*m.
 52	 *
 53	 * Since only the fractional part of the mantissa is stored and
 54	 * the integer part is assumed to be one, we place a 1 or 2 into
 55	 * the fixed point representation.
 56	 */
 57	exp = dest->exp;
 58	dest->exp = 0x3FFF;
 59	if (!(exp & 1))		/* lowest bit of exponent is set */
 60		dest->exp++;
 61	fp_copy_ext(&src2, dest);
 62
 63	/*
 64	 * The taylor row around a for sqrt(x) is:
 65	 *	sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
 66	 * With a=1 this gives:
 67	 *	sqrt(x) = 1 + 1/2*(x-1)
 68	 *		= 1/2*(1+x)
 69	 */
 70	/* It is safe to cast away the constness, as fp_one is normalized */
 71	fp_fadd(dest, (struct fp_ext *)&fp_one);
 72	dest->exp--;		/* * 1/2 */
 73
 74	/*
 75	 * We now apply the newton rule to the function
 76	 *	f(x) := x^2 - r
 77	 * which has a null point on x = sqrt(r).
 78	 *
 79	 * It gives:
 80	 *	x' := x - f(x)/f'(x)
 81	 *	    = x - (x^2 -r)/(2*x)
 82	 *	    = x - (x - r/x)/2
 83	 *          = (2*x - x + r/x)/2
 84	 *	    = (x + r/x)/2
 85	 */
 86	for (i = 0; i < 9; i++) {
 87		fp_copy_ext(&tmp, &src2);
 88
 89		fp_fdiv(&tmp, dest);
 90		fp_fadd(dest, &tmp);
 91		dest->exp--;
 92	}
 93
 94	dest->exp += (exp - 0x3FFF) / 2;
 95
 96	return dest;
 97}
 98
 99struct fp_ext *fp_fetoxm1(struct fp_ext *dest, struct fp_ext *src)
100{
101	uprint("fetoxm1\n");
102
103	fp_monadic_check(dest, src);
104
105	return dest;
106}
107
108struct fp_ext *fp_fetox(struct fp_ext *dest, struct fp_ext *src)
109{
110	uprint("fetox\n");
111
112	fp_monadic_check(dest, src);
113
114	return dest;
115}
116
117struct fp_ext *fp_ftwotox(struct fp_ext *dest, struct fp_ext *src)
118{
119	uprint("ftwotox\n");
120
121	fp_monadic_check(dest, src);
122
123	return dest;
124}
125
126struct fp_ext *fp_ftentox(struct fp_ext *dest, struct fp_ext *src)
127{
128	uprint("ftentox\n");
129
130	fp_monadic_check(dest, src);
131
132	return dest;
133}
134
135struct fp_ext *fp_flogn(struct fp_ext *dest, struct fp_ext *src)
136{
137	uprint("flogn\n");
138
139	fp_monadic_check(dest, src);
140
141	return dest;
142}
143
144struct fp_ext *fp_flognp1(struct fp_ext *dest, struct fp_ext *src)
145{
146	uprint("flognp1\n");
147
148	fp_monadic_check(dest, src);
149
150	return dest;
151}
152
153struct fp_ext *fp_flog10(struct fp_ext *dest, struct fp_ext *src)
154{
155	uprint("flog10\n");
156
157	fp_monadic_check(dest, src);
158
159	return dest;
160}
161
162struct fp_ext *fp_flog2(struct fp_ext *dest, struct fp_ext *src)
163{
164	uprint("flog2\n");
165
166	fp_monadic_check(dest, src);
167
168	return dest;
169}
170
171struct fp_ext *fp_fgetexp(struct fp_ext *dest, struct fp_ext *src)
172{
173	dprint(PINSTR, "fgetexp\n");
174
175	fp_monadic_check(dest, src);
176
177	if (IS_INF(dest)) {
178		fp_set_nan(dest);
179		return dest;
180	}
181	if (IS_ZERO(dest))
182		return dest;
183
184	fp_conv_long2ext(dest, (int)dest->exp - 0x3FFF);
185
186	fp_normalize_ext(dest);
187
188	return dest;
189}
190
191struct fp_ext *fp_fgetman(struct fp_ext *dest, struct fp_ext *src)
192{
193	dprint(PINSTR, "fgetman\n");
194
195	fp_monadic_check(dest, src);
196
197	if (IS_ZERO(dest))
198		return dest;
199
200	if (IS_INF(dest))
201		return dest;
202
203	dest->exp = 0x3FFF;
204
205	return dest;
206}
207