[Scilab-users] Problems arising from truncation of %pi

Dear all,

In order to test the FFT on a periodic signal whose period is an exact submultiple of the FFT length I found a frequency fo satisfying this for the chosen sample rate and window length, and generated the following signal:

x = sin(2*%pi*fo*t);

where t is a time vector. This should be a perfectly periodic discrete signal but it isn't because the sin() function has a (virtually) exact period of pi, while %pi is exact to 16 digits only.

For instance, we have (23 digits shown)

--> sin(%pi)
 ans  =
   0.0000000000000001224647

--> sin(1e10*%pi)
 ans  =
  -0.0000022393627619559233

--> sin(1e15*%pi)
 ans  =
  -0.2362090532517409080526

The Wolfram Alpha site yields the correct value 0 in all cases (using their own pi).

Questions:

1) How is the sin() function extended to very large values of the argument? My first guess would be to compute a quarter cycle using Taylor and then extend it by symmetry and periodicity, but with which period? The best approximation available is 2*%pi. Or it is possible to use extended precision internally?

2) Is there a way to get a periodic discrete sine other than extending it periodically with the desired period?
Wouldn't this create a slight glitch at the frontier between cycles?

Regards,

Federico Miyara  

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Hello Federico,

I'm not sure to understand exactly your needs but:

> De : Federico Miyara
> Envoyé : mardi 5 janvier 2021 09:20
>
> --> sin(%pi)
> ans  =
>   0.0000000000000001224647

seems to be something we cannot get rid off unless using symbolic calculation i.e. having a specific variable for pi which is recognised and gives the exact result with given functions.

> 1) How is the sin() function extended to very large values of the argument?

I don’t know about Scilab but as it is initially based on Fortran, it probably uses the related Fortran library (or maybe a more modern library, this is not documented).
You might have some answers here:

https://en.wikipedia.org/wiki/Sine#Software_implementations

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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Hello,

Le 05/01/2021 à 09:19, Federico Miyara a écrit :
> --> sin(%pi)
>   ans  =
>     0.0000000000000001224647


You face the limited precision of all numerical calculus. See
--> help %eps

It has no sense to use 23 digits with a precision of 2.22E-16, the 7th
last digits are noise.


> --> sin(1e10*%pi)
>   ans  =
>    -0.0000022393627619559233
>
> --> sin(1e15*%pi)
>   ans  =
>    -0.2362090532517409080526


All is consistent with the definition of the precision. A toy example:
--> a=2;
--> b=sqrt(2);
--> a-b^2

You expect zero because you do symbolic calculus, not Scilab.


> The Wolfram Alpha site yields the correct value 0 in all cases (using
> their own pi).


Because it uses symbolic calculus. So, if you want more precision with
Scilab you should change the standard used, but maybe some tricks in
your code can solve the problem... Or maybe you need symbolic calculus tool.

--Jean-Yves
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Hello Frederico,

Like Christophe, I am not sure this has anything to do with the implementation of sin().
It seems to be a known limitation of numerical calculations using floating numbers.
In particular, even with en hypothetical ideal value of %pi, because of the conversion to a double, %pi*1e15 is converted to the nearest available double and this "shift" or error as compared to the ideal infinitely-precise value is increasing with the size of your number.
Thus you make a x1e15 bigger error when using "%pi*1e15" than when using "%pi".
As Christophe as said, there is not much you can do, apart from resorting to symbolic calculation (what Alpha does).

You can see that by using nearfloat to get the distance between two doubles:

nearfloat("succ", %pi)-%pi
 ans  =   4.441D-16
nearfloat("succ", %pi*1e15)-%pi*1e15
 ans  =   0.5

As I said, the error on the sin() argument is getting x1e15 bigger!

I'm sorry if my explanation of floating point numbers is not really good.
This one was a revelation for me: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html .
It helped me understand why it's usually a good idea to use "reduced units" for your calculations (ie keep everything close to 1 when possible).

I hope it helps,

Cheers,

Antoine

Dear all,

