Formerly ieee754/cabs.c.~2~
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / drem.c
1 /* Copyright (C) 1992 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3
4 The GNU C Library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Library General Public License as
6 published by the Free Software Foundation; either version 2 of the
7 License, or (at your option) any later version.
8
9 The GNU C Library is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12 Library General Public License for more details.
13
14 You should have received a copy of the GNU Library General Public
15 License along with the GNU C Library; see the file COPYING.LIB.  If
16 not, write to the Free Software Foundation, Inc., 675 Mass Ave,
17 Cambridge, MA 02139, USA.  */
18
19 /*
20  * Copyright (c) 1985 Regents of the University of California.
21  * All rights reserved.
22  *
23  * Redistribution and use in source and binary forms are permitted provided
24  * that: (1) source distributions retain this entire copyright notice and
25  * comment, and (2) distributions including binaries display the following
26  * acknowledgement:  ``This product includes software developed by the
27  * University of California, Berkeley and its contributors'' in the
28  * documentation or other materials provided with the distribution and in
29  * all advertising materials mentioning features or use of this software.
30  * Neither the name of the University nor the names of its contributors may
31  * be used to endorse or promote products derived from this software without
32  * specific prior written permission.
33  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
34  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
35  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
36  */
37
38 #include <ansidecl.h>
39 #include <math.h>
40 #include <float.h>
41 #include "ieee754.h"
42
43 /* Return the remainder of X/Y.  */
44 double
45 DEFUN(__drem, (x, y),
46       double x AND double y)
47 {
48   union ieee754_double ux, uy;
49
50   ux.d = x;
51   uy.d = y;
52 #define x ux.d
53 #define y uy.d
54
55   uy.ieee.negative = 0;
56
57   if (!__finite (x) || y == 0.0)
58     return NAN;
59   else if (__isnan (y))
60     return y;
61   else if (__isinf (y))
62     return x;
63   else if (uy.ieee.exponent <= 1)
64     {
65       /* Subnormal (or almost subnormal) Y value.  */
66       double b = __scalb (1.0, 54);
67       y *= b;
68       x = __drem (x, y);
69       x *= b;
70       return __drem (x, y) / b;
71     }
72   else if (y >= 1.7e308 / 2)
73     {
74       y /= 2;
75       x /= 2;
76       return __drem (x, y) * 2;
77     }
78   else
79     {
80       union ieee754_double a;
81       double b;
82       unsigned int negative = ux.ieee.negative;
83       a.d = y + y;
84       b = y / 2;
85       ux.ieee.negative = 0;
86       while (x > a.d)
87         {
88           unsigned short int k = ux.ieee.exponent - a.ieee.exponent;
89           union ieee754_double tmp;
90           tmp.d = a.d;
91           tmp.ieee.exponent += k;
92           if (x < tmp.d)
93             --tmp.ieee.exponent;
94           x -= tmp.d;
95         }
96       if (x > b)
97         {
98           x -= y;
99           if (x >= b)
100             x -= y;
101         }
102       ux.ieee.negative ^= negative;
103       return x;
104     }
105 }