phpprecisionfloating-accuracyepsilon

Best way to calculate machine epsilon in PHP?


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ε0 and other algo gives ε = 1/4ε0. I'm not sure how epsilon should be calculated right.


Solution

  • 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.