Machine epsilon explained

\varepsilon .

There are two prevailing definitions, denoted here as rounding machine epsilon and interval machine epsilon. The machine epsilon is dependent on the type of rounding used and is also called unit roundoff, which has the symbol bold Roman u. However, in the commonly used definition, machine epsilon is independent of rounding method, and may be equivalent to u or 2u.

Values for standard hardware arithmetics

The following table lists machine epsilon values for standard floating-point formats.

IEEE 754 - 2008	Common name	C++ data type	Base b	Precision p	Rounding machine epsilon b ^{-(p-1)} / 2	Interval machine epsilon b ^{-(p-1)}
binary16	half precision		2	11 (one bit is implicit)	2^−11 ≈ 4.88e-04	2^−10 ≈ 9.77e-04
binary32	single precision		2	24 (one bit is implicit)	2^−24 ≈ 5.96e-08	2^−23 ≈ 1.19e-07
binary64	double precision		2	53 (one bit is implicit)	2^−53 ≈ 1.11e-16	2^−52 ≈ 2.22e-16
	extended precision, long double	^[1]	2	64	2^−64 ≈ 5.42e-20	2^−63 ≈ 1.08e-19
binary128	quad(ruple) precision		2	113 (one bit is implicit)	2^−113 ≈ 9.63e-35	2^−112 ≈ 1.93e-34
decimal32	single precision decimal	^[2]	10	7	5 × 10^−7	10^−6
decimal64	double precision decimal		10	16	5 × 10^−16	10^−15
decimal128	quad(ruple) precision decimal		10	34	5 × 10^−34	10^−33

Formal definition (Rounding machine epsilon)

Rounding is a procedure for choosing the representation of a real number in a floating point number system. For a number system and a rounding procedure, machine epsilon is the maximum relative error of the chosen rounding procedure.

Some background is needed to determine a value from this definition. A floating point number system is characterized by a radix which is also called the base,

b , and by the precision p , i.e. the number of radix b digits of the significand (including any leading implicit bit). All the numbers with the same exponent, e , have the spacing, b ^{e-(p-1)} . The spacing changes at the numbers that are perfect powers of b ; the spacing on the side of larger magnitude is b times larger than the spacing on the side of smaller magnitude.

Since machine epsilon is a bound for relative error, it suffices to consider numbers with exponent

e =0 . It also suffices to consider positive numbers. For the usual round-to-nearest kind of rounding, the absolute rounding error is at most half the spacing, or b ^{-(p-1)} / 2 . This value is the biggest possible numerator for the relative error. The denominator in the relative error is the number being rounded, which should be as small as possible to make the relative error large. The worst relative error therefore happens when rounding is applied to numbers of the form 1+ a where a is between 0 and b ^{-(p-1)} / 2 . All these numbers round to 1 with relative error a/( 1+ a) . The maximum occurs when a is at the upper end of its range. The 1+ a in the denominator is negligible compared to the numerator, so it is left off for expediency, and just b ^{-(p-1)} / 2 is taken as machine epsilon. As has been shown here, the relative error is worst for numbers that round to 1 , so machine epsilon also is called unit roundoff meaning roughly "the maximum error that can occur when rounding to the unit value".

Thus, the maximum spacing between a normalised floating point number,

x , and an adjacent normalised number is 2 \varepsilon |x| .^[3]

Arithmetic model

Numerical analysis uses machine epsilon to study the effects of rounding error. The actual errors of machine arithmetic are far too complicated to be studied directly, so instead, the following simple model is used. The IEEE arithmetic standard says all floating-point operations are done as if it were possible to perform the infinite-precision operation, and then, the result is rounded to a floating-point number. Suppose (1)

x , y are floating-point numbers, (2) \bullet is an arithmetic operation on floating-point numbers such as addition or multiplication, and (3) \circ is the infinite precision operation. According to the standard, the computer calculates:x \bullet y = round (x \circ y)

By the meaning of machine epsilon, the relative error of the rounding is at most machine epsilon in magnitude, so:

x \bullet y = (x \circ y)( 1 + z)

where

z in absolute magnitude is at most \varepsilon or u. The books by Demmel and Higham in the references can be consulted to see how this model is used to analyze the errors of, say, Gaussian elimination.

Mainstream definition (Interval machine epsilon)

The IEEE standard does not define the terms machine epsilon and unit roundoff, so differing definitions of these terms are in use, which can cause some confusion.

