blob: ba6f75a288d6e4409d2945f646c2431b248e1020 (
plain) (
blame)
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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <mach.h>
# FP_EQUALD -- The following procedure is used to compare two double precision
# numbers for equality to within the machine precision for doubles. A simple
# comparison of the difference of the two numbers with the machine epsilon
# does not suffice unless the numbers are first normalized to near 1.0, the
# constant used to compute the machine epsilon (epsilon is the smallest number
# such that 1.0 + epsilon > 1.0).
bool procedure fp_equald (x, y)
double x, y
double x1, x2, normx, normy, tol
int ex, ey
begin
# Check for the obvious first.
if (x == y)
return (true)
# We can't normalize zero, so handle the zero operand cases first.
# Note that the case 0 equals 0 is handled above.
if (x == 0.0D0 || y == 0.0D0)
return (false)
# Normalize operands and do an epsilon compare.
call fp_normd (x, normx, ex)
call fp_normd (y, normy, ey)
if (ex != ey)
return (false)
else {
tol = EPSILOND * 32.0D0
x1 = 1.0D0 + abs (normx - normy)
x2 = 1.0D0 + tol
return (x1 <= x2)
}
end
|
