Best way to calculate machine epsilon in PHP?

Question

What is best / most correct way of directly calculating machine epsilon (floating point rounding error) in PHP, not using built-in constant PHP_FLOAT_EPSILON ? Currently, I've managed to research two ways, "standard" and asymptotically approaching epsilon :

    // Standard way
    $x=1.0;
    $dx=2.0; 
    while (($x/=$dx)+1.0!=1.0); 
    echo "EPSILON : $x".", DEFINED/CALC : ".round(PHP_FLOAT_EPSILON/$x)."x\n";

    // Asymptotically approaching
    $x=1.0;
    $dx=2.0; 
    while (($x/=$dx)+1.0!=1.0) $dx/=1.0+10**-5;
    echo "EPSILON : $x".", DEFINED/CALC : ".round(PHP_FLOAT_EPSILON/$x)."x\n";

The problem is, that they gives answers with different errors :

EPSILON : 1.1102230246252E-16, DEFINED/CALC : 2x
EPSILON : 5.6311222663283E-17, DEFINED/CALC : 4x

So standard gives ε = 1/2ε₀ and other algo gives ε = 1/4ε₀. I'm not sure how epsilon should be calculated right.

It may be better to ask this on anothe stack site, maybe https://math.stackexchange.com/ — RiggsFolly, Oct 20 '20 at 12:05
Maybe, or maybe not. It may depend also on language implementation of floating point format. In any case, it's not just math floating point format implementation is computer stuff. — Agnius Vasiliauskas, Oct 20 '20 at 12:07
@RiggsFolly: No, the issues involved here are computing issues. In particular, rounding causes `($x/=$dx)+1.0` not to equal `1.0` even when `($x/=$dx)` is less than the “machine epsilon” (but above half it). This is a computing issue, not a pure mathematical issue. — Eric Postpischil, Oct 20 '20 at 15:33
If you have seen a definition of the “machine epsilon” as the smallest value x for which `1+x` is not equal to `1`, that is incorrect. The correct definition of “machine epsilon” is that it is the difference between 1 and the next representable number. Call this u. So 1 and 1+e are representable numbers. Suppose x is some number in (½u, u). Then, when `1+x` is computed, the result is 1+u because it is rounded up—the real arithmetic result is between 1 and 1+u, which are the two nearest representable values. Since it is closer to 1+u, it is rounded up to 1+u. — Eric Postpischil, Oct 20 '20 at 15:37
So, that is why the first definition is wrong. A second is that your code shows the value of `$x` **after** the loop terminates, meaning it shows the first value of `$x` for which the two expressions were equal instead of showing the last value of `$x` for which the two expressions were unequal. — Eric Postpischil, Oct 20 '20 at 15:40
@EricPostpischil Thanks, showing one value after $x, was main cause of fault. When I print previous value which was !=1, those both cases reports `1ε` and `2ε`. It's ok for me, at least first case corresponds to an epsilon definition in PHP documentation. Thanks ! However it would be nice to see your comments as a full-scale answer with your php or maybe at least pseudo code, how to calculate epsilon by your suggested definition / approach. — Agnius Vasiliauskas, Oct 20 '20 at 17:57

Agnius Vasiliauskas · Accepted Answer · 2020-10-21T08:09:38.003

Thanks to @EricPostpischil, main cause of error was printing first value of x for which x+1=1 and instead of that I should be printing last value of x for which x+1≠1. Fixed code :

// Standard way
$x=1.0;
$dx=2.0; 
while (true) {
  $px = $x;
  $x/=$dx;
  if ($x+1.0==1.0) 
    break;
} 
printf ("ε = $x ≈ %.1f ε₀ \n", $px/PHP_FLOAT_EPSILON);

// Asymptotically approaching
$x=1.0;
$dx=2.0; 
while (true) {
  $px = $x;
  $x/=$dx;
  $dx/=1.0+10**-5;
  if ($x+1.0==1.0) 
    break;
} 
printf ("ε = $x ≈ %.1f ε₀ \n", $px/PHP_FLOAT_EPSILON);

Reports =>

ε = 1.1102230246252E-16 ≈ 1.0 ε₀
ε = 5.6311222663283E-17 ≈ 0.5 ε₀

It's good enough now, because first -standard case- reports epsilon to be the same as in PHP documentation.As about second algo, it maybe due to rounding issues and can be interpreted somewhat differently. I.e. if ±ε gets you to the next representable float number, then representation error is half epsilon, because any change greater than |0.5 ε₀| will be rounded to next representable number. It's similar like physicists counts for measurement errors. Say that you have a ruler with smallest measurement unit of 1 mm. Then adding ±0.5 mm to the current reading and rounding will get you to the next representable reading on the ruler. Thus some says that real measurement error is 0.5 mm, and some that it's 1 mm, depends on definition. Same here.

Best way to calculate machine epsilon in PHP?

1 Answers1