The formal definition for machine epsilon is the one used by Prof. James Demmel in lecture scripts,^[4] the LAPACK linear algebra package,^[5] numerics research papers^[6] and some scientific computing software.^[7] Most numerical analysts use the words machine epsilon and unit roundoff interchangeably with this meaning.

This alternative definition is significantly more widespread: machine epsilon is the difference between 1 and the next larger floating point number.

By this definition, ε equals the value of the unit in the last place relative to 1, i.e.

b ^{-(p-1)} (where is the base of the floating point system and is the precision) and the unit roundoff is u = ε / 2, assuming round-to-nearest mode, and u = ε, assuming round-by-chop.

The prevalence of this definition is rooted in its use in the ISO C Standard for constants relating to floating-point types^{[8] [9]} and corresponding constants in other programming languages.^{[10] [11] [12]} It is also widely used in scientific computing software^{[13] [14] [15]} and in the numerics and computing literature.^{[16] [17] [18] [19]}

How to determine machine epsilon

Where standard libraries do not provide precomputed values (as does with FLT_EPSILON, DBL_EPSILON and LDBL_EPSILON for C and does with std::numeric_limits<<var>T</var>>::epsilon in C++), the best way to determine machine epsilon is to refer to the table, above, and use the appropriate power formula. Computing machine epsilon is often given as a textbook exercise. The following examples compute interval machine epsilon in the sense of the spacing of the floating point numbers at 1 rather than in the sense of the unit roundoff.

Note that results depend on the particular floating-point format used, such as float, double, long double, or similar as supported by the programming language, the compiler, and the runtime library for the actual platform.

Some formats supported by the processor might not be supported by the chosen compiler and operating system. Other formats might be emulated by the runtime library, including arbitrary-precision arithmetic available in some languages and libraries.

In a strict sense the term machine epsilon means the

1 + \varepsilon accuracy directly supported by the processor (or coprocessor), not some 1 + \varepsilon accuracy supported by a specific compiler for a specific operating system, unless it's known to use the best format.

IEEE 754 floating-point formats have the property that, when reinterpreted as a two's complement integer of the same width, they monotonically increase over positive values and monotonically decrease over negative values (see the binary representation of 32 bit floats). They also have the property that

0 < |f(x)| < inf ty , and |f(x +1 ) - f(x)| \ge |f(x) - f(x -1 )| (where f(x) is the aforementioned integer reinterpretation of x ). In languages that allow type punning and always use IEEE 754–1985, we can exploit this to compute a machine epsilon in constant time. For example, in C:

typedef union dbl_64;

double machine_eps (double value)This will give a result of the same sign as value. If a positive result is always desired, the return statement of machine_eps can be replaced with: return (s.i64 < 0 ? value - s.d64 : s.d64 - value);

Example in Python:def machineEpsilon(func=float): machine_epsilon = func(1) while func(1) + machine_epsilon != func(1): machine_epsilon_last = machine_epsilon machine_epsilon = func(machine_epsilon) / func(2) return machine_epsilon_last

64-bit doubles give 2.220446e-16, which is 2⁻⁵² as expected.

Approximation

The following simple algorithm can be used to approximate the machine epsilon, to within a factor of two (one order of magnitude) of its true value, using a linear search.

epsilon = 1.0; while (1.0 + 0.5 * epsilon) ≠ 1.0: epsilon = 0.5 * epsilon

The machine epsilon, $\backslash varepsilon\_\backslash text$ can also simply be calculated as two to the negative power of the number of bits used for the mantissa.

$$\backslash varepsilon\_\backslash text\backslash \; =\backslash \; 2^$$

Relationship to absolute relative error

If $y$ is the machine representation of a number $x$ then the absolute relative error in the representation is $\backslash left|\; \backslash dfrac\; \backslash right|\; \backslash leq\; \backslash varepsilon\_\backslash text.$^[20]

Proof

The following proof is limited to positive numbers and machine representations using round-by-chop.

If $x$ is a positive number we want to represent, it will be between a machine number $x\_$ below $x$ and a machine number $x\_$ above $x$.

If $x\_\; =\; \backslash left(1.b\_b\_\; \backslash ldots\; b\_\; \backslash right)\_\; \backslash times\; 2^$, where $m$ is the number of bits used for the magnitude of the significand, then:

Since the representation of $x$ will be either $x\_b$ or $x\_u$,

