c - 在 Windows 10 Pro 上使用 CodeBlocks 在 C 中计算伯努利 nTh 数的最佳方法

Question

晚安，我正在从事一个需要计算第 n 阶伯努利数的项目。我在互联网上尝试了很多算法，大部分都是 C++ 的，对我来说没有用。我总是会遇到编译错误或数字结果错误!最快的计算方法是什么？奇数的数字始终是 0.00000，我需要计算任何偶数...我只需要我放在函数中的数字的结果，不需要像我在互联网上看到的每个算法一样列出直到第 nTh 的数字做。我最后一次尝试有编译错误，修复后给了我上面错误的答案...是的，对于那些会问我是否将库放在代码中的人，是的，我做到了!问题不在于库，而在于错误的算法。我在 Windows 10 Pro 的 Code::Blocks 上的 GCC mingw 32 位上使用 C。

#include <float.h>
#include <math.h>

void bernoulli_B( int iMax, double* dB )
{
  dB[0] = 1.0;
  dB[1] = -0.5;


  for( int j = 3; j <= iMax; j += 2 )
    dB[j] = 0.0;

  const double eps = DBL_EPSILON;  
  const double TwoPi = 6.2831853071795860;
  double dCoeff = 2.0 / (TwoPi * TwoPi);
  double d2 = 2.0;


  for( int n = 1; n <= iMax/2; n++ )
  {
    double   g1 = 1.0,     
    g2 = 1.0;

    for( int j = 0; j < n; j++ )
    {
      g1 *= 4.0;    
      g2 *= 9.0;                    
    }

    double   S1 = 1.0 - 1.0/g1,         
    S2 = S1 + 1.0/g2,  S3;

    double   T1 = S1 + 1.0/(g2 + g1),   
    T2;

    long r = 4;
    double s = -1.0;      
    int nSuccess = 0;

    while( !nSuccess )
    {

      double r2 = double(r*r);
      double g3 = 1.0;

      for( int j = 0; j < n; j++ )
        g3 *= r2;


      S3 = S2 + s/g3;
      T2 = S2 + s/(g3 + g2);

      if( fabs(T2-T1) > eps*fabs(T2) ) 
      {
        g2 = g3;
        S2 = S3;
        T1 = T2;
        s = -s;
        r++;
      }
      else       
      {
        nSuccess = 1;
      }
    }
    d2 /= 4.0;  
    dB[2*n] = 2.0 * dCoeff / (1.0-d2) * T2;
    dCoeff *= -double((2*n+1)*(2*n+2)) / (TwoPi * TwoPi);
  }
}

我以前从未处理过此类问题，但我正在处理的系列需要伯努利数。所以我不太确定我正在做什么来找到这些数字。这不是我的领域。 可能我在这里做了一些愚蠢的事情。

我会告诉你我是如何陷入伯努利这个问题的，我最初是在研究黎曼泽塔函数。我编写了C代码，但它只适用于>1，所以我开始研究如何计算负奇数，我看到公式中包含Bn(N阶伯努利数)!我不知道如何计算伯努利数，当我开始编写 Zeta 函数代码时，我对伯努利一无所知!

Answer 1

我的建议是使用库或包来执行此计算；为处理简单情况的数学函数编写一些代码很容易，但正确处理所有情况则非常困难且耗时。想必伯努利数的计算正是您在真正感兴趣的主题上取得进展所需要的。如果是这样，您最好找到现有的库或包。 (即使您在使用库时遇到问题，解决该问题仍然要容易得多，而不必重新实现算法。)

Sage ( https://sagemath.org ) 可以计算伯努利数，并且可能还有很多其他数论内容。另请参阅 Maxima (http://maxima.sourceforge.net)，也许还有 GAP 和 PARI-GP(网络搜索即可找到这些)。