/src/igraph/src/math/utils.c

Source (jump to first uncovered line)
/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2007-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program 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 General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301 USA

*/

#include "igraph_complex.h"
#include "igraph_nongraph.h"
#include "igraph_types.h"

#include <math.h>
#include <float.h>

/**
 * \function igraph_almost_equals
 * \brief Compare two double-precision floats with a tolerance.
 *
 * Determines whether two double-precision floats are "almost equal"
 * to each other with a given level of tolerance on the relative error.
 *
 * \param  a  The first float.
 * \param  b  The second float.
 * \param  eps  The level of tolerance on the relative error. The relative
 *         error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
 *         two numbers are considered equal if this is less than \c eps.
 *
 * \return True if the two floats are nearly equal to each other within
 *         the given level of tolerance, false otherwise.
 */
igraph_bool_t igraph_almost_equals(double a, double b, double eps) {
    return igraph_cmp_epsilon(a, b, eps) == 0;
}

/* Use value-safe floating point math for igraph_cmp_epsilon() with
 * the Intel compiler.
 *
 * The Intel compiler rewrites arithmetic expressions for faster
 * evaluation by default. In the below function, it will evaluate
 * (eps * fabs(a) + eps * fabs(b)) as eps*(fabs(a) + fabs(b)).
 * However, this code path is taken precisely when fabs(a) + fabs(b)
 * overflows, thus this rearrangement of the expression causes
 * the function to return incorrect results, and some test failures.
 * To avoid this, we switch the Intel compiler to "precise" mode.
 */
#ifdef __INTEL_COMPILER
#pragma float_control(push)
#pragma float_control (precise, on)
#endif

/**
 * \function igraph_cmp_epsilon
 * \brief Compare two double-precision floats with a tolerance.
 *
 * Determines whether two double-precision floats are "almost equal"
 * to each other with a given level of tolerance on the relative error.
 *
 * </para><para>
 * The function supports infinities and NaN values. NaN values are considered
 * not equal to any other value (even another NaN), but the ordering is
 * arbitrary; in other words, we only guarantee that comparing a NaN with
 * any other value will not return zero. Positive infinity is considered to
 * be greater than any finite value with any tolerance. Negative infinity is
 * considered to be smaller than any finite value with any tolerance.
 * Positive infinity is considered to be equal to another positive infinity
 * with any tolerance. Negative infinity is considered to be equal to another
 * negative infinity with any tolerance.
 *
 * \param  a  The first float.
 * \param  b  The second float.
 * \param  eps  The level of tolerance on the relative error. The relative
 *         error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
 *         two numbers are considered equal if this is less than \c eps.
 *         Negative epsilon values are not allowed; the returned value will
 *         be undefined in this case. Zero means to do an exact comparison
 *         without tolerance.
 *
 * \return Zero if the two floats are nearly equal to each other within
 *         the given level of tolerance, positive number if the first float is
 *         larger, negative number if the second float is larger.
 */
int igraph_cmp_epsilon(double a, double b, double eps) {
    double diff;
    double abs_diff;
    double sum;

    if (a == b) {
        /* shortcut, handles infinities */
        return 0;
    }

    diff = a - b;
    abs_diff = fabs(diff);
    sum = fabs(a) + fabs(b);

    if (a == 0 || b == 0 || sum < DBL_MIN) {
        /* a or b is zero or both are extremely close to it; relative
         * error is less meaningful here so just compare it with
         * epsilon */
        return abs_diff < (eps * DBL_MIN) ? 0 : (diff < 0 ? -1 : 1);
    } else if (!isfinite(sum)) {
        /* addition overflow, so presumably |a| and |b| are both large; use a
         * different formulation */
        return (abs_diff < (eps * fabs(a) + eps * fabs(b))) ? 0 : (diff < 0 ? -1 : 1);
    } else {
        return (abs_diff / sum < eps) ? 0 : (diff < 0 ? -1 : 1);
    }
}

