谐波平均计算与浮点精度

Question

我在PHP中实现了Pythagorean方法,算术和几何手段是小菜一碟,但我很难想出一个可靠的调和平均实现.

这是WolframAlpha的定义:

---

这是PHP中的等效实现:

function harmonicMeanV1()
{
    $result = 0;
    $arguments = func_get_args();

    foreach ($arguments as $argument)
    {
        $result += 1 / $argument;
    }

    return func_num_args() / $result;
}

现在,如果任何参数是0这将抛出除0的警告,但由于1 / n与n ^-1相同并且pow(0, -1)优雅地返回INF常量而不抛出任何错误,我可以将其重写为以下(它仍将抛出错误,如果没有参数,但让我们暂时忽略它:

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

    return count($arguments) / array_sum($arguments);
}

两种实现在大多数情况下都能正常工作(例如v1,v2和WolframAlpha),但是如果 1/n _i 系列的总和为0,它们会失败,我应该得到另一个除0的警告,但我不...

考虑下面的一组:-2, 3, 6(WolframAlpha的说,这是一个复杂的无穷大):

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

= 0

但是,我的两个实现都返回-2.7755575615629E-17为总和(v1,v2)而不是0.

而在键盘返回的结果是-108086391056890000我的dev的机器(32位)说,这是-1.0808639105689E+17,它仍然是不一样的0或INF我期待.我甚至尝试调用is_infinite()返回值,但它false按预期返回.

我还发现了PECL扩展stats_harmonic_mean()的一部分stats功能,但令我惊讶的是我得到了完全相同的错误结果:-1.0808639105689E+17如果有任何参数0,0则返回但是没有对系列的总和进行检查,正如你所看到的在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    /* }}} */

这看起来像一个典型的浮动精度错误,但我不能真正理解为什么因为个别计算非常精确:

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

是否可以解决此问题而无需还原到gmp/ bcmathextensions？

Answer 1

你是对的。您找到的数字是浮点运算特性的产物。

增加更多的精度对你没有帮助。你所做的只是移动球门柱。

底线是计算是以有限的精度完成的。这意味着在某个时刻，中间结果将被四舍五入。那么中间结果就不再准确了。错误会通过计算传播，并最终进入最终结果。当精确结果为零时，您通常会得到大约 1e-16 的双精度数字结果。

每当您的计算涉及分母不是 2 的幂的分数时，就会发生这种情况。

解决这个问题的唯一方法是用整数或有理数（如果可以的话）来表达计算，并使用任意精度的整数包来进行计算。这就是 Wolfram|Alpha 所做的事情。

请注意，几何平均值的计算也不是微不足道的。尝试 20 次 1e20 的序列。由于数字都相同，因此结果应该是 1e20。但你会发现结果是无限的。原因是这 20 个数字的乘积 (10e400) 超出了双精度浮点数的范围，因此它被设置为无穷大。无穷大的 20 次方根仍然是无穷大。

最后，一个元观察：毕达哥拉斯的意思实际上只对正数有意义。3 和-3 的几何平均值是多少？是虚幻的吗？？您链接到的维基百科页面上的不等式链仅在所有值均为正数时才有效。