Although this proof is limited to positive numbers and round-by-chop, the same method can be used to prove the inequality in relation to negative numbers and round-to-nearest machine representations.

See also

• Floating point, general discussion of accuracy issues in floating point arithmetic
• Unit in the last place (ULP)

Notes and references

• Anderson, E.; LAPACK Users' Guide, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, third edition, 1999.
• Cody, William J.; MACHAR: A Soubroutine to Dynamically Determine Machine Parameters, ACM Transactions on Mathematical Software, Vol. 14(4), 1988, 303–311.
• Besset, Didier H.; Object-Oriented Implementation of Numerical Methods, Morgan & Kaufmann, San Francisco, CA, 2000.
• Demmel, James W., Applied Numerical Linear Algebra, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.
• Higham, Nicholas J.; Accuracy and Stability of Numerical Algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2002.
• Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; and Flannery, Brian P.; Numerical Recipes in Fortran 77, 2nd ed., Chap. 20.2, pp. 881–886
• Forsythe, George E.; Malcolm, Michael A.; Moler, Cleve B.; "Computer Methods for Mathematical Computations", Prentice-Hall,, 1977

External links

• MACHAR, a routine (in C and Fortran) to "dynamically compute machine constants" (ACM algorithm 722)
• Diagnosing floating point calculations precision, Implementation of MACHAR in Component Pascal and Oberon based on the Fortran 77 version of MACHAR published in Numerical Recipes (Press et al., 1992).

Notes and References

1. https://gcc.gnu.org/onlinedocs/gcc/Floating-Types.html#Floating-Types Floating Types - Using the GNU Compiler Collection (GCC)
2. https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html#Decimal-Float Decimal Float - Using the GNU Compiler Collection (GCC)
3. Web site: Basic Issues in Floating Point Arithmetic and Error Analysis . 21 October 1999 . University of California, Berkeley . 11 June 2022 . The distance between 1 and the next larger floating point number is 2*macheps..
4. Web site: Basic Issues in Floating Point Arithmetic and Error Analysis . 21 Oct 1999 . 11 Apr 2013.
5. Web site: LAPACK Users' Guide Third Edition . 22 August 1999 . 9 March 2012.
6. Web site: David Goldberg: What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys, Vol 23, No 1, March 1991 . 11 Apr 2013.
7. Web site: Scilab documentation - number_properties - determine floating-point parameters . 11 Apr 2013.
8. Book: Jones, Derek M. . 2009 . The New C Standard - An Economic and Cultural Commentary . 377.
9. Web site: float.h reference at cplusplus.com . 11 Apr 2013.
10. Web site: std::numeric_limits reference at cplusplus.com . 11 Apr 2013.
11. Web site: Python documentation - System-specific parameters and functions . 11 Apr 2013.
12. Extended Pascal ISO 10206:1990. The value of epsreal shall be the result of subtracting 1.0 from the smallest value of real-type that is greater than 1.0..
13. Web site: Mathematica documentation: $MachineEpsilon . 11 Apr 2013.
14. Web site: Matlab documentation - eps - Floating-point relative accuracy . 11 Apr 2013 . dead . https://web.archive.org/web/20130807071540/http://www.mathworks.de/de/help/matlab/ref/eps.html. 2013-08-07.
15. Web site: Octave documentation - eps function . 11 Apr 2013.
16. Book: Higham, Nicholas . Accuracy and Stability of Numerical Algorithms (2 ed) . SIAM . 2002 . 27–28.
17. Book: Numerical Mathematics . Alfio Quarteroni. Quarteroni . Alfio . Sacco . Riccardo . Saleri . Fausto . Springer . 2000 . 0-387-98959-5 . 49 . 2013-04-11 . https://web.archive.org/web/20171114040621/http://www.techmat.vgtu.lt/~inga/Files/Quarteroni-SkaitMetod.pdf . 2017-11-14 . dead.
18. Book: Numerical Recipes . 890 . Press . William H. . Teukolsky. Saul A. . Vetterling . William T. . Flannery . Brian P..
19. Book: Numerik-Algorithmen . Gisela . Engeln-Müllges . Fritz . Reutter . 1996 . 3-18-401539-4 . 6.
20. Web site: Machine Epsilon Value for IEEE Double Precision Standard Alternative Proof Using Relative Error . 12 October 2020 . 5 May 2022.