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 |
Dang Ngoc Chan, Christophe |
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 |
In reply to this post by fmiyara
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 |
Antoine Monmayrant |
In reply to this post by fmiyara
Hello Frederico, Like Christophe, I am not sure this has anything to do with the
implementation of sin(). You can see that by using nearfloat to get the distance between two doubles:
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. I hope it helps, Cheers, On 05/01/2021 09:19, Federico Miyara
wrote:
_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users |
In reply to this post by jbaudais
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 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, -- 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); I like new year's jokes ! S. Le 06/01/2021 à 09:49, Stéphane
Mottelet a écrit :
-- 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 |
In reply to this post by jbaudais
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 > > > >> --> 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 > > -- 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 |
In reply to this post by Antoine Monmayrant
Antoine, Thanks for the link, very interesting. Regards, Federico Miyara On 06/01/2021 05:22, Antoine Monmayrant
wrote:
_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users |
In reply to this post by mottelet
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 On 06/01/2021 05:49, Stéphane Mottelet
wrote:
x = (10^[1:16])*%pi _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users |
In reply to this post by fmiyara
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) 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 piHum... 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) which means that 3.141592653589793 and
3.141592653589793115997963468544185161590576171875 have the same
internal IEEE754 representation. S. Le 08/01/2021 à 12:20, Stéphane
Mottelet a écrit :
The nearest double-precision IEEE-754 binary floating-point number for the decimal number PI -- 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 |
In reply to this post by mottelet
----- 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 |
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 -- 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 |
In reply to this post by jbaudais
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 |
In reply to this post by mottelet
Stéphane, Thanks for your comments, I didn't know bitstring(). Regards, Federico Miyara On 08/01/2021 10:48, Stéphane Mottelet
wrote:
_______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users |
Samuel GOUGEON |
In reply to this post by fmiyara
Hello Federico,
Le 09/01/2021 à 01:32, Federico Miyara
a écrit :
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. _______________________________________________ users mailing list [hidden email] http://lists.scilab.org/mailman/listinfo/users |
Samuel GOUGEON |
Le 09/01/2021 à 13:04, Samuel Gougeon a
écrit :
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 |
Free forum by Nabble | Edit this page |