1 | /******************************************************************\ |
---|
2 | * * |
---|
3 | * <math-68881.h> last modified: 18 May 1989. * |
---|
4 | * * |
---|
5 | * Copyright (C) 1989 by Matthew Self. * |
---|
6 | * You may freely distribute verbatim copies of this software * |
---|
7 | * provided that this copyright notice is retained in all copies. * |
---|
8 | * You may distribute modifications to this software under the * |
---|
9 | * conditions above if you also clearly note such modifications * |
---|
10 | * with their author and date. * |
---|
11 | * * |
---|
12 | * Note: errno is not set to EDOM when domain errors occur for * |
---|
13 | * most of these functions. Rather, it is assumed that the * |
---|
14 | * 68881's OPERR exception will be enabled and handled * |
---|
15 | * appropriately by the operating system. Similarly, overflow * |
---|
16 | * and underflow do not set errno to ERANGE. * |
---|
17 | * * |
---|
18 | * Send bugs to Matthew Self (self@bayes.arc.nasa.gov). * |
---|
19 | * * |
---|
20 | \******************************************************************/ |
---|
21 | |
---|
22 | #include <errno.h> |
---|
23 | |
---|
24 | #undef HUGE_VAL |
---|
25 | #define HUGE_VAL \ |
---|
26 | ({ \ |
---|
27 | double huge_val; \ |
---|
28 | \ |
---|
29 | __asm ("fmove%.d %#0x7ff0000000000000,%0" /* Infinity */ \ |
---|
30 | : "=f" (huge_val) \ |
---|
31 | : /* no inputs */); \ |
---|
32 | huge_val; \ |
---|
33 | }) |
---|
34 | |
---|
35 | __inline static const double sin (double x) |
---|
36 | { |
---|
37 | double value; |
---|
38 | |
---|
39 | __asm ("fsin%.x %1,%0" |
---|
40 | : "=f" (value) |
---|
41 | : "f" (x)); |
---|
42 | return value; |
---|
43 | } |
---|
44 | |
---|
45 | __inline static const double cos (double x) |
---|
46 | { |
---|
47 | double value; |
---|
48 | |
---|
49 | __asm ("fcos%.x %1,%0" |
---|
50 | : "=f" (value) |
---|
51 | : "f" (x)); |
---|
52 | return value; |
---|
53 | } |
---|
54 | |
---|
55 | __inline static const double tan (double x) |
---|
56 | { |
---|
57 | double value; |
---|
58 | |
---|
59 | __asm ("ftan%.x %1,%0" |
---|
60 | : "=f" (value) |
---|
61 | : "f" (x)); |
---|
62 | return value; |
---|
63 | } |
---|
64 | |
---|
65 | __inline static const double asin (double x) |
---|
66 | { |
---|
67 | double value; |
---|
68 | |
---|
69 | __asm ("fasin%.x %1,%0" |
---|
70 | : "=f" (value) |
---|
71 | : "f" (x)); |
---|
72 | return value; |
---|
73 | } |
---|
74 | |
---|
75 | __inline static const double acos (double x) |
---|
76 | { |
---|
77 | double value; |
---|
78 | |
---|
79 | __asm ("facos%.x %1,%0" |
---|
80 | : "=f" (value) |
---|
81 | : "f" (x)); |
---|
82 | return value; |
---|
83 | } |
---|
84 | |
---|
85 | __inline static const double atan (double x) |
---|
86 | { |
---|
87 | double value; |
---|
88 | |
---|
89 | __asm ("fatan%.x %1,%0" |
---|
90 | : "=f" (value) |
---|
91 | : "f" (x)); |
---|
92 | return value; |
---|
93 | } |
---|
94 | |
---|
95 | __inline static const double atan2 (double y, double x) |
---|
96 | { |
---|
97 | double pi, pi_over_2; |
---|
98 | |
---|
99 | __asm ("fmovecr%.x %#0,%0" /* extended precision pi */ |
---|
100 | : "=f" (pi) |
---|
101 | : /* no inputs */ ); |
---|
102 | __asm ("fscale%.b %#-1,%0" /* no loss of accuracy */ |
---|
103 | : "=f" (pi_over_2) |
---|
104 | : "0" (pi)); |
---|
105 | if (x > 0) |
---|
106 | { |
---|
107 | if (y > 0) |
---|
108 | { |
---|
109 | if (x > y) |
---|
110 | return atan (y / x); |
---|
111 | else |
---|
112 | return pi_over_2 - atan (x / y); |
---|
113 | } |
---|
114 | else |
---|
115 | { |
---|
116 | if (x > -y) |
---|
117 | return atan (y / x); |
---|
118 | else |
---|
119 | return - pi_over_2 - atan (x / y); |
---|
120 | } |
---|
121 | } |
---|
122 | else |
---|
123 | { |
---|
124 | if (y > 0) |
---|
125 | { |
---|
126 | if (-x > y) |
---|
127 | return pi + atan (y / x); |
---|
128 | else |
---|
129 | return pi_over_2 - atan (x / y); |
---|
130 | } |
---|
131 | else |
---|
132 | { |
---|
133 | if (-x > -y) |
---|
134 | return - pi + atan (y / x); |
---|
135 | else if (y < 0) |
---|
136 | return - pi_over_2 - atan (x / y); |
---|
137 | else |
---|
138 | { |
---|
139 | double value; |
---|
140 | |
---|
141 | errno = EDOM; |
---|
142 | __asm ("fmove%.d %#0x7fffffffffffffff,%0" /* quiet NaN */ |
---|
143 | : "=f" (value) |
---|
144 | : /* no inputs */); |
---|
145 | return value; |
---|
146 | } |
---|
147 | } |
---|
148 | } |
---|
149 | } |
---|
150 | |
---|
151 | __inline static const double sinh (double x) |
---|
152 | { |
---|
153 | double value; |
---|
154 | |
---|
155 | __asm ("fsinh%.x %1,%0" |
---|
156 | : "=f" (value) |
---|
157 | : "f" (x)); |
---|
158 | return value; |
---|
159 | } |
---|
160 | |
---|
161 | __inline static const double cosh (double x) |
---|
162 | { |
---|
163 | double value; |
---|
164 | |
---|
165 | __asm ("fcosh%.x %1,%0" |
---|
166 | : "=f" (value) |
---|
167 | : "f" (x)); |
---|
168 | return value; |
---|
169 | } |
---|
170 | |
---|
171 | __inline static const double tanh (double x) |
---|
172 | { |
---|
173 | double value; |
---|
174 | |
---|
175 | __asm ("ftanh%.x %1,%0" |
---|
176 | : "=f" (value) |
---|
177 | : "f" (x)); |
---|
178 | return value; |
---|
179 | } |
---|
180 | |
---|
181 | __inline static const double atanh (double x) |
---|
182 | { |
---|
183 | double value; |
---|
184 | |
---|
185 | __asm ("fatanh%.x %1,%0" |
---|
186 | : "=f" (value) |
---|
187 | : "f" (x)); |
---|
188 | return value; |
---|
189 | } |
---|
190 | |
---|
191 | __inline static const double exp (double x) |
---|
192 | { |
---|
193 | double value; |
---|
194 | |
---|
195 | __asm ("fetox%.x %1,%0" |
---|
196 | : "=f" (value) |
---|
197 | : "f" (x)); |
---|
198 | return value; |
---|
199 | } |
---|
200 | |
---|
201 | __inline static const double expm1 (double x) |
---|
202 | { |
---|
203 | double value; |
---|
204 | |
---|
205 | __asm ("fetoxm1%.x %1,%0" |
---|
206 | : "=f" (value) |
---|
207 | : "f" (x)); |
---|
208 | return value; |
---|
209 | } |
---|
210 | |
---|
211 | __inline static const double log (double x) |
---|
212 | { |
---|
213 | double value; |
---|
214 | |
---|
215 | __asm ("flogn%.x %1,%0" |
---|
216 | : "=f" (value) |
---|
217 | : "f" (x)); |
---|
218 | return value; |
---|
219 | } |
---|
220 | |
---|
221 | __inline static const double log1p (double x) |
---|
222 | { |
---|
223 | double value; |
---|
224 | |
---|
225 | __asm ("flognp1%.x %1,%0" |
---|
226 | : "=f" (value) |
---|
227 | : "f" (x)); |
---|
228 | return value; |
---|
229 | } |
---|
230 | |
---|
231 | __inline static const double log10 (double x) |
---|
232 | { |
---|
233 | double value; |
---|
234 | |
---|
235 | __asm ("flog10%.x %1,%0" |
---|
236 | : "=f" (value) |
---|
237 | : "f" (x)); |
---|
238 | return value; |
---|
239 | } |
---|
240 | |
---|
241 | __inline static const double sqrt (double x) |
---|
242 | { |
---|
243 | double value; |
---|
244 | |
---|
245 | __asm ("fsqrt%.x %1,%0" |
---|
246 | : "=f" (value) |
---|
247 | : "f" (x)); |
---|
248 | return value; |
---|
249 | } |
---|
250 | |
---|
251 | __inline static const double pow (const double x, const double y) |
---|
252 | { |
---|
253 | if (x > 0) |
---|
254 | return exp (y * log (x)); |
---|
255 | else if (x == 0) |
---|
256 | { |
---|
257 | if (y > 0) |
---|
258 | return 0.0; |
---|
259 | else |
---|
260 | { |
---|
261 | double value; |
---|
262 | |
---|
263 | errno = EDOM; |
---|
264 | __asm ("fmove%.d %#0x7fffffffffffffff,%0" /* quiet NaN */ |
---|
265 | : "=f" (value) |
---|
266 | : /* no inputs */); |
---|
267 | return value; |
---|
268 | } |
---|
269 | } |
---|
270 | else |
---|
271 | { |
---|
272 | double temp; |
---|
273 | |
---|
274 | __asm ("fintrz%.x %1,%0" |
---|
275 | : "=f" (temp) /* integer-valued float */ |
---|
276 | : "f" (y)); |
---|
277 | if (y == temp) |
---|
278 | { |
---|
279 | int i = (int) y; |
---|
280 | |
---|
281 | if (i & 1 == 0) /* even */ |
---|
282 | return exp (y * log (x)); |
---|
283 | else |
---|
284 | return - exp (y * log (x)); |
---|
285 | } |
---|
286 | else |
---|
287 | { |
---|
288 | double value; |
---|
289 | |
---|
290 | errno = EDOM; |
---|
291 | __asm ("fmove%.d %#0x7fffffffffffffff,%0" /* quiet NaN */ |
---|
292 | : "=f" (value) |
---|
293 | : /* no inputs */); |
---|
294 | return value; |
---|
295 | } |
---|
296 | } |
---|
297 | } |
---|
298 | |
---|
299 | __inline static const double fabs (double x) |
---|
300 | { |
---|
301 | double value; |
---|
302 | |
---|
303 | __asm ("fabs%.x %1,%0" |
---|
304 | : "=f" (value) |
---|
305 | : "f" (x)); |
---|
306 | return value; |
---|
307 | } |
---|
308 | |
---|
309 | __inline static const double ceil (double x) |
---|
310 | { |
---|
311 | int rounding_mode, round_up; |
---|
312 | double value; |
---|
313 | |
---|
314 | __asm volatile ("fmove%.l %%fpcr,%0" |
---|
315 | : "=dm" (rounding_mode) |
---|
316 | : /* no inputs */ ); |
---|
317 | round_up = rounding_mode | 0x30; |
---|
318 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
319 | : /* no outputs */ |
---|
320 | : "dmi" (round_up)); |
---|
321 | __asm volatile ("fint%.x %1,%0" |
---|
322 | : "=f" (value) |
---|
323 | : "f" (x)); |
---|
324 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
325 | : /* no outputs */ |
---|
326 | : "dmi" (rounding_mode)); |
---|
327 | return value; |
---|
328 | } |
---|
329 | |
---|
330 | __inline static const double floor (double x) |
---|
331 | { |
---|
332 | int rounding_mode, round_down; |
---|
333 | double value; |
---|
334 | |
---|
335 | __asm volatile ("fmove%.l %%fpcr,%0" |
---|
336 | : "=dm" (rounding_mode) |
---|
337 | : /* no inputs */ ); |
---|
338 | round_down = (rounding_mode & ~0x10) |
---|
339 | | 0x20; |
---|
340 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
341 | : /* no outputs */ |
---|
342 | : "dmi" (round_down)); |
---|
343 | __asm volatile ("fint%.x %1,%0" |
---|
344 | : "=f" (value) |
---|
345 | : "f" (x)); |
---|
346 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
347 | : /* no outputs */ |
---|
348 | : "dmi" (rounding_mode)); |
---|
349 | return value; |
---|
350 | } |
---|
351 | |
---|
352 | __inline static const double rint (double x) |
---|
353 | { |
---|
354 | int rounding_mode, round_nearest; |
---|
355 | double value; |
---|
356 | |
---|
357 | __asm volatile ("fmove%.l %%fpcr,%0" |
---|
358 | : "=dm" (rounding_mode) |
---|
359 | : /* no inputs */ ); |
---|
360 | round_nearest = rounding_mode & ~0x30; |
---|
361 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
362 | : /* no outputs */ |
---|
363 | : "dmi" (round_nearest)); |
---|
364 | __asm volatile ("fint%.x %1,%0" |
---|
365 | : "=f" (value) |
---|
366 | : "f" (x)); |
---|
367 | __asm volatile ("fmove%.l %0,%%fpcr" |
---|
368 | : /* no outputs */ |
---|
369 | : "dmi" (rounding_mode)); |
---|
370 | return value; |
---|
371 | } |
---|
372 | |
---|
373 | __inline static const double fmod (double x, double y) |
---|
374 | { |
---|
375 | double value; |
---|
376 | |
---|
377 | __asm ("fmod%.x %2,%0" |
---|
378 | : "=f" (value) |
---|
379 | : "0" (x), |
---|
380 | "f" (y)); |
---|
381 | return value; |
---|
382 | } |
---|
383 | |
---|
384 | __inline static const double drem (double x, double y) |
---|
385 | { |
---|
386 | double value; |
---|
387 | |
---|
388 | __asm ("frem%.x %2,%0" |
---|
389 | : "=f" (value) |
---|
390 | : "0" (x), |
---|
391 | "f" (y)); |
---|
392 | return value; |
---|
393 | } |
---|
394 | |
---|
395 | __inline static const double scalb (double x, int n) |
---|
396 | { |
---|
397 | double value; |
---|
398 | |
---|
399 | __asm ("fscale%.l %2,%0" |
---|
400 | : "=f" (value) |
---|
401 | : "0" (x), |
---|
402 | "dmi" (n)); |
---|
403 | return value; |
---|
404 | } |
---|
405 | |
---|
406 | __inline static double logb (double x) |
---|
407 | { |
---|
408 | double exponent; |
---|
409 | |
---|
410 | __asm ("fgetexp%.x %1,%0" |
---|
411 | : "=f" (exponent) |
---|
412 | : "f" (x)); |
---|
413 | return exponent; |
---|
414 | } |
---|
415 | |
---|
416 | __inline static const double ldexp (double x, int n) |
---|
417 | { |
---|
418 | double value; |
---|
419 | |
---|
420 | __asm ("fscale%.l %2,%0" |
---|
421 | : "=f" (value) |
---|
422 | : "0" (x), |
---|
423 | "dmi" (n)); |
---|
424 | return value; |
---|
425 | } |
---|
426 | |
---|
427 | __inline static double frexp (double x, int *exp) |
---|
428 | { |
---|
429 | double float_exponent; |
---|
430 | int int_exponent; |
---|
431 | double mantissa; |
---|
432 | |
---|
433 | __asm ("fgetexp%.x %1,%0" |
---|
434 | : "=f" (float_exponent) /* integer-valued float */ |
---|
435 | : "f" (x)); |
---|
436 | int_exponent = (int) float_exponent; |
---|
437 | __asm ("fgetman%.x %1,%0" |
---|
438 | : "=f" (mantissa) /* 1.0 <= mantissa < 2.0 */ |
---|
439 | : "f" (x)); |
---|
440 | if (mantissa != 0) |
---|
441 | { |
---|
442 | __asm ("fscale%.b %#-1,%0" |
---|
443 | : "=f" (mantissa) /* mantissa /= 2.0 */ |
---|
444 | : "0" (mantissa)); |
---|
445 | int_exponent += 1; |
---|
446 | } |
---|
447 | *exp = int_exponent; |
---|
448 | return mantissa; |
---|
449 | } |
---|
450 | |
---|
451 | __inline static double modf (double x, double *ip) |
---|
452 | { |
---|
453 | double temp; |
---|
454 | |
---|
455 | __asm ("fintrz%.x %1,%0" |
---|
456 | : "=f" (temp) /* integer-valued float */ |
---|
457 | : "f" (x)); |
---|
458 | *ip = temp; |
---|
459 | return x - temp; |
---|
460 | } |
---|
461 | |
---|