In order to test the FFT on a periodic signal whose period is an exact submultiple of the FFT length I found a frequency fo satisfying this for the chosen sample rate and window length, and generated the following signal:

x = sin(2*%pi*fo*t);

where t is a time vector. This should be a perfectly periodic discrete signal but it isn't because the sin() function has a (virtually) exact period of pi, while %pi is exact to 16 digits only.

For instance, we have (23 digits shown)

--> sin(%pi)
 ans  =
   0.0000000000000001224647

--> sin(1e10*%pi)
 ans  =
  -0.0000022393627619559233

--> sin(1e15*%pi)
 ans  =
  -0.2362090532517409080526

The Wolfram Alpha site yields the correct value 0 in all cases (using their own pi).

Questions:

1) How is the sin() function extended to very large values of the argument? My first guess would be to compute a quarter cycle using Taylor and then extend it by symmetry and periodicity, but with which period? The best approximation available is 2*%pi. Or it is possible to use extended precision internally?

2) Is there a way to get a periodic discrete sine other than extending it periodically with the desired period?
Wouldn't this create a slight glitch at the frontier between cycles?

Regards,

Federico Miyara  

Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

I do not agree with your answer. Doing *by hand* the shift to a principal period  of sin () solves the problem. For example, here is a plot of sin(x) vs sin(x-floor(x/2/%pi)*%pi*2) for large values of x:

x = (10^[1:16])*%pi
sx = sin(x-floor(x/2/%pi)*%pi*2)
plot("ln",x,sin(x),'-o',x,sx,'-o')
legend("$\Large\sin x$","$\Large\sin \left(x-2\pi\left\lfloor\frac{x}{2\pi}\right\rfloor\right)$",3)

I suppose that the hardware implementation used by the compiler when using the standard math library does not use this simple trick because it would decrease the floating point performance.

Hello,

Le 05/01/2021 à 09:19, Federico Miyara a écrit :

--> sin(%pi)
  ans  =
    0.0000000000000001224647

You face the limited precision of all numerical calculus. See
--> help %eps

It has no sense to use 23 digits with a precision of 2.22E-16, the 7th last digits are noise.

--> sin(1e10*%pi)
  ans  =
   -0.0000022393627619559233

--> sin(1e15*%pi)
  ans  =
   -0.2362090532517409080526

All is consistent with the definition of the precision. A toy example:
--> a=2;
--> b=sqrt(2);
--> a-b^2

You expect zero because you do symbolic calculus, not Scilab.

The Wolfram Alpha site yields the correct value 0 in all cases (using their own pi).

Because it uses symbolic calculus. So, if you want more precision with Scilab you should change the standard used, but maybe some tricks in your code can solve the problem... Or maybe you need symbolic calculus tool.

--Jean-Yves
_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Of course, this does only work for integer multiples of %pi. This plot illustrates well Antoine's remark:

x=(10^15*%pi)+(0:0.001:2*%pi);
plot(sin(x));

I like new year's jokes !

S.

I do not agree with your answer. Doing *by hand* the shift to a principal period  of sin () solves the problem. For example, here is a plot of sin(x) vs sin(x-floor(x/2/%pi)*%pi*2) for large values of x:

x = (10^[1:16])*%pi
sx = sin(x-floor(x/2/%pi)*%pi*2)
plot("ln",x,sin(x),'-o',x,sx,'-o')
legend("$\Large\sin x$","$\Large\sin \left(x-2\pi\left\lfloor\frac{x}{2\pi}\right\rfloor\right)$",3)

I suppose that the hardware implementation used by the compiler when using the standard math library does not use this simple trick because it would decrease the floating point performance.

S.

Le 06/01/2021 à 09:14, Jean-Yves Baudais a écrit :

Hello,

Le 05/01/2021 à 09:19, Federico Miyara a écrit :

--> sin(%pi)
  ans  =
    0.0000000000000001224647

You face the limited precision of all numerical calculus. See
--> help %eps

