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
68sin (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
79cos (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
90tan (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
101asin (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
112acos (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
123atan (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
134atan2 (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
191sinh (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
202cosh (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
213tanh (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
224atanh (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
235exp (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
246expm1 (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
257log (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
268log1p (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
279log10 (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
290sqrt (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
301hypot (double x, double y)
302{
303  return sqrt (x*x + y*y);
304}
305
306__inline extern double
307pow (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
356fabs (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
367ceil (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
389floor (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
412rint (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
434fmod (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
446drem (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
458scalb (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
470logb (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
481ldexp (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
493frexp (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
518modf (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