Floting point "difference" between FPC vs OCTAVE lazarus / 21-General ---------------------------------------------------------------------- juanirias · 2022-11-18 14:21 Hi! I think there is something wrong with the arithmetic of FPC or OCTAVE; but the result is different. I attach an image and the extracts of the codes. Octave code: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i % -------- gstime       - greenwich sidereal time    year = 1992;    mon = 8;    day = 20;    hr = 12;    min = 14;    sec = 0.00;    [jdut1, jdut1frac] = jday( year,mon,day,hr,min,sec );     fprintf(1,'jd %18.10f  %18.10f   %18.10f \n\n',jdut1, jdut1frac, jdut1+jdut1frac);     % ---- gst = gstime(jdut1);    tut1= ( jdut1+jdut1frac - 2451545.0 ) / 36525.0        temp = - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1  ...        + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 67310.54841     % 360/86400 = 1/240, to deg, to rad    %       temp = rem( temp*rad/240.0,twopi )    temp = rem( temp,86400 )     temp = temp / 240.0     % ------------------------ check quadrants --------------------    if ( temp < 0.0 )        temp = temp + 360;  end     gst = temp      fprintf(1,'gst = %16.12f\n', gst );   PASCAL code Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program fallo;uses asttime,astmath;   CONST     //Pi      : EXTENDED =     Pi;  { 6.28318530717959; }     Deg2Rad    : EXTENDED =     Pi/180.0;     Rad    : EXTENDED =     180.0/Pi;   VAR     Temp, TUT1 : EXTENDED;   Year,Mon,Day,Hr,MIN                   : INTEGER;      SEC, JDUT1, Lon, LST, GST             : EXTENDED;   BEGIN          Year   := 1992;       Mon    := 08;       Day    := 20;       Hr     := 12;     MIN    := 14;       SEC    := 0.00;              Lon    := -104/Rad; // 104 deg W              writeln('year ',year:4);       writeln('mon  ',mon:2);       writeln('day  ',day:2);       writeln(hr:2,' : ',min:2,' : ',sec:2:2);              JDAY( Year,Mon,Day,Hr,MIN,SEC, JDUT1);              writeln;       writeln('JD:   ',JDUT1:9:10);                {*    tut1= ( jdut1+jdut1frac - 2451545.0 ) / 36525.0     temp = - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1  ...      + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 67310.54841        temp = rem( temp,86400 )     temp = temp / 240.0 *}        TUT1:= ( JDUT1 - 2451545.0 ) / 36525.0;     writeln('tut1: ',TUT1:9:10);     Temp:= - 6.2E-6*TUT1*TUT1*TUT1            + 0.093104*TUT1*TUT1          + (876600.0*3600 + 8640184.812866)*TUT1            + 67310.54841;  {sec }     writeln('Temp: ',Temp:9:10);;     ;  // Temp:= MODFUNC( Temp*Deg2Rad/240.0,TwoPi ); { 360/86400 = 1/240, to deg, to rad }     Temp := modfunc(temp,86400);     writeln('Temp: ',Temp:9:10);     Temp := Temp/240;     writeln('Temp: ',Temp:9:10);     if ( temp < 0.0 ) then temp := temp + 360;     writeln('Temp: ',Temp:9:10);       { ------------------------ Check quadrants --------------------- }//     IF Temp < 0.0 THEN//         Temp:= Temp + TwoPi;        GST:= Temp*Deg2Rad;                                                    writeln;       //writeln(' LST: ',LST*Rad:3:9);       writeln(' GST: ',GST*Rad:3:9);         END.    The "Book" show same values than OCTAVE !!! see the image below Thank you! ---------------------------------------------------------------------- marcov · 2022-11-18 15:41 The JD is different, and off by an half.  The original definition of Julian starts at noon, while many computer julian variants place it at midnight (to be more in line with normal date transitions). ---------------------------------------------------------------------- juanirias · 2022-11-18 15:59 Hi Marco, maybe can review the equations: JDAY pascal Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i PROCEDURE JDAY          ( Year,Mon,Day,Hr,MIN                   : INTEGER;                             SEC                                   : EXTENDED;                               VAR JD                              : EXTENDED );   BEGIN     JD:= 367.0 * Year          - TRUNC( (7* (Year+TRUNC ( (Mon+9)/12) ) ) * 0.25 )          + TRUNC( 275*Mon / 9 )          + Day + 1721013.5          + ( (SEC/60.0 + MIN ) / 60.0 + Hr ) / 24.0;  { UT in days }          { - 0.5*SGN(100.0*Year + Mon - 190002.5) + 0.5; }   END;  { PROCEDURE JDAY }{|    Vallado       2007, 190} function [jd, jdfrac] = jday(yr, mon, day, hr, min, sec)         % ------------------------  implementation   ------------------      jd = 367.0 * yr  ...             - floor( (7 * (yr + floor( (mon + 9) / 12.0) ) ) * 0.25 )   ...             + floor( 275 * mon / 9.0 ) ...             + day + 1721013.5;   % use - 678987.0 to go to mjd directly        jdfrac = (sec + min * 60.0 + hr *3600.0) / 86400.0;              % check jdfrac        if jdfrac > 1.0             jd = jd + floor(jdfrac);            jdfrac = jdfrac - floor(jdfrac);        end  Octave JDAY Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function [jd, jdfrac] = jday(yr, mon, day, hr, min, sec)         % ------------------------  implementation   ------------------        jd = 367.0 * yr  ...             - floor( (7 * (yr + floor( (mon + 9) / 12.0) ) ) * 0.25 )   ...             + floor( 275 * mon / 9.0 ) ...             + day + 1721013.5;   % use - 678987.0 to go to mjd directly        jdfrac = (sec + min * 60.0 + hr *3600.0) / 86400.0;              % check jdfrac        if jdfrac > 1.0               jd = jd + floor(jdfrac);            jdfrac = jdfrac - floor(jdfrac);      end  please review the image regards, Juan ---------------------------------------------------------------------- Thaddy · 2022-11-18 16:12 Marco is right, so is freepascal, but Octave is indeed wrong. See https://en.wikipedia.org/wiki/Julian_day ---------------------------------------------------------------------- juanirias · 2022-11-18 16:23 Thanks for the comments. One remark: Im not using Julian Day from Octave, the functions used are shown above, they are basically the same; maybe you can review it Regards Juan ---------------------------------------------------------------------- Thaddy · 2022-11-18 16:43 The Julian date from Octave is actually Modified Julian Date (MJD). https://core2.gsfc.nasa.gov/time/ and https://en.wikipedia.org/wiki/Julian_day Functions to handle MJD are declared but not implemented - yet -in Dateutils. There are also well tested routines to handle julian dates, you might prefer those over your own. https://www.freepascal.org/docs-html/rtl/dateutils/index-5.html ---------------------------------------------------------------------- juanirias · 2022-11-18 19:15 I'm not using Octave function, and is nor MJD are JD. ---------------------------------------------------------------------- juanirias · 2022-11-18 19:38 The functions are now the "same", please review: PASCAL Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i PROCEDURE JDAY          ( Year,Mon,Day,Hr,MIN                   : INTEGER;                             SEC                                   : EXTENDED;                               VAR JD                       : EXTENDED );var jdfrac: extended;   BEGIN     JD:= 367.0 * Year        - TRUNC( (7* (Year+TRUNC ( (Mon+9)/12) ) ) * 0.25 )          + TRUNC( 275*Mon / 9 )          + Day + 1721013.5;          jdfrac := ( sec + min*60 + hr*3600 )/86400;          if (jdfrac > 1.0) then          begin            jd := jd + trunc(jdfrac);            jdfrac := jdfrac - trunc(jdfrac);            end;            jd:=jd+jdfrac;                // + ( (SEC/60.0 + MIN ) / 60.0 + Hr ) / 24.0;  { UT in days }          { - 0.5*SGN(100.0*Year + Mon - 190002.5) + 0.5; }     writeln('JD:',JD);   END;  { PROCEDURE JDAY }  Octave Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i function [jd, jdfrac] = jday(yr, mon, day, hr, min, sec)         % ------------------------  implementation   ------------------        jd = 367.0 * yr  ...             - floor( (7 * (yr + floor( (mon + 9) / 12.0) ) ) * 0.25 )   ...             + floor( 275 * mon / 9.0 ) ...             + day + 1721013.5;   % use - 678987.0 to go to mjd directly        jdfrac = (sec + min * 60.0 + hr *3600.0) / 86400.0;              % check jdfrac        if jdfrac > 1.0             jd = jd + floor(jdfrac);            jdfrac = jdfrac - floor(jdfrac);        end        please review the Image above... regards, Juan ---------------------------------------------------------------------- juanirias · 2022-11-18 20:21 I think Pascal has better precision! In octave: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i >> epsans = 2.220446049250313e-16  In pascal: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i eps: 1.08420217248550443401E-0019 ---------------------------------------------------------------------- MarkMLl · 2022-11-18 20:36 Quote from: juanirias on November 18, 2022, 09:21:26 pmI think Pascal has better precision! There's lots of ways that this numeric stuff can go wrong. The transcript at https://github.com/MarkMLl/NOAA_AstroDemo might interest. I don't claim to have good enough Numerical Analysis skills to have worked through it in detail, but I can say that in practice it produced results that were a good match to observed phenomena. MarkMLl ---------------------------------------------------------------------- juanirias · 2022-11-18 23:57 Thanks! Maybe can review the last pascal code, and run for share your machine eps. regards, Juan ---------------------------------------------------------------------- tetrastes · 2022-11-19 00:00 For purity I changed Trunc to Floor and extended to double. Also I don't know where to find asttime and astmath units, so code stops before MODFUNC, but the point of interest is just before it: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program project1; uses math;//, asttime,astmath; PROCEDURE JDAY          ( Year,Mon,Day,Hr,MIN : INTEGER;                               SEC  : double; //EXTENDED;                               VAR JD : double); //EXTENDED );varjdfrac: double; //extended;   BEGIN     JD:= 367.0 * Year          - Floor( (7* (Year+ Floor( (Mon+9)/12) ) ) * 0.25 )        + Floor( 275.0*Mon / 9 )          + Day + 1721013.5;          jdfrac := ( sec + min*60.0 + hr*3600.0 )/86400.0;          if (jdfrac > 1.0) then            begin            jd := jd + Floor(jdfrac);          jdfrac := jdfrac - Floor(jdfrac);            end;            jd:=jd+jdfrac;           // + ( (SEC/60.0 + MIN ) / 60.0 + Hr ) / 24.0;  { UT in days }          { - 0.5*SGN(100.0*Year + Mon - 190002.5) + 0.5; }     writeln('JD:',JD);   END;  { PROCEDURE JDAY } CONST  //Pi      : EXTENDED =     Pi;  { 6.28318530717959; }  Deg2Rad    : {EXTENDED} double =     Pi/180.0;  Rad    : {EXTENDED} double =     180.0/Pi;VAR  Temp, TUT1 : double; //EXTENDED;  Year,Mon,Day,Hr,MIN                   : INTEGER;   SEC, JDUT1, Lon, LST, GST             : double;  //EXTENDED;BEGIN     Year   := 1992;  Mon    := 08;    Day    := 20;    Hr     := 12;    MIN    := 14;    SEC    := 0.00;     Lon    := -104/Rad; // 104 deg W     writeln('year ',year:4);    writeln('mon  ',mon:2);    writeln('day  ',day:2);    writeln(hr:2,' : ',min:2,' : ',sec:2:2);     JDAY( Year,Mon,Day,Hr,MIN,SEC, JDUT1);     writeln;    writeln('JD:   ',JDUT1:9:10);      TUT1:= ( JDUT1 - 2451545.0 ) / 36525.0;  writeln('tut1: ',TUT1{:9:10});  Temp:= - 6.2E-6*TUT1*TUT1*TUT1         + 0.093104*TUT1*TUT1         + (876600.0*3600 + 8640184.812866)*TUT1         + 67310.54841;  {sec }  writeln('Temp: ',Temp:9:10); // Temp:= MODFUNC( Temp*Deg2Rad/240.0,TwoPi ); { 360/86400 = 1/240, to deg, to rad }// Temp := modfunc(temp,86400);     readln;END.  Compiled with fpc 3.2.2 and the output is like OP's (see Temp): Code: Bash  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i year 1992mon   8day  2012 : 14 : 0.00JD: 2.4488550097222221E+006 JD:   2448855.0097222221tut1: -7.3647919993919009E-002Temp: -232984190.5178592500  I rewrote this code to C Code: C  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i #include #define _USE_MATH_DEFINES#include  void jday(int year, int mon, int day, int hr, int min, double sec, double *jd){      double jdfrac; //extended;      *jd = 367.0 * year          - floor( (7* (year+ floor( (mon+9)/12) ) ) * 0.25 )          + floor( 275*mon / 9 )          + day + 1721013.5;          jdfrac = ( sec + min*60 + hr*3600 )/86400;          if (jdfrac > 1.0)          {            *jd = *jd + floor(jdfrac);            jdfrac = jdfrac - floor(jdfrac);        }          *jd = *jd + jdfrac;      printf("jd: %.10lf\n", *jd);} const double deg2rad = M_PI/180.0;const double rad = 180.0/M_PI; void main(void){  double temp, tut1, dummy;  int year,mon,day,hr,min;  double sec, jdut1, lon, gst;     year   = 1992;    mon    = 8;    day    = 20;    hr     = 12;    min    = 14;    sec  = 0.0;     lon    = -104/rad; // 104 deg w     printf("year %d\n",year);    printf("mon  %d\n",mon);    printf("day  %d\n",day);  printf("%d : %d : %2.2f\n", hr, min, sec);     jday( year,mon,day,hr,min,sec, &jdut1);    printf("\njd:   %.10lf\n",jdut1);   tut1 = ( jdut1 - 2451545.0 ) / 36525.0;  printf("tut1: %.10lf\n",tut1);  temp = - 6.2e-6*tut1*tut1*tut1         + 0.093104*tut1*tut1         + (876600.0*3600 + 8640184.812866)*tut1         + 67310.54841;  //sec   printf("temp: %.10lf\n",temp);   // temp:= modfunc( temp*deg2rad/240.0,twopi ); { 360/86400 = 1/240, to deg, to rad }  // Temp := modfunc(temp,86400);    // modf(temp/86400.0, &temp); - I tried this, but it's not modfunc }  and compiled with MS VS 2022 cl.exe. The result is as Octave's (see temp): Code: Bash  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i year 1992mon  8day  2012 : 14 : 0.00jd: 2448855.0097222221 jd:   2448855.0097222221tut1: -0.0736479200temp: -232984181.0909254849  So there is smth wrong with fpc. ---------------------------------------------------------------------- juanirias · 2022-11-19 00:05 Dear Tetrastes MODFUNC Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i FUNCTION MODFUNC             ( XVal,Modby                            : EXTENDED ) : EXTENDED;   VAR       TempValue: EXTENDED;   BEGIN     TempValue := XVal;     WHILE ABS(TempValue) > ModBy DO         TempValue:= TempValue - INT(XVal/ModBy)*ModBy;     MODFUNC:= TempValue;   END;  { FUNCTION MODFUNC }  Have same issue C++ same as Octave (but big eps). Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i # include # include # include  templatefloat_t machine_eps(){        union        {              float_t f;                int_t   i;        } one, one_plus, little, last_little;         one.f    = 1.0;        little.f = 1.0;        last_little.f = little.f;         while(true)        {              one_plus.f = one.f;                one_plus.f += little.f;                 if( one.i != one_plus.i )                {                      last_little.f = little.f;                        little.f /= 2.0;                }                else                {                        return last_little.f;                }      }} int main(){        std::cout << "machine epsilon:\n";        std::cout << "float: " << std::setprecision(18)<< machine_eps() << std::endl;        std::cout << "double: " << std::setprecision(18) << machine_eps() << std::endl;}   Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program epsilon_test;varepsilon : extended;begin epsilon := 1.0; while (1.0 + 0.5 * epsilon) <> 1.0 do    epsilon := 0.5 * epsilon;    writeln('eps:',epsilon);    end.   regards, Juan ---------------------------------------------------------------------- tetrastes · 2022-11-19 00:34 For Extended eps: 1.08420217248550443401E-0019 For double eps: 2.2204460492503131E-016 But results for Temp by fpc are not different: Temp: -232984190.5178493540 for extended Temp: -232984190.5178592500 for double while by Octave and cl.exe temp: -232984181.0909255 And I trust cl.exe... ---------------------------------------------------------------------- juanirias · 2022-11-19 00:36 Exactly! I feeling wire with this.... thanks for your time ---------------------------------------------------------------------- uart · 2022-11-19 00:46 Back to the original problem, where your "temp" variables differ after only 7 significant digits, I've tested the relevant code and it does appear that  your fpc code is the one with the precision issue here. Hopefully someone else can explain specifically why the problem occurs, but I've isolated it very specifically to the literal (876600.0 * 3600.0) in your first "Temp" expression. Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i temp := - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1        + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 67310.54841; Modifying the code as follows results in correct precision. Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i var ci1 : Int64 .... ci1 := 876600 * 3600; temp := - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1        + (ci1 + 8640184.812866) * tut1 + 67310.54841; ---------------------------------------------------------------------- uart · 2022-11-19 00:52 On further testing, simply dropping the decimal points from the integer part of the expression is enough to fix it. Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i     temp := - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1            + (876600 * 3600 + 8640184.812866) * tut1 + 67310.54841; The above gives the expected precision. ---------------------------------------------------------------------- juanirias · 2022-11-19 01:09 Dear UART Thanks a lot! ---------------------------------------------------------------------- uart · 2022-11-19 01:10 Here is the simplest code to demonstrate the problem. Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i var temp1, temp2 : double;...  temp1 := 876600*3600 + 8640184.812866;  temp2 := 876600.0*3600.0 + 8640184.812866;  writeln(temp1,'  ',temp2); Variable "temp2" ends up with only 7 significant digits of precision. ---------------------------------------------------------------------- uart · 2022-11-19 01:38 Quote from: juanirias on November 19, 2022, 02:09:03 amThanks a lot! No problems. I know that's fixed your immediate problem, but I'm still interested in why (is it a bug or is it expected behavior)? BTW. The 7 significant digit thing sounds suspiciously like single precision, and indeed when I check it does appear that the "3600.0*876600.0" part of the expression is evaluated at only single precision. In the following code, both temp2 and temp3 result in the same level of reduced precision. Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i  var temp1, temp2, temp3 : double;      cs1 : single;    ...      temp1 := 876600*3600 + 8640184.812866;      temp2 := 876600.0*3600.0 + 8640184.812866;      cs1 := 876600*3600;      temp3 := cs1 + 8640184.812866;      writeln(temp1,'  ',temp2,'  ',temp3); ---------------------------------------------------------------------- uart · 2022-11-19 02:31 So I'm assuming that this happens because the literals in the expression get the smallest representation that can hold them exactly (if there is a smaller representation that can hold them exactly). So both 3600.0 and 876600.0 can be represented exactly in single precision, but unfortunately their product cannot, so the problem arises when (3600.0 * 876600.0) is evaluated in single precision. Hopefully when the gurus come back online someone can tell us if this is expected behavior or not. ---------------------------------------------------------------------- juanirias · 2022-11-19 19:58 Thanks a lot! the final answer is very similar now. regards, Juan ---------------------------------------------------------------------- MarkMLl · 2022-11-19 20:34 On a whim I've just tried this: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program test; {var temp1, temp2, temp3 : double;      cs1 : single;    ...    temp1 := 876600*3600 + 8640184.812866;      temp2 := 876600.0*3600.0 + 8640184.812866;      cs1 := 876600*3600;      temp3 := cs1 + 8640184.812866;      writeln(temp1,'  ',temp2,'  ',temp3);} var  cs1: single; begin  WriteLn(SizeOf(876600*3600 + 8640184.812866));  WriteLn(SizeOf(876600.0*3600.0 + 8640184.812866));  cs1 := 876600*3600;  WriteLn(SizeOf(cs1 + 8640184.812866))end.  The result, on an x86_64 system, is straight 10s. I'm not entirely sure what this is telling us, and whether it's in any way useful. But it would appear to suggest that the loss of precision is more subtle than a "least size" choice of types. MarkMLl ---------------------------------------------------------------------- tetrastes · 2022-11-19 21:20 It becomes funny... Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i program test; var  cs1: single; begin  WriteLn(SizeOf(8640184.1));  WriteLn(SizeOf(8640184.0));  WriteLn(SizeOf(3.1));  WriteLn(SizeOf(3.0));  cs1 := 876600*3600;  WriteLn(SizeOf(cs1 + 8640184.0));  readln;end.  results in Code: Bash  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i 1041044 ---------------------------------------------------------------------- PascalDragon · 2022-11-19 21:55 FPC uses the smallest floating point type that can represent a floating point number. This can influence calculations as you noticed (in addition to the floating point type with the highest precision differing between all x86 targets except x86_64-win64 and any non-x86 target plus x86_64-win64). If you want these constants to use the highest precision type you need to use {$ExcessPrecision On}. ---------------------------------------------------------------------- tetrastes · 2022-11-19 22:08 Quote from: PascalDragon on November 19, 2022, 10:55:01 pmFPC uses the smallest floating point type that can represent a floating point number. So why SizeOf(3.1) = 10 ? Quote from: PascalDragon on November 19, 2022, 10:55:01 pmIf you want these constants to use the highest precision type you need to use {$ExcessPrecision On}. Added {$ExcessPrecision On} to code in my previous post. Result hasn't change. ---------------------------------------------------------------------- PascalDragon · 2022-11-19 23:03 Quote from: tetrastes on November 19, 2022, 11:08:17 pmQuote from: PascalDragon on November 19, 2022, 10:55:01 pmFPC uses the smallest floating point type that can represent a floating point number. So why SizeOf(3.1) = 10 ? Because 3.1 can't be represented as a Single. If you add a Writeln(3.1) you'll see something like this in the assembly code: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i .section .rodata.n__$TFLOAT$_Ld1,"d"        .balign 16.globl  _$TFLOAT$_Ld1_$TFLOAT$_Ld1:# value: 0d+3.09999999999999999991E+0000      .byte   102,102,102,102,102,102,102,198,0,64  Quote from: tetrastes on November 19, 2022, 11:08:17 pmQuote from: PascalDragon on November 19, 2022, 10:55:01 pmIf you want these constants to use the highest precision type you need to use {$ExcessPrecision On}. Added {$ExcessPrecision On} to code in my previous post. Result hasn't change. Please try with FPC 3.3.1, there seems to be a bug in older versions. ---------------------------------------------------------------------- tetrastes · 2022-11-19 23:32 Quote from: PascalDragon on November 20, 2022, 12:03:23 amBecause 3.1 can't be represented as a Single. Of course it can. Only not exactly. But it can't be represented exactly as double or extended also, as you can see from your assembly code. The difference is only precision: Write(3.1) gives 3.09999999999999999991E+0000 (extended). Write(single(3.1)) gives  3.099999905E+00. I have another question: why 31.0 is represented as 3.100000000E+01, and not 3.09999999999999999991E+0001 ? ---------------------------------------------------------------------- uart · 2022-11-20 00:02 Quote from: MarkMLl on November 19, 2022, 09:34:26 pmCode: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i WriteLn(SizeOf(8876600.0*3600.0 + 8640184.812866));  The result, on an x86_64 system, is straight 10s Yes, the 8640184.812866 cannot be represented exactly as a single precision, so the expression gets promoted to a higher precision at that point in the calculation, but the 876600.*3600.0 part appears to be done at single precision (4 bytes). ---------------------------------------------------------------------- uart · 2022-11-20 00:13 Quote from: PascalDragon on November 19, 2022, 10:55:01 pmFPC uses the smallest floating point type that can represent a floating point number. Thanks for clarifying that PD, that's what I was suspecting. ---------------------------------------------------------------------- uart · 2022-11-20 00:15 Quote from: tetrastes on November 20, 2022, 12:32:39 amOf course it can. Only not exactly. Yes, lots of decimal numbers cannot be represent exactly with any sized mantissa, but they will be handled by whatever is the default precision. The claim is that IF (the if is important here) a literal can be represented EXACTLY in single precision, then that's what it will use. ---------------------------------------------------------------------- tetrastes · 2022-11-20 07:53 Quote from: uart on November 20, 2022, 01:15:15 amThe claim is that IF (the if is important here) a literal can be represented EXACTLY in single precision, then that's what it will use. Yes, but this is not the same as Quote from: PascalDragon on November 19, 2022, 10:55:01 pmFPC uses the smallest floating point type that can represent a floating point number. ---------------------------------------------------------------------- MarkMLl · 2022-11-20 10:41 Quote from: tetrastes on November 20, 2022, 08:53:25 amYes, but this is not the same as Cutting to the chase, all of which emphasises the importance of carefully-chosen algorithms. I'm concerned by the complexity of the code quoted in #2 Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i PROCEDURE JDAY          ( Year,Mon,Day,Hr,MIN                   : INTEGER;                             SEC                                   : EXTENDED;                               VAR JD                              : EXTENDED );   BEGIN     JD:= 367.0 * Year          - TRUNC( (7* (Year+TRUNC ( (Mon+9)/12) ) ) * 0.25 )          + TRUNC( 275*Mon / 9 )          + Day + 1721013.5          + ( (SEC/60.0 + MIN ) / 60.0 + Hr ) / 24.0;  { UT in days }          { - 0.5*SGN(100.0*Year + Mon - 190002.5) + 0.5; }   END;  { PROCEDURE JDAY }  when compared with this which is derived from Meeus via NOAA: Code: Pascal  [Select][+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i FUNCTION CalcJD(year, month, day: INTEGER): DOUBLE; //************************* **********************************************///* Name:    CalcJD                                                   *///* Type:    Function                                                   *///* Purpose: Julian day from calendar day                               *///* Arguments:                                                        *///*   year : 4 digit year                                             *///*   month: January = 1                                              *///*   day  : 1 - 31                                                   *///* Return value:                                                     *///*   The Julian day corresponding to the date                        *///* Note:                                                             *///*   Number is returned for start of day.  Fractional days should be   *///*   added later.                                                    *///********************************** *************************************/ VAR     A, B: INTEGER;        JD: DOUBLE; BEGIN  IF month <= 2 THEN BEGIN    year -= 1;    month += 12  END;  A:= Floor(year / 100);  B:= 2 - A + Floor(A / 4);  JD:= Floor(365.25 * (year + 4716)) + Floor(30.6001 * (month + 1)) + day + B - 1524.5;  RESULT:= JDEND { CalcJD } ;  The bottom line is that because floating point numbers cannot be relied upon to be precise, numerical methods have to be chosen carefully. MarkMLl ---------------------------------------------------------------------- Thaddy · 2022-11-20 10:56 I totally agree with your remark about choosing the right algorithm. But I would like to highlight that there is a difference between FPU - any - supported Floating points, for speed, and the more exact - but still arbitrary - precision you can often achieve with bigintegers and scaling at the cost of speed. ---------------------------------------------------------------------- MarkMLl · 2022-11-20 11:05 Absolutely, but this is not a case where speed is significant. Speed /might/ be significant if a large number of observations etc. were to be reduced as a batch, but at that point things like astronomical time and coordinates of the observer (both with multiple appropriate variants) would be precomputed. Can things like AVX-512 be used to improve precision (rather than throughput)? MarkMLl ---------------------------------------------------------------------- Thaddy · 2022-11-20 11:20 Yes, but it still stays arbitrary to infinity. That is the crux. Scaled bigintegers can grow to any precision above 512. But still stay approximations for some values that end in infinity, like E or PI just with a higher precision than 512 if needed. I know of some FPGA designs with extremely high precision using scaled bigintegers. (E.g. NASA has at least one public domain design for that. I also know ASML has at least one, but that is not public) ---------------------------------------------------------------------- MarkMLl · 2022-11-20 11:31 True. But OTOH, for calculations involving astronomy hence time we know the maximum achievable precision (currently around 10 microarcsecs for stellar positions), so have something finite to aspire to. MarkMLl ---------------------------------------------------------------------- Thaddy · 2022-11-20 11:39 Yes. the same goes for - like ASML does - for micro. See https://en.wikipedia.org/wiki/Extreme_ultraviolet_lithography which in its current implementation can't be expressed in 512 accurate enough. I may add that the 5 microns in the wiki for Samsung only means they are launch customers for ASML. My source about NASA is https://usa.nissannews.com/en-US/releases/dr- maarten-sierhuis-bio who is a very close friend of mine. If you are interested, contact him. It is perfectly OK to mention my name as reference as long as it is about computer science. He is always busy, but he always answers. ---------------------------------------------------------------------- MarkMLl · 2022-11-20 13:06 Well... 10 microarcsec is of the order of 10e-12 of the circle, so in principle IEEE 64-bit representation isn't too bad. PROVIDED that the chosen numerical method doesn't degrade this beyond four or five bits of precision, so IEEE 80-bit (as implemented on x86) is arguably OK PROVIDED that calculations are actually being done (rather than intermediates being merely stored) to this precision. Now I don't think Meeus is challenged as an authority, and I'm fairly confident in NOAA as an intermediate source (the link I posted earlier was my Pascal transcription of their Javascript, complete with their methodical comments). But I've mentioned before that when I've used this for Sunset calculations I've observed a seasonally-varying difference of roughly three minutes between that and a fairly well- regarded Android GPS (etc.) app: I never investigated whether this variation was smooth or was a step-change at various points of the year. Since the Android app was closed-source, I've got no idea of where they got their algorithms, but it's not inconceivable that this was down to an imprecise Julian date/time calculation. Going back to OP's position, GNU Octave is of course highly-regarded although my recollection is that this depends on the selection of backend computational software. And FPC can, by now, be considered to have fairly robust numerical performance. So this really boils down to how the required compilation is expressed, and the possibly mistaken assumptions that any compiler can make about the precision required of intermediate steps. Meeus, as a professional, was undoubtedly aware of this, and almost certainly benefited from feedback he received from people implementing his algorithms on disparate architectures with a much wider variety of compilers than we are used to seeing these days. The same does not necessarily apply to other people who attempt to reimplement his work from the ground up. MarkMLl ---------------------------------------------------------------------- PascalDragon · 2022-11-20 21:53 Quote from: tetrastes on November 20, 2022, 08:53:25 amQuote from: uart on November 20, 2022, 01:15:15 amThe claim is that IF (the if is important here) a literal can be represented EXACTLY in single precision, then that's what it will use. Yes, but this is not the same as Quote from: PascalDragon on November 19, 2022, 10:55:01 pmFPC uses the smallest floating point type that can represent a floating point number. I should have been a bit more precise FPC uses the smallest floating point type that can exactly represent a floating point number or the one with the highest precision if it can't be exactly represented by any.