It has no sense to use 23 digits with a precision of 2.22E-16, the 7th last digits are noise.

--> sin(1e10*%pi)
  ans  =
   -0.0000022393627619559233

--> sin(1e15*%pi)
  ans  =
   -0.2362090532517409080526

All is consistent with the definition of the precision. A toy example:
--> a=2;
--> b=sqrt(2);
--> a-b^2

You expect zero because you do symbolic calculus, not Scilab.

The Wolfram Alpha site yields the correct value 0 in all cases (using their own pi).

Because it uses symbolic calculus. So, if you want more precision with Scilab you should change the standard used, but maybe some tricks in your code can solve the problem... Or maybe you need symbolic calculus tool.

--Jean-Yves
_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Jean-Yves,

Thanks for your answer.

> You face the limited precision of all numerical calculus. See --> help
> %eps

Yes, I'm aware of this limitation, but I think an improved version of
the sine function could be provided for its use in certain cases, in a
similar fashion as the log1p() function, which takes advantage of
foating point capabilitties when the input argument is close to 1.

The function could be sinpi() or similar, with two arguments: the main
argument x and an integer argument n, being its result equivalent to

sin(x - n*pi)

where pi is the exact value of pi

> It has no sense to use 23 digits with a precision of 2.22E-16, the 7th
> last digits are noise.

I've just reported the digits provided when using format(25)

Probably it is not proper to refer to it as noise, since the difference
is deterministic as can be demonstrated by using the function nearfloat()

Regards,

Federico Miyara


--
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Antoine,

Thanks for the link, very interesting.

Regards,

Federico Miyara

Hello Frederico,

Like Christophe, I am not sure this has anything to do with the implementation of sin().
It seems to be a known limitation of numerical calculations using floating numbers.
In particular, even with en hypothetical ideal value of %pi, because of the conversion to a double, %pi*1e15 is converted to the nearest available double and this "shift" or error as compared to the ideal infinitely-precise value is increasing with the size of your number.
Thus you make a x1e15 bigger error when using "%pi*1e15" than when using "%pi".
As Christophe as said, there is not much you can do, apart from resorting to symbolic calculation (what Alpha does).

You can see that by using nearfloat to get the distance between two doubles:

nearfloat("succ", %pi)-%pi
 ans  =   4.441D-16
nearfloat("succ", %pi*1e15)-%pi*1e15
 ans  =   0.5

As I said, the error on the sin() argument is getting x1e15 bigger!

I'm sorry if my explanation of floating point numbers is not really good.
This one was a revelation for me: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html .
It helped me understand why it's usually a good idea to use "reduced units" for your calculations (ie keep everything close to 1 when possible).

I hope it helps,

Cheers,

Antoine

On 05/01/2021 09:19, Federico Miyara wrote:

Dear all,

In order to test the FFT on a periodic signal whose period is an exact submultiple of the FFT length I found a frequency fo satisfying this for the chosen sample rate and window length, and generated the following signal:

x = sin(2*%pi*fo*t);

where t is a time vector. This should be a perfectly periodic discrete signal but it isn't because the sin() function has a (virtually) exact period of pi, while %pi is exact to 16 digits only.

For instance, we have (23 digits shown)

--> sin(%pi)
 ans  =
   0.0000000000000001224647

--> sin(1e10*%pi)
 ans  =
  -0.0000022393627619559233

--> sin(1e15*%pi)
 ans  =
  -0.2362090532517409080526

The Wolfram Alpha site yields the correct value 0 in all cases (using their own pi).

Questions:

1) How is the sin() function extended to very large values of the argument? My first guess would be to compute a quarter cycle using Taylor and then extend it by symmetry and periodicity, but with which period? The best approximation available is 2*%pi. Or it is possible to use extended precision internally?

2) Is there a way to get a periodic discrete sine other than extending it periodically with the desired period?
Wouldn't this create a slight glitch at the frontier between cycles?

Regards,

Federico Miyara  

Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Stéphane,

This would be great if it worked for any x. When x is close to 1e16 there are few valid numbers per cycle. For instance

