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

Source (jump to first uncovered line)
/*
   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-07-01 07:00

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