1 /******************************************************************\
2 * *
3 * <math-68881.h> last modified: 23 May 1992. *
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 /* This file is NOT a part of GCC, just distributed with it. */
23
24 /* If you find this in GCC,
25 please send bug reports to bug-gcc@prep.ai.mit.edu. */
26
27 /* Changed by Richard Stallman:
28 May 1993, add conditional to prevent multiple inclusion.
29 % inserted before a #.
30 New function `hypot' added.
31 Nans written in hex to avoid 0rnan.
32 May 1992, use %! for fpcr register. Break lines before function names.
33 December 1989, add parens around `&' in pow.
34 November 1990, added alternate definition of HUGE_VAL for Sun. */
35
36 /* Changed by Jim Wilson:
37 September 1993, Use #undef before HUGE_VAL instead of #ifdef/#endif. */
38
39 /* Changed by Ian Lance Taylor:
40 September 1994, use extern inline instead of static inline. */
41
42 #ifndef __math_68881
43 #define __math_68881
44
45 #include <errno.h>
46
47 #undef HUGE_VAL
48 #ifdef __sun__
49 /* The Sun assembler fails to handle the hex constant in the usual defn. */
50 #define HUGE_VAL \
51 ({ \
52 static union { int i[2]; double d; } u = { {0x7ff00000, 0} }; \
53 u.d; \
54 })
55 #else
56 #define HUGE_VAL \
57 ({ \
58 double huge_val; \
59 \
60 __asm ("fmove%.d #0x7ff0000000000000,%0" /* Infinity */ \
61 : "=f" (huge_val) \
62 : /* no inputs */); \
63 huge_val; \
64 })
65 #endif
66
67 __inline extern double
sin(double x)68 sin (double x)
69 {
70 double value;
71
72 __asm ("fsin%.x %1,%0"
73 : "=f" (value)
74 : "f" (x));
75 return value;
76 }
77
78 __inline extern double
cos(double x)79 cos (double x)
80 {
81 double value;
82
83 __asm ("fcos%.x %1,%0"
84 : "=f" (value)
85 : "f" (x));
86 return value;
87 }
88
89 __inline extern double
tan(double x)90 tan (double x)
91 {
92 double value;
93
94 __asm ("ftan%.x %1,%0"
95 : "=f" (value)
96 : "f" (x));
97 return value;
98 }
99
100 __inline extern double
asin(double x)101 asin (double x)
102 {
103 double value;
104
105 __asm ("fasin%.x %1,%0"
106 : "=f" (value)
107 : "f" (x));
108 return value;
109 }
110
111 __inline extern double
acos(double x)112 acos (double x)
113 {
114 double value;
115
116 __asm ("facos%.x %1,%0"
117 : "=f" (value)
118 : "f" (x));
119 return value;
120 }
121
122 __inline extern double
atan(double x)123 atan (double x)
124 {
125 double value;
126
127 __asm ("fatan%.x %1,%0"
128 : "=f" (value)
129 : "f" (x));
130 return value;
131 }
132
133 __inline extern double
atan2(double y,double x)134 atan2 (double y, double x)
135 {
136 double pi, pi_over_2;
137
138 __asm ("fmovecr%.x #0,%0" /* extended precision pi */
139 : "=f" (pi)
140 : /* no inputs */ );
141 __asm ("fscale%.b #-1,%0" /* no loss of accuracy */
142 : "=f" (pi_over_2)
143 : "0" (pi));
144 if (x > 0)
145 {
146 if (y > 0)
147 {
148 if (x > y)
149 return atan (y / x);
150 else
151 return pi_over_2 - atan (x / y);
152 }
153 else
154 {
155 if (x > -y)
156 return atan (y / x);
157 else
158 return - pi_over_2 - atan (x / y);
159 }
160 }
161 else
162 {
163 if (y < 0)
164 {
165 if (-x > -y)
166 return - pi + atan (y / x);
167 else
168 return - pi_over_2 - atan (x / y);
169 }
170 else
171 {
172 if (-x > y)
173 return pi + atan (y / x);
174 else if (y > 0)
175 return pi_over_2 - atan (x / y);
176 else
177 {
178 double value;
179
180 errno = EDOM;
181 __asm ("fmove%.d #0x7fffffffffffffff,%0" /* quiet NaN */
182 : "=f" (value)
183 : /* no inputs */);
184 return value;
185 }
186 }
187 }
188 }
189
190 __inline extern double
sinh(double x)191 sinh (double x)
192 {
193 double value;
194
195 __asm ("fsinh%.x %1,%0"
196 : "=f" (value)
197 : "f" (x));
198 return value;
199 }
200
201 __inline extern double
cosh(double x)202 cosh (double x)
203 {
204 double value;
205
206 __asm ("fcosh%.x %1,%0"
207 : "=f" (value)
208 : "f" (x));
209 return value;
210 }
211
212 __inline extern double
tanh(double x)213 tanh (double x)
214 {
215 double value;
216
217 __asm ("ftanh%.x %1,%0"
218 : "=f" (value)
219 : "f" (x));
220 return value;
221 }
222
223 __inline extern double
atanh(double x)224 atanh (double x)
225 {
226 double value;
227
228 __asm ("fatanh%.x %1,%0"
229 : "=f" (value)
230 : "f" (x));
231 return value;
232 }
233
234 __inline extern double
exp(double x)235 exp (double x)
236 {
237 double value;
238
239 __asm ("fetox%.x %1,%0"
240 : "=f" (value)
241 : "f" (x));
242 return value;
243 }
244
245 __inline extern double
expm1(double x)246 expm1 (double x)
247 {
248 double value;
249
250 __asm ("fetoxm1%.x %1,%0"
251 : "=f" (value)
252 : "f" (x));
253 return value;
254 }
255
256 __inline extern double
log(double x)257 log (double x)
258 {
259 double value;
260
261 __asm ("flogn%.x %1,%0"
262 : "=f" (value)
263 : "f" (x));
264 return value;
265 }
266
267 __inline extern double
log1p(double x)268 log1p (double x)
269 {
270 double value;
271
272 __asm ("flognp1%.x %1,%0"
273 : "=f" (value)
274 : "f" (x));
275 return value;
276 }
277
278 __inline extern double
log10(double x)279 log10 (double x)
280 {
281 double value;
282
283 __asm ("flog10%.x %1,%0"
284 : "=f" (value)
285 : "f" (x));
286 return value;
287 }
288
289 __inline extern double
sqrt(double x)290 sqrt (double x)
291 {
292 double value;
293
294 __asm ("fsqrt%.x %1,%0"
295 : "=f" (value)
296 : "f" (x));
297 return value;
298 }
299
300 __inline extern double
hypot(double x,double y)301 hypot (double x, double y)
302 {
303 return sqrt (x*x + y*y);
304 }
305
306 __inline extern double
pow(double x,double y)307 pow (double x, double y)
308 {
309 if (x > 0)
310 return exp (y * log (x));
311 else if (x == 0)
312 {
313 if (y > 0)
314 return 0.0;
315 else
316 {
317 double value;
318
319 errno = EDOM;
320 __asm ("fmove%.d #0x7fffffffffffffff,%0" /* quiet NaN */
321 : "=f" (value)
322 : /* no inputs */);
323 return value;
324 }
325 }
326 else
327 {
328 double temp;
329
330 __asm ("fintrz%.x %1,%0"
331 : "=f" (temp) /* integer-valued float */
332 : "f" (y));
333 if (y == temp)
334 {
335 int i = (int) y;
336
337 if ((i & 1) == 0) /* even */
338 return exp (y * log (-x));
339 else
340 return - exp (y * log (-x));
341 }
342 else
343 {
344 double value;
345
346 errno = EDOM;
347 __asm ("fmove%.d #0x7fffffffffffffff,%0" /* quiet NaN */
348 : "=f" (value)
349 : /* no inputs */);
350 return value;
351 }
352 }
353 }
354
355 __inline extern double
fabs(double x)356 fabs (double x)
357 {
358 double value;
359
360 __asm ("fabs%.x %1,%0"
361 : "=f" (value)
362 : "f" (x));
363 return value;
364 }
365
366 __inline extern double
ceil(double x)367 ceil (double x)
368 {
369 int rounding_mode, round_up;
370 double value;
371
372 __asm volatile ("fmove%.l %!,%0"
373 : "=dm" (rounding_mode)
374 : /* no inputs */ );
375 round_up = rounding_mode | 0x30;
376 __asm volatile ("fmove%.l %0,%!"
377 : /* no outputs */
378 : "dmi" (round_up));
379 __asm volatile ("fint%.x %1,%0"
380 : "=f" (value)
381 : "f" (x));
382 __asm volatile ("fmove%.l %0,%!"
383 : /* no outputs */
384 : "dmi" (rounding_mode));
385 return value;
386 }
387
388 __inline extern double
floor(double x)389 floor (double x)
390 {
391 int rounding_mode, round_down;
392 double value;
393
394 __asm volatile ("fmove%.l %!,%0"
395 : "=dm" (rounding_mode)
396 : /* no inputs */ );
397 round_down = (rounding_mode & ~0x10)
398 | 0x20;
399 __asm volatile ("fmove%.l %0,%!"
400 : /* no outputs */
401 : "dmi" (round_down));
402 __asm volatile ("fint%.x %1,%0"
403 : "=f" (value)
404 : "f" (x));
405 __asm volatile ("fmove%.l %0,%!"
406 : /* no outputs */
407 : "dmi" (rounding_mode));
408 return value;
409 }
410
411 __inline extern double
rint(double x)412 rint (double x)
413 {
414 int rounding_mode, round_nearest;
415 double value;
416
417 __asm volatile ("fmove%.l %!,%0"
418 : "=dm" (rounding_mode)
419 : /* no inputs */ );
420 round_nearest = rounding_mode & ~0x30;
421 __asm volatile ("fmove%.l %0,%!"
422 : /* no outputs */
423 : "dmi" (round_nearest));
424 __asm volatile ("fint%.x %1,%0"
425 : "=f" (value)
426 : "f" (x));
427 __asm volatile ("fmove%.l %0,%!"
428 : /* no outputs */
429 : "dmi" (rounding_mode));
430 return value;
431 }
432
433 __inline extern double
fmod(double x,double y)434 fmod (double x, double y)
435 {
436 double value;
437
438 __asm ("fmod%.x %2,%0"
439 : "=f" (value)
440 : "0" (x),
441 "f" (y));
442 return value;
443 }
444
445 __inline extern double
drem(double x,double y)446 drem (double x, double y)
447 {
448 double value;
449
450 __asm ("frem%.x %2,%0"
451 : "=f" (value)
452 : "0" (x),
453 "f" (y));
454 return value;
455 }
456
457 __inline extern double
scalb(double x,int n)458 scalb (double x, int n)
459 {
460 double value;
461
462 __asm ("fscale%.l %2,%0"
463 : "=f" (value)
464 : "0" (x),
465 "dmi" (n));
466 return value;
467 }
468
469 __inline extern double
logb(double x)470 logb (double x)
471 {
472 double exponent;
473
474 __asm ("fgetexp%.x %1,%0"
475 : "=f" (exponent)
476 : "f" (x));
477 return exponent;
478 }
479
480 __inline extern double
ldexp(double x,int n)481 ldexp (double x, int n)
482 {
483 double value;
484
485 __asm ("fscale%.l %2,%0"
486 : "=f" (value)
487 : "0" (x),
488 "dmi" (n));
489 return value;
490 }
491
492 __inline extern double
frexp(double x,int * exp)493 frexp (double x, int *exp)
494 {
495 double float_exponent;
496 int int_exponent;
497 double mantissa;
498
499 __asm ("fgetexp%.x %1,%0"
500 : "=f" (float_exponent) /* integer-valued float */
501 : "f" (x));
502 int_exponent = (int) float_exponent;
503 __asm ("fgetman%.x %1,%0"
504 : "=f" (mantissa) /* 1.0 <= mantissa < 2.0 */
505 : "f" (x));
506 if (mantissa != 0)
507 {
508 __asm ("fscale%.b #-1,%0"
509 : "=f" (mantissa) /* mantissa /= 2.0 */
510 : "0" (mantissa));
511 int_exponent += 1;
512 }
513 *exp = int_exponent;
514 return mantissa;
515 }
516
517 __inline extern double
modf(double x,double * ip)518 modf (double x, double *ip)
519 {
520 double temp;
521
522 __asm ("fintrz%.x %1,%0"
523 : "=f" (temp) /* integer-valued float */
524 : "f" (x));
525 *ip = temp;
526 return x - temp;
527 }
528
529 #endif /* not __math_68881 */
530