用于计算谐波系列的Python程序

Question

有谁知道如何用Python编写程序来计算谐波系列的加法.即1 + 1/2 +1/3 +1/4 ......

Answer 1

@Kiv的答案是正确的,但如果你不需要无限精度,它对于大n来说很慢.在这种情况下,最好使用渐近公式:

#!/usr/bin/env python
from math import log

def H(n):
    """Returns an approximate value of n-th harmonic number.

       http://en.wikipedia.org/wiki/Harmonic_number
    """
    # Euler-Mascheroni constant
    gamma = 0.57721566490153286060651209008240243104215933593992
    return gamma + log(n) + 0.5/n - 1./(12*n**2) + 1./(120*n**4)

@Kiv对Python 2.6 的回答:

from fractions import Fraction

harmonic_number = lambda n: sum(Fraction(1, d) for d in xrange(1, n+1))

例:

>>> N = 100
>>> h_exact = harmonic_number(N)
>>> h = H(N)
>>> rel_err = (abs(h - h_exact) / h_exact)
>>> print n, "%r" % h, "%.2g" % rel_err
100 5.1873775176396242 6.8e-16

在N = 100相对误差小于1e-15.

Answer 2

@ recursive的解决方案对于浮点近似是正确的.如果您愿意,可以使用fractions模块在Python 3.0中获得确切的答案:

>>> from fractions import Fraction
>>> def calc_harmonic(n):
...   return sum(Fraction(1, d) for d in range(1, n + 1))
...
>>> calc_harmonic(20) # sum of the first 20 terms
Fraction(55835135, 15519504)

请注意,位数会快速增长,因此对于大n来说,这将需要大量内存.如果你想得到真正的幻想,你也可以使用一个生成器来查看一系列的部分和.

Answer 3

关于使用浮点的其他答案的脚注; 首先是最大除数和迭代向下(朝向与最大值的倒数)将推迟积累的舍入误差尽可能.

Answer 4

一个快速的H功能的准确，流畅，复值版本可以使用双伽玛函数的解释来计算这里。numpy和scipy库中分别提供了Euler-Mascheroni（γ）常数和digamma函数。

from numpy import euler_gamma
from scipy.special import digamma

def digamma_H(s):
    """ If s is complex the result becomes complex. """
    return digamma(s + 1) + euler_gamma

from fractions import Fraction

def Kiv_H(n):
    return sum(Fraction(1, d) for d in xrange(1, n + 1))

def J_F_Sebastian_H(n):
    return euler_gamma + log(n) + 0.5/n - 1./(12*n**2) + 1./(120*n**4)

这是三种速度和精度方法的比较（以Kiv_H为参考）：

Kiv_H(x) J_F_Sebastian_H(x) digamma_H(x) x seconds bits seconds bits seconds bits 1 5.06e-05 exact 2.47e-06 8.8 1.16e-05 exact 10 4.45e-04 exact 3.25e-06 29.5 1.17e-05 52.6 100 7.64e-03 exact 3.65e-06 50.4 1.17e-05 exact 1000 7.62e-01 exact 5.92e-06 52.9 1.19e-05 exact

Answer 5

这应该可以解决问题。

def calc_harmonic(n):
    return sum(1.0/d for d in range(2,n+1))