<p>对于real<em>x</em>,您可以使用log的系列扩展(Ei(<em>x</em>)来表示<em>x</em>→ ∞ (参见<a href="https://en.wikipedia.org/wiki/Exponential_integral#Asymptotic_(divergent)_series" rel="nofollow noreferrer">here</a>),这对于大型<code>x</code>非常准确</p>
<p>通过一些快速实验,当x>;=50(及
这就是<code>scipy</code>开始失去全部精度的地方)。而且系列的扩展也很不错,
系数是阶乘数,因此我们可以使用精确的计算,而不会产生灾难性的后果
取消:</p>
<pre><code>def _exp1ln_taylor(x, k):
r = 0; s = -1
for i in range(k, 0, -1):
r = (r + s) * i / x
s = -s
return r*s
def exp1ln(x):
x = np.array(x)
use_approx = x >= 50
# Avoid infinity/underflows from outside domain.
ax = np.where(use_approx, x, 100)
sx = np.where(use_approx, 1, x)
approx = -x + -np.log(x) + np.log1p(_exp1ln_taylor(ax, 18))
sp = np.log(scipy.special.exp1(sx))
return np.where(use_approx, approx, sp) * 1
</code></pre>
<p>与WolframAlpha(可提供数百位数字)的比较:</p>
<pre><code>x exp1ln(x) WolframAlpha
0.1 0.6004417824948862 0.6004417824948862
1 -1.5169319590020440 -1.5169319590020456
10 -12.390724371937408 -12.390724371937408
100 -104.61502435050535 -104.61502435050535
1000 -1006.9087537832978 -1006.9087537832978
1000000000 -1000000020.7232659 -1000000020.7232658
</code></pre>