Best way to calculate Bernoulli's nTh numbers in C with CodeBlocks on Windows 10 Pro

Question

Good night, I'm working in a project that needs to calculate the Bernoulli numbers for nTh order. I tried exhaustively a lot of algorithms on internet, mostly it's in C++ what's not usefull for me. And always I'd got compilations erros or the wrong results to the numbers! What is the fastest way to calculate it? Odd's numbers is always 0.00000 and I need to calculate for any even numbers... I just need the result of the number I put on the function, don't need to list the numbers until the nTh like every algorithm I saw on the internet does. The last I tried that had compilation erros and after fixed give me wrong answers above... Yes, for people that will ask me if I put the libraries on the code, yes, I did it! The problem is not libraries, it's wrong algorithms. I'm using C on GCC mingw 32-bits on Code::Blocks for Windows 10 Pro.

#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);
  }
}

I've never worked with this type of stuff before, but on the series I'am working requires the Bernoulli numbers. So I don't have many sure what I'm doing to find those numbers. It's not my area. Probably I made some stupid thing here.

I'll tell you how I fall in this problem of Bernoulli, I'm originally working on Riemann Zeta's function. I made the C code but it only worked for >1, So I started to study how to calculate for Negative odd's, and I saw that Bn(Bernoulli numbers of N order) are in the formulae! I don't know how to calculate Bernoulli Numbers, and when I started to code Zeta's function I didn't know nothing about Bernoulli!

Answer 1

My advice is to use a library or package to carry out this computation; it is easy to write some code for a mathematical function which handles simple cases, and very difficult and time-consuming to handle all cases correctly. Presumably the calculation of Bernoulli numbers is just something you need to make progress on your real topic of interest. If so, you're better off finding an existing library or package. (Even you are having trouble with a library, it is still far easier to solve that problem instead of having to reimplement the algorithm.)

Sage (https://sagemath.org) can calculate Bernoulli numbers, and probably has a lot of other number theory stuff. See also Maxima (http://maxima.sourceforge.net) and maybe also GAP and PARI-GP (a web search will find those).