1

I'm trying to code a function that calculates the N'th value of the positive Bernouille numbers. I want to do it with this recursive formula:

What I've tried so far:

   function Get_B(N : in Integer) return Float is

     X,Bn,Bk:Float;
   begin
     if N = 0 then
       return 1.0;
     else
       Bn:=0.0;
       for K in  0..(N-1) loop
        Bk:=Get_B(K);
        X:=1.0-Float(F(N))/(Float(F(K))*Float(F(N-K))) * 
           Bk/(Float(N)-Float(K)+1.0);
        Bn:=Bn+X;
     end loop;
     return Bn;
      end if;

   end Get_B;

where F is a factorial function (that is, F(N) means N!). I think there's something wrong with the loop, but I don't know what it is.

I can't post images but here's a link to the equation (the bottom one): https://wikimedia.org/api/rest_v1/media/math/render/svg/2571d0a7024741c0a0ad0eb32cc9333f6024c528

Ludd
  • 53
  • 5
  • 1
    Why do you think there’s something wrong with the loop? – Simon Wright Oct 24 '17 at 14:07
  • I get a value for each N i put into the function, but only N=0 , N=1 and N=6 gives med the right answer. The others give wrong answers. I can't find an error in how i wrote the equation so I'm figuring it might be the loop where the problem's at. – Ludd Oct 24 '17 at 16:34
  • A quick search on wikipedia suggests the `1.0-` part should be outside the loop – egilhh Oct 24 '17 at 18:21
  • Wow, thanks Egilhh, that actually did it! Where did you find that info? – Ludd Oct 24 '17 at 20:15
  • @Ludd, https://en.wikipedia.org/wiki/Bernoulli_number search for "recursive formulas" – egilhh Oct 24 '17 at 20:21
  • I compared the results of your code to the Wikipedia results and your results for 0, 1, and all the even numbers except 2 are correct. Your result for 2 is 1.166666, should be 0.166666, and your results for odd numbers greater than 1 aren’t zero. Something strange going on ... – Simon Wright Oct 25 '17 at 16:19

1 Answers1

2

The 1.0- should be outside the loop:

function Get_B(N : in Integer) 
               return Float is

 X,Bn,Bk:Float;
 begin
   if N = 0 then
     return 1.0;
   else
     Bn:=0.0;
     for K in  0..(N-1) loop
      Bk:=Get_B(K);
      X:=Float(F(N))/(Float(F(K))*Float(F(N-K))) * 
         Bk/(Float(N)-Float(K)+1.0);
      Bn:=Bn+X;
    end loop;
    return 1.0 - Bn;
  end if;

 end Get_B;
Ludd
  • 53
  • 5
  • btw, when your input is positive numbers (including 0), you should tell that to the compiler: `function Get_B(N : in Natural)`. And if your factorial function would return Float in the first place, the expression becomes a lot more readable: `X:=F(N)/(F(K)*F(N-K)) * Bk/Float(N-K+1);` – egilhh Oct 24 '17 at 20:29
  • Thanks, I will change that! – Ludd Oct 25 '17 at 06:15