我想评估总和
在哪里和具有 8 位精度。
如果我愿意每次计算最多花费一秒钟的 CPU 时间,我应该如何进行?
假设我想要一个非常快速的评估方案并且我愿意花费大量的工作。我应该如何进行?
任何线索、线索或建议将不胜感激。
我想评估总和
如果我愿意每次计算最多花费一秒钟的 CPU 时间,我应该如何进行?
假设我想要一个非常快速的评估方案并且我愿意花费大量的工作。我应该如何进行?
任何线索、线索或建议将不胜感激。
你试过什么了?这个完全幼稚的实现设法在我的笔记本电脑上在 2.5 秒内计算出 7 个(可能是 7.5 个)数字:
#include <iostream>
#include <complex>
#include <cmath>
#include <iomanip>
int main ()
{
const double alpha = 1;
std::cout.precision(16);
std::complex<double> sum = 0;
for (unsigned int k=1; k<10000000; ++k)
{
sum += std::pow(std::complex<double>(1,1)/std::sqrt(2.), k)
*
std::pow(k, -alpha);
if (k % 1000000 == 0)
std::cout << k << ' ' << sum << std::endl;
}
}
结果是:
g/x> c++ -O3 x.cc
g/x> time ./a.out
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)
real 0m2.553s
user 0m2.548s
sys 0m0.001s
我确信这可以大大改善,例如通过使用以下事实只能获得 8 个不同的值,可以制成表格,即可以很容易地计算出来,等等。我毫不怀疑,通过 20 分钟的工作,我可以在一秒钟内完成 8 位数的计算。对于较小的值可能更难做,然而,收敛速度较慢。
由于该问题有一个 Fortran 标签,因此这是 Wolfgang 在 Fortran 中的解决方案:
implicit none
integer, parameter :: dp = kind(0.d0)
complex(dp), parameter :: i_ = (0, 1)
real(dp) :: alpha = 1
complex(dp) :: s = 0
integer :: k
do k = 1, 10000000
s = s + ((i_+1)/sqrt(2._dp))**k * k**(-alpha)
if (modulo(k, 1000000) == 0) print *, k, s
end do
end
结果是:
$ gfortran -O3 x.f90
$ time ./a.out
1000000 ( 0.26740049836809593 , 1.1780960379893362 )
2000000 ( 0.26740024836936921 , 1.1780966415427796 )
3000000 ( 0.26740016503629038 , 1.1780968427272547 )
4000000 ( 0.26740012336968905 , 1.1780969433194985 )
5000000 ( 0.26740009836973161 , 1.1780970036749265 )
6000000 ( 0.26740008170307927 , 1.1780970439117682 )
7000000 ( 0.26740006979833991 , 1.1780970726523770 )
8000000 ( 0.26740006086978074 , 1.1780970942078899 )
9000000 ( 0.26740005392532235 , 1.1780971109732488 )
10000000 ( 0.26740004836978204 , 1.1780971243856078 )
real 0m1.477s
user 0m1.472s
sys 0m0.000s
为了比较,在我的计算机上,C++ 代码采用:
$ g++ -O3 x.cc
$ time ./a.out
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)
real 0m2.123s
user 0m2.116s
sys 0m0.004s
我使用 gcc 4.6.3。
mpmath 有这个
>>> 导入 mpmath >>> mpmath.polylog(1, (1+1j)/mpmath.sqrt(2)) mpc(真实='0.26739999836978523',图像='1.1780972450961724') >>> 导入 numpy >>> for x in numpy.linspace(0.75, 1.0, 11): print mpmath.polylog(x, mpmath.sqrt(1j)) ... (0.13205298513153 + 1.22168160828646j) (0.147070653672424 + 1.21783198879211j) (0.161744021745461 + 1.21384547018142j) (0.17607856182228 + 1.20973163332567j) (0.190079791459314 + 1.20549967730407j) (0.2037532636373 + 1.20115843051452j) (0.217104557681381 + 1.19671636161604j) (0.230139270736483 + 1.19218159029609j) (0.242863009773464 + 1.18756189785626j) (0.25528138410235 + 1.18286473761121j) (0.267399998369785 + 1.17809724509617j) >>>
如果您绘制该系列的部分总和,您将得到以下结果:
这是真正的部分对于 10 到 100 项的部分总和。显然,这里有一个模式,虚部也遵循类似的模式。似乎您可以取 8 个连续的部分和并对它们进行平均以获得更好的近似值。或者至少您可以从这些部分和中获得渐近数值行为。
更具分析性,您可以将总和分解为 8 个不同单位根上的单独总和:
每个内部和都是一个 Hurwitz zeta 函数,其中有特殊的评估例程(如GSL中)。以这种方式绘制总和作为给
蓝线是实部,红线是虚部。这似乎表现得很好。