3 fp_trig.c: floating-point math routines for the Linux-m68k
4 floating point emulator.
6 Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.
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
13 THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
20 static const struct fp_ext fp_one
=
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 extern struct fp_ext
*fp_fmul(struct fp_ext
*dest
, const struct fp_ext
*src
);
30 fp_fsqrt(struct fp_ext
*dest
, struct fp_ext
*src
)
32 struct fp_ext tmp
, src2
;
35 dprint(PINSTR
, "fsqrt\n");
37 fp_monadic_check(dest
, src
);
50 * sqrt(m) * 2^(p) , if e = 2*p
52 * sqrt(2*m) * 2^(p) , if e = 2*p + 1
54 * So we use the last bit of the exponent to decide wether to
57 * Since only the fractional part of the mantissa is stored and
58 * the integer part is assumed to be one, we place a 1 or 2 into
59 * the fixed point representation.
63 if (!(exp
& 1)) /* lowest bit of exponent is set */
65 fp_copy_ext(&src2
, dest
);
68 * The taylor row around a for sqrt(x) is:
69 * sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
70 * With a=1 this gives:
71 * sqrt(x) = 1 + 1/2*(x-1)
74 fp_fadd(dest
, &fp_one
);
75 dest
->exp
--; /* * 1/2 */
78 * We now apply the newton rule to the function
80 * which has a null point on x = sqrt(r).
83 * x' := x - f(x)/f'(x)
84 * = x - (x^2 -r)/(2*x)
89 for (i
= 0; i
< 9; i
++) {
90 fp_copy_ext(&tmp
, &src2
);
97 dest
->exp
+= (exp
- 0x3FFF) / 2;
103 fp_fetoxm1(struct fp_ext
*dest
, struct fp_ext
*src
)
107 fp_monadic_check(dest
, src
);
116 fp_fetox(struct fp_ext
*dest
, struct fp_ext
*src
)
120 fp_monadic_check(dest
, src
);
126 fp_ftwotox(struct fp_ext
*dest
, struct fp_ext
*src
)
130 fp_monadic_check(dest
, src
);
136 fp_ftentox(struct fp_ext
*dest
, struct fp_ext
*src
)
140 fp_monadic_check(dest
, src
);
146 fp_flogn(struct fp_ext
*dest
, struct fp_ext
*src
)
150 fp_monadic_check(dest
, src
);
156 fp_flognp1(struct fp_ext
*dest
, struct fp_ext
*src
)
160 fp_monadic_check(dest
, src
);
166 fp_flog10(struct fp_ext
*dest
, struct fp_ext
*src
)
170 fp_monadic_check(dest
, src
);
176 fp_flog2(struct fp_ext
*dest
, struct fp_ext
*src
)
180 fp_monadic_check(dest
, src
);
186 fp_fgetexp(struct fp_ext
*dest
, struct fp_ext
*src
)
188 dprint(PINSTR
, "fgetexp\n");
190 fp_monadic_check(dest
, src
);
199 fp_conv_long2ext(dest
, (int)dest
->exp
- 0x3FFF);
201 fp_normalize_ext(dest
);
207 fp_fgetman(struct fp_ext
*dest
, struct fp_ext
*src
)
209 dprint(PINSTR
, "fgetman\n");
211 fp_monadic_check(dest
, src
);