From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- sys/gio/fpequald.x | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 sys/gio/fpequald.x (limited to 'sys/gio/fpequald.x') diff --git a/sys/gio/fpequald.x b/sys/gio/fpequald.x new file mode 100644 index 00000000..ba6f75a2 --- /dev/null +++ b/sys/gio/fpequald.x @@ -0,0 +1,41 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include + +# 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 -- cgit