Loading...
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
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