/**
 * \function igraph_complex_almost_equals
 * \brief Compare two complex numbers with a tolerance.
 *
 * Determines whether two complex numbers are "almost equal"
 * to each other with a given level of tolerance on the relative error.
 *
 * \param  a  The first complex number.
 * \param  b  The second complex number.
 * \param  eps  The level of tolerance on the relative error. The relative
 *         error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
 *         two numbers are considered equal if this is less than \c eps.
 *
 * \return True if the two complex numbers are nearly equal to each other within
 *         the given level of tolerance, false otherwise.
 */
igraph_bool_t igraph_complex_almost_equals(igraph_complex_t a,
                                           igraph_complex_t b,
                                           igraph_real_t eps) {

    igraph_real_t a_abs = igraph_complex_abs(a);
    igraph_real_t b_abs = igraph_complex_abs(b);
    igraph_real_t sum = a_abs + b_abs;
    igraph_real_t abs_diff = igraph_complex_abs(igraph_complex_sub(a, b));

    if (a_abs == 0 || b_abs == 0 || sum < DBL_MIN) {
        return abs_diff < eps * DBL_MIN;
    } else if (! isfinite(sum)) {
        return abs_diff < (eps * a_abs + eps * b_abs);
    } else {
        return abs_diff/ sum < eps;
    }
}

#ifdef __INTEL_COMPILER
#pragma float_control(pop)
#endif

Coverage Report

Created: 2025-04-11 06:16

