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