x(1) = 1e16
for k=1:5
    x(k+1) = nearfloat("succ",x(k));
end
x

yields (with format(25))

   10000000000000000. 
   10000000000000002.
   10000000000000004.
   10000000000000006.
   10000000000000008.
   10000000000000010.

Those few numbers have each a definite exact value and sin(x) has also a definite exact value which seemingly is already well approximated, unlike what I had previously supposed. For instance

[x sin(x)]

yields

   10000000000000000.   0.77968800660697878957 
   10000000000000002.  -0.8938378287657305909519
   10000000000000004.  -0.0357524369529285263036
   10000000000000006.   0.9235943558393552299535
   10000000000000008.  -0.7329493019177583112977
   10000000000000010.  -0.3135652891543322384749

I have underlined the last correct digit compared with Alpha, so it is the expected result with 16 digits.

Regards,

Federico Miyara

x = (10^[1:16])*%pi
sx = sin(x-floor(x/2/%pi)*%pi*2)
plot("ln",x,sin(x),'-o',x,sx,'-o')
legend("$\Large\sin x$","$\Large\sin \left(x-2\pi\left\lfloor\frac{x}{2\pi}\right\rfloor\right)$",3)

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Hello,

----- Original Message -----
> The function could be sinpi() or similar, with two arguments: the main
> argument x and an integer argument n, being its result equivalent to
>
> sin(x - n*pi)

So now the problem can be how these large numbers are obtained
--> a=1e16+1
--> a-1e16
of course equals zero.

> where pi is the exact value of pi

