Harmonic calculation and float accuracy

Question

I implement Pythagorean means in PHP, arithmetic and geometric means is a piece of cake, but I have a very difficult time with a reliable harmonic mean .

This is the definition of WolframAlpha :

---

And this is the equivalent implementation in PHP:

function harmonicMeanV1() { $result = 0; $arguments = func_get_args(); foreach ($arguments as $argument) { $result += 1 / $argument; } return func_num_args() / $result; } 

Now, if any of the arguments is 0 , this will cause division by 0, but since 1 / n same as n ^-1 and pow(0, -1) gracefully returns the INF constant without any errors, I could rewrite this as follows ( it will still cause errors if there are no arguments, but this time ignore):

 function harmonicMeanV2() { $arguments = func_get_args(); $arguments = array_map('pow', $arguments, array_fill(0, count($arguments), -1)); return count($arguments) / array_sum($arguments); } 

Both implementations work fine in most cases (example v1 , v2 and WolframAlpha ), but they look inefficient if the sum 1 / n _i is 0, I should get another division by 0, but I don’t ...

Consider the following set: -2, 3, 6 ( WolframAlpha says it's complex infinite):

  1 / -2 // -0.5 + 1 / 3 // 0.33333333333333333333333333333333 + 1 / 6 // 0.16666666666666666666666666666667 = 0 

However, both of my implementations return -2.7755575615629E-17 as the sum ( v1 , v2 ) instead of 0 .

While the result of returning to CodePad is -108086391056890000 , my dev machine (32-bit version) says it is -1.0808639105689E+17 , it still doesn't look like 0 or INF , which I expected. I even tried calling is_infinite() on the return value, but it returned as false as expected.

I also found the stats_harmonic_mean() function, which is part of the PECL stats extension, but, to my surprise, I got exactly the same error result: -1.0808639105689E+17 , if any of the arguments are 0 , it returns 0 , but not the series amount is checked, as you can see on line line 3585 :

 3557 /* {{{ proto float stats_harmonic_mean(array a) 3558 Returns the harmonic mean of an array of values */ 3559 PHP_FUNCTION(stats_harmonic_mean) 3560 { 3561 zval *arr; 3562 double sum = 0.0; 3563 zval **entry; 3564 HashPosition pos; 3565 int elements_num; 3566 3567 if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { 3568 return; 3569 } 3570 if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { 3571 php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); 3572 RETURN_FALSE; 3573 } 3574 3575 zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); 3576 while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { 3577 convert_to_double_ex(entry); 3578 if (Z_DVAL_PP(entry) == 0) { 3579 RETURN_LONG(0); 3580 } 3581 sum += 1 / Z_DVAL_PP(entry); 3582 zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); 3583 } 3584 3585 RETURN_DOUBLE(elements_num / sum); 3586 } 3587 /* }}} */ 

This seems like a typical floating-point error, but I can’t understand the reason, because the individual calculations are accurate enough:

 Array ( [0] => -0.5 [1] => 0.33333333333333 [2] => 0.16666666666667 ) 

Can I get around this problem without returning to the gmp / bcmath ?

Answer 1

You're right. The numbers you find are an artifact of the features of floating point arithmetic.

Adding more accuracy will not help. All you do is move target posts.

The bottom line is that the calculations are performed with finite precision. This means that at some point, the intermediate result will be rounded. Then this intermediate result is no longer accurate. The error propagates through computation and ultimately turns it into the final result. When the exact result is zero, you usually get a numerical result around 1e-16 with double precision numbers.

This happens every time your calculation includes a fraction with a denominator that is not a power of 2.

The only way is to express the calculations in terms of integers or rational numbers (if possible) and use an arbitrary array of integers to perform the calculations. This is what Wolfram Alpha does.

Note that calculating the geometric mean is also not trivial. Try a sequence of 20 times 1e20. Since the numbers are all the same, the result should be 1e20. But you will find that the result is endless. The reason is that the product of these 20 numbers (10e400) goes beyond the number of floating-point numbers with double precision and therefore is set to infinity. The 20th root of infinity is still infinite.

Finally, a meta-observation: Pythogarian really means only positive numbers. What is the geometric mean of 3 and -3? Is this imaginary? The inequality chain on the Wikipedia page that you link to is valid only if all values ​​are positive.

Answer 2

Yes, this is a floating point precision issue. -1/2 can be represented exactly, but 1/3 and 1/6 cannot. That way, when you add them, you are not getting zero.

You can use the use-a-common-denominator approach you mentioned (the H2 and H3 formulas that you published), but it just hits the box office a bit along the way, you still get inaccurate results, the final product starts to round.

Why do you accept the average of the harmonics of numbers, which may be negative? That an inherently unstable calculation (H (-2.3.6 + epsilon) varies widely for very small epsilon).