Re: [Scilab-users] Bernoulli numbers calculation

2021-12-30 Thread Samuel Gougeon 

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

2021-12-30 Thread Lester Anderson 
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

2021-12-29 Thread Samuel Gougeon 

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

2021-12-29 Thread Dang Ngoc Chan, Christophe 
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

2021-12-29 Thread Dang Ngoc Chan, Christophe 
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

2021-12-29 Thread Dang Ngoc Chan, Christophe 
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