Android Question Math problem (arithmetic on integers)[solved]

[Didier9]

I am working on a project where I need to convert back and forth between date in month-day-year and in Modified Julian Date formats.
I have conversion routines that I have used and verified in C under both Keil C51 compiler and gcc.
I use an on-line converter that has also been checked:
http://www.csgnetwork.com/julianmodifdateconv.html
The routines were obtained from reputable sources (see references in code comments)

I adapted the routines for B4A and B4J. The changes are only the way variables are declared. Otherwise it is plain vanilla arithmetic.

'=============================================================================
' These functions are adapted to B4A from Tom VanBaak's web site, leapsecond.com
' For background on Julian Day And Modified Julian Day, see http://tycho.usno.navy.mil/mjd.html
'=============================================================================

' Convert year/month/day calendar date to MJD (Modified Julian Day).
' - year is (4-digit calendar year), month is (1-12), day is (1-31)
' - valid For Gregorian dates from 17-Nov-1858 (adapted from sci.astro FAQ)
Sub ymd_to_mjd( sd As DateType ) As Long
    Dim year, month, day As Int, mjd As Long

    year = sd.y
    month = sd.m
    day = sd.d
    mjd = _
        367 * year _
        - 7 * (year + (month + 9) / 12) / 4 _
        - 3 * ((year + (month - 9) / 7) / 100 + 1) / 4 _
        + 275 * month / 9 _
        + day + 1721028 - 2400000
    Return( mjd )
End Sub ' ymd_to_mjd()

' Convert Modified Julian day (mjd) to year/month/day calendar date.
' - year is (4-digit calendar year), month is (1-12), day is (1-31)
' - adapted from Fliegel/van Flandern ACM 11/#10 p 657 Oct 1968
Sub mjd_to_ymd( mjd As Long ) As DateType
    Dim J, C, Y, M, d, mo, y As Long
    Dim sd As DateType

    J = mjd + 2400001 + 68569
    C = 4 * J / 146097
    J = J - (146097 * C + 3) / 4
    Y = 4000 * (J + 1) / 1461001
    J = J - 1461 * Y / 4 + 31
    M = 80 * J / 2447
    d = J - 2447 * M / 80
    J = M / 11
    sd.d = d
    mo = M + 2 - (12 * J)
    sd.m = mo
    y = 100 * (C - 49) + Y + J
    sd.y = y
    Return( sd )
End Sub ' mjd_to_ymd()

When I run the code under B4A or B4J, the calculations are off by a day or two both ways but in the same direction so they do not even cancel each other.
Here is the log:

date in: 4-15-2000
mjd is 51647
date out: 4-12-2000

According to the web calculator, MJD 51647 is April 13, 2000.

I have been working on this for a little while and can't seem to find what's wrong.
I suspect it is either an order of precedence issue or the way variable types may be handled, or maybe the way Java optimizes calculations that may result in overflow?
The gcc code uses the same variable types (int and long) and I believe they actually match (ints are signed 32 bit integers and longs are signed 64 bit integers).

Any pointer to potential differences between Java and C would be appreciated.

 

[agraham]

Where is the original code that you translated?

http://tycho.usno.navy.mil/mjd.html is not a valid web site.

 

[agraham]

If you are sure that those Subs are correctly translated from C then I suspect the problem is that those algorithms are not in fact exact mathematical transformations but may rely upon implicit rounding to Int and/or Long values in the intermediate calculations. Inspection of the the generated code (under the Objects\src folder) shows that B4X casts divide operations to operate as Doubles whose fractional values may explain the discrepancies. I've tried using Floor to round the intermediate results in ymd_to_mjd() and got closer but I'm still one day out but maybe I've not got it quite right

Also note that B4X is case-insensitive so in
Dim J, C, Y, M, d, mo, y As Long
Y and y are the same Long variable - though in this case it doesn't matter.

 

[mc73]

Agraham is correct. The original formula uses lots of ints. One way to achieve this is by casting your elements to int and then inserting in the final step. Another way is to use the built-in floor function as Agraham suggested.
Here you'll find the correct formula for Julian Date from which you can always obtain the modified one.
http://scienceworld.wolfram.com/astronomy/JulianDate.html

 

[Marcus Araujo]

Hello,

This would be a direct conversion from the javascript function in the site you referred:
http://www.csgnetwork.com/julianmodifdateconv.html

Sub mjd_to_ymd2(mjd As Int) As dateType

    ' Julian day
    Dim jd As Float = Floor (mjd) + 2400000.5

    ' Integer Julian day
    Dim jdi As Int = Floor (jd)
       
    ' Fractional part of day
    Dim jdf As Float = jd - jdi + 0.5
       
    ' Really the Next calendar day?
    If (jdf >= 1.0) Then
        jdf = jdf - 1.0
        jdi  = jdi + 1
    End If


    Dim l As Int = jdi + 68569
    Dim n As Int = Floor (4 * l / 146097)
       
    l = Floor (l) - Floor ((146097 * n + 3) / 4)
    Dim year As Int = Floor (4000 * (l + 1) / 1461001)
       
    l = l - (Floor (1461 * year / 4)) + 31
    Dim month As Int = Floor (80 * l / 2447)
       
    Dim day As Int = l - Floor (2447 * month / 80)
       
    l = Floor (month / 11)
       
    month = Floor (month + 2 - 12 * l)
    year = Floor (100 * (n - 49) + year + l)
       
    Dim sd As dateType
    sd.d = day
    sd.m = month
    sd.y = year
    Return sd