Line	Count	Source (jump to first uncovered line)
1		/* -- mode: C -- */
2		/*
3		IGraph library.
4		Copyright (C) 2007-2012 Gabor Csardi <csardi.gabor@gmail.com>
5		334 Harvard street, Cambridge, MA 02139 USA
6
7		This program is free software; you can redistribute it and/or modify
8		it under the terms of the GNU General Public License as published by
9		the Free Software Foundation; either version 2 of the License, or
10		(at your option) any later version.
11
12		This program is distributed in the hope that it will be useful,
13		but WITHOUT ANY WARRANTY; without even the implied warranty of
14		MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15		GNU General Public License for more details.
16
17		You should have received a copy of the GNU General Public License
18		along with this program; if not, write to the Free Software
19		Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20		02110-1301 USA
21
22		*/
23
24		#include "igraph_complex.h"
25		#include "igraph_nongraph.h"
26		#include "igraph_types.h"
27
28		#include <math.h>
29		#include <float.h>
30
31		/**
32		* \function igraph_almost_equals
33		* \brief Compare two double-precision floats with a tolerance.
34		*
35		* Determines whether two double-precision floats are "almost equal"
36		* to each other with a given level of tolerance on the relative error.
37		*
38		* \param a The first float.
39		* \param b The second float.
40		* \param eps The level of tolerance on the relative error. The relative
41		* error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
42		* two numbers are considered equal if this is less than \c eps.
43		*
44		* \return True if the two floats are nearly equal to each other within
45		* the given level of tolerance, false otherwise.
46		*/
47	0	igraph_bool_t igraph_almost_equals(double a, double b, double eps) {
48	0	return igraph_cmp_epsilon(a, b, eps) == 0;
49	0	}
50
51		/* Use value-safe floating point math for igraph_cmp_epsilon() with
52		* the Intel compiler.
53		*
54		* The Intel compiler rewrites arithmetic expressions for faster
55		* evaluation by default. In the below function, it will evaluate
56		* (eps * fabs(a) + eps * fabs(b)) as eps*(fabs(a) + fabs(b)).
57		* However, this code path is taken precisely when fabs(a) + fabs(b)
58		* overflows, thus this rearrangement of the expression causes
59		* the function to return incorrect results, and some test failures.
60		* To avoid this, we switch the Intel compiler to "precise" mode.
61		*/
62		#ifdef __INTEL_COMPILER
63		#pragma float_control(push)
64		#pragma float_control (precise, on)
65		#endif
66
67		/**
68		* \function igraph_cmp_epsilon
69		* \brief Compare two double-precision floats with a tolerance.
70		*
71		* Determines whether two double-precision floats are "almost equal"
72		* to each other with a given level of tolerance on the relative error.
73		*
74		* </para><para>
75		* The function supports infinities and NaN values. NaN values are considered
76		* not equal to any other value (even another NaN), but the ordering is
77		* arbitrary; in other words, we only guarantee that comparing a NaN with
78		* any other value will not return zero. Positive infinity is considered to
79		* be greater than any finite value with any tolerance. Negative infinity is
80		* considered to be smaller than any finite value with any tolerance.
81		* Positive infinity is considered to be equal to another positive infinity
82		* with any tolerance. Negative infinity is considered to be equal to another
83		* negative infinity with any tolerance.
84		*
85		* \param a The first float.
86		* \param b The second float.
87		* \param eps The level of tolerance on the relative error. The relative
88		* error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
89		* two numbers are considered equal if this is less than \c eps.
90		* Negative epsilon values are not allowed; the returned value will
91		* be undefined in this case. Zero means to do an exact comparison
92		* without tolerance.
93		*
94		* \return Zero if the two floats are nearly equal to each other within
95		* the given level of tolerance, positive number if the first float is
96		* larger, negative number if the second float is larger.
97		*/
98	0	int igraph_cmp_epsilon(double a, double b, double eps) {
99	0	double diff;
100	0	double abs_diff;
101	0	double sum;
102
103	0	if (a == b) {
104		/* shortcut, handles infinities */
105	0	return 0;
106	0	}
107
108	0	diff = a - b;
109	0	abs_diff = fabs(diff);
110	0	sum = fabs(a) + fabs(b);
111
112	0	if (a == 0 \|\| b == 0 \|\| sum < DBL_MIN) {
113		/* a or b is zero or both are extremely close to it; relative
114		* error is less meaningful here so just compare it with
115		* epsilon */
116	0	return abs_diff < (eps * DBL_MIN) ? 0 : (diff < 0 ? -1 : 1);
117	0	} else if (!isfinite(sum)) {
118		/* addition overflow, so presumably \|a\| and \|b\| are both large; use a
119		* different formulation */
120	0	return (abs_diff < (eps * fabs(a) + eps * fabs(b))) ? 0 : (diff < 0 ? -1 : 1);
121	0	} else {
122	0	return (abs_diff / sum < eps) ? 0 : (diff < 0 ? -1 : 1);
123	0	}
124	0	}
125
126		/**
127		* \function igraph_complex_almost_equals
128		* \brief Compare two complex numbers with a tolerance.
129		*
130		* Determines whether two complex numbers are "almost equal"
131		* to each other with a given level of tolerance on the relative error.
132		*
133		* \param a The first complex number.
134		* \param b The second complex number.
135		* \param eps The level of tolerance on the relative error. The relative
136		* error is defined as <code>abs(a-b) / (abs(a) + abs(b))</code>. The
137		* two numbers are considered equal if this is less than \c eps.
138		*
139		* \return True if the two complex numbers are nearly equal to each other within
140		* the given level of tolerance, false otherwise.
141		*/
142		igraph_bool_t igraph_complex_almost_equals(igraph_complex_t a,
143		igraph_complex_t b,
144	0	igraph_real_t eps) {
145
146	0	igraph_real_t a_abs = igraph_complex_abs(a);
147	0	igraph_real_t b_abs = igraph_complex_abs(b);
148	0	igraph_real_t sum = a_abs + b_abs;
149	0	igraph_real_t abs_diff = igraph_complex_abs(igraph_complex_sub(a, b));
150
151	0	if (a_abs == 0 \|\| b_abs == 0 \|\| sum < DBL_MIN) {
152	0	return abs_diff < eps * DBL_MIN;
153	0	} else if (! isfinite(sum)) {
154	0	return abs_diff < (eps * a_abs + eps * b_abs);
155	0	} else {
156	0	return abs_diff/ sum < eps;
157	0	}
158	0	}
159
160		#ifdef __INTEL_COMPILER
161		#pragma float_control(pop)
162		#endif