Hum... What does "exact" mean in numerical calculus? (In symbolic one, it's simpler :-)

--> format(25)
--> %pi
3.141592653589793115998

the 6 last digits are wrong (cf. https://oeis.org/A000796) because of float64 format.

> Probably it is not proper to refer to it as noise, since the difference
> is deterministic as can be demonstrated by using the function nearfloat()

Maybe noise is not the right word. But, the difference is fixed and %pi give all the time the same digit in format(25). It's a bit weird (up to my ignorance) to have these wrong deterministic digits... Hum, maybe I need to read "What Every Computer Scientist Should Know about Floating-Point Arithmetic", David Goldberg, ACM Computing Surveys, vol. 23, no  1, mars 1991 recommanded by https://en.wikipedia.org/wiki/IEEE_754

--Jean-Yves

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

The nearest double-precision IEEE-754 binary floating-point number for the decimal number PI
 
3.141592653589793 23846264338327950288419716939937510582097494459230 78164 06286 ....
 
is

3.141592653589793 115997963468544185161590576171875

It can be shown this way: its internal base 2 representation is (I added spaces to separate the sign, exponent and mantissa) is given in Scilab by

--> bitstring(%pi)
 ans  =

  "0 10000000000 1001001000011111101101010100010001000010110100011000"

i.e. a positive number, with exponent 2^(1024-1023) = 1 and mantissa

1.1001001000011111101101010100010001000010110100011000

i.e. the leading "1."  is implicit (normalized convention) and has to be added. Hence, the final numer, after mutiplication by 2^exponent, is

11.001001000011111101101010100010001000010110100011000

You won't be able to convert this back to decimal by using floating point arithmetic. There is an efficient converter at  https://www.exploringbinary.com/binary-converter/ which gives the actual result

S.

Hello,

----- Original Message -----

The function could be sinpi() or similar, with two arguments: the main
argument x and an integer argument n, being its result equivalent to

sin(x - n*pi)

So now the problem can be how these large numbers are obtained
--> a=1e16+1
--> a-1e16
of course equals zero. 

where pi is the exact value of pi

Hum... What does "exact" mean in numerical calculus? (In symbolic one, it's simpler :-)

--> format(25)
--> %pi
3.141592653589793115998

the 6 last digits are wrong (cf. https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/oeis.org/A000796) because of float64 format.

Probably it is not proper to refer to it as noise, since the difference
is deterministic as can be demonstrated by using the function nearfloat()

Maybe noise is not the right word. But, the difference is fixed and %pi give all the time the same digit in format(25). It's a bit weird (up to my ignorance) to have these wrong deterministic digits... Hum, maybe I need to read "What Every Computer Scientist Should Know about Floating-Point Arithmetic", David Goldberg, ACM Computing Surveys, vol. 23, no  1, mars 1991 recommanded by https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/en.wikipedia.org/wiki/IEEE_754

--Jean-Yves

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

But, of course, we have

--> bitstring(3.141592653589793) == bitstring(%pi)
 ans  =

  T

which means that 3.141592653589793 and 3.141592653589793115997963468544185161590576171875 have the same internal IEEE754 representation.

S.

The nearest double-precision IEEE-754 binary floating-point number for the decimal number PI
 
3.141592653589793 23846264338327950288419716939937510582097494459230 78164 06286 ....
 
is

3.141592653589793 115997963468544185161590576171875

It can be shown this way: its internal base 2 representation is (I added spaces to separate the sign, exponent and mantissa) is given in Scilab by

--> bitstring(%pi)
 ans  =

  "0 10000000000 1001001000011111101101010100010001000010110100011000"

i.e. a positive number, with exponent 2^(1024-1023) = 1 and mantissa

1.1001001000011111101101010100010001000010110100011000

i.e. the leading "1."  is implicit (normalized convention) and has to be added. Hence, the final numer, after mutiplication by 2^exponent, is

11.001001000011111101101010100010001000010110100011000

You won't be able to convert this back to decimal by using floating point arithmetic. There is an efficient converter at  https://www.exploringbinary.com/binary-converter/ which gives the actual result

S.

Le 08/01/2021 à 10:48, Jean-Yves Baudais a écrit :

Hello,

----- Original Message -----

The function could be sinpi() or similar, with two arguments: the main
argument x and an integer argument n, being its result equivalent to

sin(x - n*pi)

So now the problem can be how these large numbers are obtained
--> a=1e16+1
--> a-1e16
of course equals zero. 

where pi is the exact value of pi

Hum... What does "exact" mean in numerical calculus? (In symbolic one, it's simpler :-)

--> format(25)
--> %pi
3.141592653589793115998

the 6 last digits are wrong (cf. https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/oeis.org/A000796) because of float64 format.

Probably it is not proper to refer to it as noise, since the difference
is deterministic as can be demonstrated by using the function nearfloat()

Maybe noise is not the right word. But, the difference is fixed and %pi give all the time the same digit in format(25). It's a bit weird (up to my ignorance) to have these wrong deterministic digits... Hum, maybe I need to read "What Every Computer Scientist Should Know about Floating-Point Arithmetic", David Goldberg, ACM Computing Surveys, vol. 23, no  1, mars 1991 recommanded by https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/en.wikipedia.org/wiki/IEEE_754

--Jean-Yves

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

----- Original Message -----
> From: "Stéphane Mottelet"
> --> bitstring(%pi)
> ans =
>
> "0 10000000000 1001001000011111101101010100010001000010110100011000"
> [...]
> https://www.exploringbinary.com/binary-converter/


Thank you for the recall on floating point representation and for the link (the French wiki page on IEEE 754 gives also a full and clear example). Of course we can convert float64 from finite binary representation to decimal one with more than 16 digits, but in that case the binary to float tranform is not a "relation binaire univoque" and the finite number representation problem is "put under the carpet". As I prefer one-to-one convertion and I hate carpet, I use 16 digits at most to write float64.

-- Jean-Yves
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

----- Original Message -----

From: "Stéphane Mottelet" 
--> bitstring(%pi)
ans =

"0 10000000000 1001001000011111101101010100010001000010110100011000"
[...]
https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/www.exploringbinary.com/binary-converter/

Thank you for the recall on floating point representation and for the link (the French wiki page on IEEE 754 gives also a full and clear example). Of course we can convert float64 from finite binary representation to decimal one with more than 16 digits, but in that case the binary to float tranform is not a "relation binaire univoque" and the finite number representation problem is "put under the carpet". As I prefer one-to-one convertion and I hate carpet, I use 16 digits at most to write float64.

There is no hiding of anything in my demonstration. I was just showing that the figures displayed by Scilab after the 16th are computed on a clear basis. This is a legacy choice of Scilab to allow the display of more than 16 decimals and this can be discussed. For example, "long" format of Matlab displays the shortest decimal number matching the internal IEEE754 representation.

S.

-- Jean-Yves

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Jean-Yves,

>> sin(x - n*pi)
> So now the problem can be how these large numbers are obtained
> --> a=1e16+1
> --> a-1e16
> of course equals zero.

Yes, I've thought about it and you are right, above 1e16 x is so sparse,
cycle-wise speaking, that my original intention doesn't make much sense.

Besides, at a sample rate of 44100, and a maximum frequency of 22050 Hz,
a phase of 1e16 is reached only  after 2288 years, so it's not a
practical problem!

I must confess, by the way, that because of a stupid interpretation of
some results on my part, I thought an error I was having was due to the
slow departure from strict discrete signal periodicity but it wasn't. Sorry!

Regards,

Federico Miyara

--
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Stéphane,

Thanks for your comments, I didn't know bitstring().

Regards,

Federico Miyara

Le 08/01/2021 à 13:28, Jean-Yves Baudais a écrit :

----- Original Message -----

From: "Stéphane Mottelet" 
--> bitstring(%pi)
ans =

"0 10000000000 1001001000011111101101010100010001000010110100011000"
[...]
https://antispam.utc.fr/proxy/2/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/www.exploringbinary.com/binary-converter/

Thank you for the recall on floating point representation and for the link (the French wiki page on IEEE 754 gives also a full and clear example). Of course we can convert float64 from finite binary representation to decimal one with more than 16 digits, but in that case the binary to float tranform is not a "relation binaire univoque" and the finite number representation problem is "put under the carpet". As I prefer one-to-one convertion and I hate carpet, I use 16 digits at most to write float64.

There is no hiding of anything in my demonstration. I was just showing that the figures displayed by Scilab after the 16th are computed on a clear basis. This is a legacy choice of Scilab to allow the display of more than 16 decimals and this can be discussed. For example, "long" format of Matlab displays the shortest decimal number matching the internal IEEE754 representation.

S.

-- Jean-Yves

_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Jean-Yves,

sin(x - n*pi)

So now the problem can be how these large numbers are obtained
--> a=1e16+1
--> a-1e16
of course equals zero.

Yes, I've thought about it and you are right, above 1e16 x is so sparse, cycle-wise speaking, that my original intention doesn't make much sense.

MPScilab is an external module available up to Scilab 5.3 -- not recompiled since then -- allowing arbitrary precision computations for a short list of the mathematical functions.
sin() is not in the list, but exp() is. So if this release supports complex input arguments, Euler formulae should be usable to comput trigonometric functions as well in arbitrary precision.

Unlike MPScilab, DD_QD is available for Scilab 6.1. It implements floating point computations for decimal numbers encoded over 2x8 or 4x8 bytes, instead of "only" 8 bytes.
sin(), cos() and tan() are among supported functions.

You may try them out, out of curiosity.

Regards
Samuel

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

Hello Federico,

Le 09/01/2021 à 01:32, Federico Miyara a écrit :

Jean-Yves,

sin(x - n*pi)

So now the problem can be how these large numbers are obtained
--> a=1e16+1
--> a-1e16
of course equals zero.

Yes, I've thought about it and you are right, above 1e16 x is so sparse, cycle-wise speaking, that my original intention doesn't make much sense.

MPScilab is an external module available up to Scilab 5.3 -- not recompiled since then -- allowing arbitrary precision computations for a short list of the mathematical functions.
sin() is not in the list, but exp() is.

Sorry: all trigonometric functions are supported. I have just missed the dedicated section where they are addressed.

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users