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!

What are you actually asking about here, the compilation errors or calculation method? One question at a time please. — meowgoesthedog, Sep 18 '18 at 22:01
The calculation method. I fixed the compilation erros, but the answers is getting wrong. I have no idea why. — Adriano rox, Sep 18 '18 at 22:13
My principal work is Riemann Zeta function, but for negative odds numbers it's need to Bernoulli numbers. So I need this to the negative real part. The positive Real part is already done. — Adriano rox, Sep 18 '18 at 22:14
"but the answers is getting wrong" isn't a very good problem statement. What answers do you get and what do you expect? — jwdonahue, Sep 18 '18 at 23:09
@jwdonahue The correct numbers of Bernoulli! Did you at least read what is the numbers? — Adriano rox, Sep 18 '18 at 23:27
@meowgoesthedog I'm having a lot of issues with this GMP library, I tried this but seems that my GMP are different. Or maybe this code was written in a old version of GMP. — Adriano rox, Sep 18 '18 at 23:29
Take the [tour] and read [Ask]. You might get better results asking about the algorithm at the [Math site](https://math.stackexchange.com/). If you want to know how your code deviates from an algorithm, tell us what algorithm you are using. For me, when I read C code, I like to know what inputs you've tried and what outputs you are getting and why they are wrong. You haven't done a very good job at explaining those things. — jwdonahue, Sep 18 '18 at 23:41
@jwdonahue OMG, please. Bernoulli numbers works with even until infinite! — Adriano rox, Sep 19 '18 at 00:15
Ok, so one or two input examples and the resulting erroneous outputs should be included in your post. I really do recommend that you take the [tour] and read [ask]. — jwdonahue, Sep 19 '18 at 00:16
I am not a mathematician, but I am a pretty damn good C coder and I've taken a lot of crappy code written by phd's in mathematics and got it working, once they properly explained the intended algorithm and provided some input/output examples. — jwdonahue, Sep 19 '18 at 00:19
@jwdonahue It's getting very wrongs results that just prove how stupid I'm! — Adriano rox, Sep 19 '18 at 00:20
LOL. Sorry, I don't think you're stupid. I think you're frustrated. Step back and review the links I provided, then if you don't like Robert's answer, [edit] your question to help us help you. — jwdonahue, Sep 19 '18 at 00:22
I'm having serious issues with GMP library on Windows. To run the @meowgoesthedog example. — Adriano rox, Sep 19 '18 at 00:43
If you just want to know the Bernoulli numbers for some other purpose, why not just look them up? [Here are the numerators](http://oeis.org/A027641/b027641.txt) and [here are the denominators](http://oeis.org/A027642/b027642.txt) — AakashM, Sep 19 '18 at 09:42

score 1 · Answer 1 · answered Sep 18 '18 at 23:42

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).

answered Sep 18 '18 at 23:42

Robert Dodier

16,905
2
31
48

While not technically the answer to the OP's question, it's strong advice. I am not sure, but I think the OP really wants to figure out how to make their code work. – jwdonahue Sep 19 '18 at 00:25
I need just to get this Bernoulli numbers working to me get back to my principal work. Riemann Zeta's function – Adriano rox Sep 19 '18 at 00:46
@Adrianorox Yes, I understand, that's why I'm recommending that you get a package that has a built-in function for Bernoulli numbers. – Robert Dodier Sep 19 '18 at 01:15
@jwdonahue I honestly don't think the code can be fixed. The Bernoulli numbers grow in magnitude, and B[260] is already too big to even be represented by a double float (ignoring any other problems). – Robert Dodier Sep 19 '18 at 01:18
Ah, nice catch. I was waiting for input/output examples before I bothered looking at the code. – jwdonahue Sep 19 '18 at 01:57
So this is one of those problems that requires infinite precision to get very far. Basically the binary equivalent of decimal notation, just add as many digits as you need? – jwdonahue Sep 19 '18 at 02:08
No, my algorithm is giving wrong answers, I need the right numbers with something like 18~20 digits precision. In C. I know it's easier in Perl or Python, but my entire work of Zeta's function are in C. It's too late to change the language. – Adriano rox Sep 19 '18 at 02:52

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

1 Answers1