End Sub

 

[Didier9]

If you are sure that those Subs are correctly translated from C

Here is the gcc code:

struct sdate {
  uint16_t y;
  uint8_t m;
  uint8_t d;
};

// Convert year/month/day calendar date to MJD (Modified Julian Day).
// - year is (4-digit calendar year), month is (1-12), day is (1-31)
// - valid for Gregorian dates from 17-Nov-1858 (adapted from sci.astro FAQ)
uint32_t ymd_to_mjd( struct sdate sd ){
    int32_t year, month, day, mjd;

    year = (int32_t)sd.y;
    month = (int32_t)sd.m;
    day = (int32_t)sd.d;
    mjd =
        367 * year
        - 7 * (year + (month + 9) / 12) / 4
        - 3 * ((year + (month - 9) / 7) / 100 + 1) / 4
        + 275 * month / 9
        + day + 1721028 - 2400000;

    return( (uint32_t)mjd );

} // ymd_to_mjd()

// Convert Modified Julian Day (MJD) to year/month/day calendar date.
// - year is (4-digit calendar year), month is (1-12), day is (1-31)
// - adapted from Fliegel/van Flandern ACM 11/#10 p 657 Oct 1968
struct sdate mjd_to_ymd( uint32_t mjd ){
    int32_t J, C, Y, M;
    struct sdate sd;

    J = mjd + 2400001 + 68569;
    C = 4 * J / 146097;
    J = J - (146097 * C + 3) / 4;
    Y = 4000 * (J + 1) / 1461001;
    J = J - 1461 * Y / 4 + 31;
    M = 80 * J / 2447;
    sd.d = (uint8_t)(J - 2447 * M / 80);
    J = M / 11;
    sd.m = (uint8_t)(M + 2 - (12 * J));
    sd.y = (uint)(100 * (C - 49) + Y + J);
    return( sd );

} // mjd_to_ymd()

I suspected something about the way Java treats complex formulas since it is so aggressive at optimizing. The original code may have taken advantage of peculiar C idiosyncrasies.

Then, if Java is converting the integer variables to double, that obviously won't help.

Let me work on that a little bit, I thought about using floor(), I need to look back exactly at what it does with negative numbers.

It is a little bit frustrating to have to do that and no be able to just do integer math...

 

[Didier9]

Here is a short B4J project that shows it

 

[Didier9]

I fixed the first routine by adding floor() for every division.

Sub ymd_to_mjd( sd As DateType ) As Long
    Dim year, month, day As Int, mjd As Long

    year = sd.y
    month = sd.m
    day = sd.d
    mjd = _
        367 * year _
        - Floor(7 * (year + Floor((month + 9) / 12)) / 4) _
        - Floor(3 * (Floor((year + Floor((month - 9) / 7)) / 100) + 1) / 4) _
        + Floor(275 * month / 9) _
        + day + 1721028 - 2400000
    Return( mjd )
End Sub ' ymd_to_mjd()

Now that returns 51649 for April 15, 2000 and it matches the web calculator!
So it seems to confirm the calculations are made on doubles.

 

[Didier9]

Did this to the other routine and it works too:

Sub mjd_to_ymd( mjd As Long ) As DateType
    Dim J, C, Y, M, d, mo, y As Long
    Dim sd As DateType

    J = mjd + 2400001 + 68569
    C = Floor(4 * J / 146097)
    J = J - Floor((146097 * C + 3) / 4)
    Y = Floor(4000 * (J + 1) / 1461001)
    J = J - Floor(1461 * Y / 4) + 31
    M = Floor(80 * J / 2447)
    d = J - Floor(2447 * M / 80)
    J = Floor(M / 11)
    sd.d = d
    mo = M + 2 - (12 * J)
    sd.m = mo
    y = 100 * (C - 49) + Y + J
    sd.y = y
    Return( sd )
End Sub ' mjd_to_ymd()

Problem fixed, thanks a lot!

 

[Didier9]

Inspection of the the generated code (under the Objects\src folder) shows that B4X casts divide operations to operate as Doubles whose fractional values may explain the discrepancies.

This was the problem right there, you nailed it on the head! Thanks again.

 

[Didier9]

Also note that B4X is case-insensitive so in
Dim J, C, Y, M, d, mo, y As Long
Y and y are the same Long variable - though in this case it doesn't matter.

Yes, I did it for clarity, or at least that's what I thought at the time. A while later, it's not very clear to me why I did it...

 

[Didier9]

For background, MJD is a very convenient format to calculate dates separated by a certain number of days/weeks.
In this application, I use it to correct for the GPS system week rollover issue which causes dates to be off by 1024 weeks, or 7168 days.
The same problem occurs again every 19.6 years.