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