### Re: [Scilab-users] Bernoulli numbers calculation


Hello Lester,

Le 30/12/2021 à 08:59, Lester Anderson a écrit :

Hello Samuel,

Thanks for the solution. As pointed out it is best to show the
equation being assessed (from www.bernoulli.org
). The one I looked at was the following:

Explicit_formula.PNG

Using nchoosek in the original code gives the same issue.

The inner sum over v is very prone to cumulative rounding errors:
The term v^n gets huge rapidly -- so is numerically truncated --,
while the (-1)^v term makes the sum alternate, which enhances residues...
that then mainly come from numerical truncations.
With n_max = 20, the maximum value of nchoosek(20,10)=184756
and is definitely not an issue. While even only 13^13
--> 13^13
ans  =
3.029D+14
is already not far from 1/%eps.
The recurrent implementation proposed earlier and based on
Bₘ =  -∑_k=0_→ m-1 (C_^k _m+1 ).B_k /(m+1)

\frac{-1}{m+1}{\sum _{{k=0}}^{{m-1}}{m+1 \choose k}B_k

has neither alternate terms nor huge power values that make a
direct computation numerically catastrophic.

Regards
Samuel
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users



### Re: [Scilab-users] Bernoulli numbers calculation

Hello Samuel,

Thanks for the solution. As pointed out it is best to show the equation
being assessed (from www.bernoulli.org). The one I looked at was the
following:

[image: Explicit_formula.PNG]

Using nchoosek in the original code gives the same issue.

Lester

On Wed, 29 Dec 2021 at 22:46, Samuel Gougeon  wrote:

> Hello Lester,
>
>
> Le 29/12/2021 à 09:00, Lester Anderson a écrit :
>
> Hello,
>
> A quick query. Have adapted existing Matlab code for Scilan to calculate
> Bernoulli numbers using an explicit formula (www.bernoulli.org).
>
>
> Not so explicit. Could you please provide the formula, before implementing
> it in Scilab language?
> The following is a "vectorized" version of your code, likely with the same
> mistake:
>
> B=[];for m = 0:20
> Sum = 0;
> for k=0:m
> v = 0:k;
> Sum = Sum + sum(((-1).^v).*nchoosek(k,v).*(v.^m)/(k+1));
> end
> B(m+1) = Sum;end
>
> A working implementation based on the recurrent formula Bm=−1m+1∑k=0m−1(m+
> 1k)Bk
>
> could be:
>
> mMax = 20;B = zeros(1,mMax);B(1:2) = [1 -1/2];for m = 2:2:mMax
> B(1,m+1) = -sum(nchoosek(m+1,0:m-1).*B(1,1:m))/(m+1);end
>
> --> B
>  ans  =
>  column 1 to 13
>1.  -0.5   0.167   0.  -0.033   0.   0.0238095   0.  -0.033
> 0.   0.0757576   0.  -0.2531136
>
>
>  column 14 to 21
>0.   1.167   0.  -7.0921569   0.   54.971178   0.  -529.12424
>
> To fasten the loop, the call to nchoosek() could be replaced with the 3
> last rows of Pascal's triangle.
>
> Regards
> Samuel
> ___
> users mailing list
> users@lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
>
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users



### Re: [Scilab-users] Bernoulli numbers calculation


Hello Lester,

Le 29/12/2021 à 09:00, Lester Anderson a écrit :

Hello,

A quick query. Have adapted existing Matlab code for Scilan to
calculate Bernoulli numbers using an explicit formula
(www.bernoulli.org ).

Not so explicit. Could you please provide the formula, before
implementing it in Scilab language?
The following is a "vectorized" version of your code, likely with the
same mistake:

B=[];
for  m  =  0:20
Sum  =  0;
for  k=0:m
v  =  0:k;
Sum  =  Sum  +  sum(((-1).^v).*nchoosek(k,v).*(v.^m)/(k+1));
end
B(m+1)  =  Sum;
end

A working implementation based on the recurrent formula
Bm=−1m+1∑k=0m−1(m+1k)Bk

could be:

mMax  =  20;
B  =  zeros(1,mMax);
B(1:2)  =  [1  -1/2];
for  m  =  2:2:mMax
B(1,m+1)  =  -sum(nchoosek(m+1,0:m-1).*B(1,1:m))/(m+1);
end--> B ans  = column 1 to 13   1.  -0.5   0.167   0.
-0.033   0.   0.0238095   0.  -0.033   0.   0.0757576   0.
-0.2531136

column 14 to 21
0.   1.167   0.  -7.0921569   0.   54.971178   0. -529.12424

To fasten the loop, the call to nchoosek() could be replaced with the 3
last rows of Pascal's triangle.

Regards
Samuel
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users



### Re: [Scilab-users] Bernoulli numbers calculation

Hello Lester,

> De : Lester Anderson
> Envoyé : mercredi 29 décembre 2021 09:00
>
> Have adapted existing Matlab code for Scilan to calculate Bernoulli
> numbers using an explicit formula [...]  11 onwards the values are incorrect:

I don't know what you already investigated but did you consider the possibility
of overflow or underflow?

You may already know these documents but if need be:

Michaël Baudin, /Scilab is not naïve/, 2010
https://www.scilab.org/sites/default/files/scilabisnotnaive.pdf

Michaël Baudin, /Floating point numbers in Scilab/, 2011
https://www.scilab.org/sites/default/files/floatingpoint-so_v0.4.pdf

Concerning the binomial coefficient, you should maybe use

binomial(p, k)(\$)

Hope this helps

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are
not the intended recipient (or have received this e-mail in error), please
notify the sender immediately and destroy this e-mail. Any unauthorized
copying, disclosure or distribution of the material in this e-mail is strictly
forbidden.
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users



### Re: [Scilab-users] Bernoulli numbers calculation

Hello,

Concerning my former message,
it seems that binomial() is no longer documented (although it works).

There is the function nchoosek() that is documented, with a diagram showing the
overflow limit:

https://help.scilab.org/docs/6.1.1/en_US/nchoosek.html

regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are
not the intended recipient (or have received this e-mail in error), please
notify the sender immediately and destroy this e-mail. Any unauthorized
copying, disclosure or distribution of the material in this e-mail is strictly
forbidden.
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users



### Re: [Scilab-users] Bernoulli numbers calculation

Sorry for the spam,

I'd better collect all the information before posting ^_^

To be clear:

The function binomial() is for the binomial distribution.

The binomial coefficient is obtained with nchoosek().

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are
not the intended recipient (or have received this e-mail in error), please
notify the sender immediately and destroy this e-mail. Any unauthorized
copying, disclosure or distribution of the material in this e-mail is strictly
forbidden.
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users