Notes on Julian and Gregorian calendar algorithms

Notes_on_Julian_and_Gregorian_calendar_algorithms

table.speeds td.fig {
text-align: right;
}

table.speeds_null_null td.valign {
text-valign: top;
}

Notes on Julian and Gregorian calendar algorithms

Contents

Foreword

These notes aim to persuade people using floating point algorithms for working
with dates in Julian and Gregorian calendars to use integer algorithms instead.
In my experience the integer algorithms are substantially faster,
and – in my reasoned opinion – they are easier to understand,
making them easier to maintain and to adapt for other uses.

I have tried to avoid errors in these notes: if anyone (correctly) points out
an error I will make a dated correction with a full credit to the person
reporting the error, adding a footnote showing the text before the correction.

If anything here is unclear, please comment:
I don’t guarantee to make changes because these notes intentionally assume
a detailed background knowledge (see Preliminaries),
I will at least think about making things clearer.

If anyone has a clearer way of explaining something, and agrees to have their
explanation posted here, I am more than willing to put their explanation here
with a full credit to them.

Caveat: while I have tried to avoid errors and have tested the code (see the
code examples below), I make no warranties that the algorithms or code or psuedo-code below
are correct, and anyone who wishes to use any of the algorithms or code or pseudo-code below
(unmodified or modified) should make their own judgment
on the correctness and suitability of the algorithms or code or pseudo-code.

I have derived the integer algorithms (including some faster algorithms, and
some shortcuts for dates in the range years 1901 to 2099) from first principles,
but I believe they are unlikely to be original: once you understand the basic
cycles of Julian and Gregorian calendars (1461 days, 146097 days, 153 days),
and the periods of days that make up those basic cycles, then the calculations
needed are fairly obvious, although you do need to be careful about the details.

(Moreover, I may well have been subconciously influenced by some aspects of code
in the Java JodaTime and ThreeTen projects, and by the Ruby Date class. That said,
the latter currently uses floating point algorithms, and I have not looked
in detail at the Java projects equivalents of the algorithms discussed here.)

So if anyone is aware of good alternative explanations of these algorithms, or,
indeed, of better implementations of them, I will be very happy to link to
(or cite references for) those explanations or implementations.

My ultimate aim is that there should be an online resource pulling together
calendar algorithms (not just for Julian and Gregorian calendars) which gives
a balanced disussion of the pros and cons of different algorithms, so that
projects like Ruby, Python, etc (but not only open source projects) can make
an informed choice about the algorithms they use.

(This aim is prompted by my looking in detail at implementations of floorDiv
and similar arithmetic functions in two major projects and in Apache commons,
and either finding very obscure bugs or – in my view – finding the code could be
improved, either in clarity or speed or both. And by finding similar bugs
in my first attempts to implement floorDiv and similar functions.)

The idea is that although date calculations and time calculations are connected,
I believe it should be possible to do date and time calculations by:

  • Date calculations work with Chronological Julian Day Numbers (CJDN) to deal with:
    • Julian and Gregorian calendar years and months and days;
    • commercial years and weeks and days;
    • other calendars.
  • Time calculations deal with:
    • hours, minutes, seconds (and possibly with leap seconds);
    • time zones;
    • Daylight Saving Time (DST) switches;
    • adding or subtracting a number of whole days to or from the time,

      including what happens if adding or subtracting the days
      crosses a DST switch:

      crossing a DST switch is the difficult part of dealing with adding
      or subtracting days;

    • dates (years, months, etc) by:
      • (a) calculate the CJDN of the time;
      • (b) use a date calculation to calculate a new CJDN from that CJDN;
      • (c) calculate the difference in days (b) – (a);
      • (d) use the difference in days to add or subtract a number of whole days
        to or from the time to get a new time.

In short, if a time system can deal with numbers of whole days and can provide a CJDN,

I believe we can use CJDNs to make date and calendar calculations,

and it should be reasonably easy to use any calendar with any time system.

At least, that’s the hope: there may be annoying “minor details”
which make it impractical.

Preliminaries

These notes assume familiarity with calculations for working with
Julian and Gregorian calendars. See, for example:

They also assume familiarity with limitations of computing using
integers and floating point numbers, including possibilities of:

  • Integer calculations resulting in numbers which are too big
    to be held in the type of integer being used, that is “overflow”.

  • Floating point calculations resulting in numbers which are too big
    to be held exactly in the type of floating point being used:
    this makes calendar algorithms using floating point calculations
    unreliable for very large dates.

    This is is unlikely to be a practical limitation on floating point
    algorithms: my point is that if someone objects to integer algorithms
    because there are limits on the range of dates that a particular
    integer type can handle, I can observe that there is a similar objection
    to floating point algorithms, and note that both these objections
    only apply to dates which are far outside the range of dates
    that any practical application will use.
    So limitations on the range of dates for an integer algorithm
    is not a valid reason for using a floating point algorithm instead.

  • Floating point calculations using numbers which cannot be held exactly
    as floating point numbers, for example the fraction 306/10, that is 30.6.
    When such numbers are used in date calculations we must be very careful,
    and we may need to “adjust” such numbers from the theoretically correct
    values to “fudged” values which correct for possible floating point
    rounding errors affecting the correct working of an algorithm.
    This “fudging” can obscure what an algorithm is really doing.

Reasons to change from floating point to integer algorithms

  • Theoretically correct algorithms for working with dates in
    Julian and Gregorian calendars use integer arithmetic:

    • to see this look at the articles linked above.
  • But floating point algorithms are a mixture of:
    • direct translations of integer arithmetic to floating point arithmetic;
    • adaptations of integer arithmetic to floating point arithmetic
      with some numbers being changed from theoretically correct values
      to “fudged” values so the floating point calculations give correct results;
      this is because there are some rational numbers which cannot be represented
      exactly by floating point numbers, for example 306/10 == 30.6.
  • Integer arithmetic algorithms are somewhat easier to understand than
    floating point algorithms (which have some changes to the theoretically
    correct numbers to avoid floating point rounding errors),
    from which it follows that integer arithmetic algorithms will be
    easier to maintain and to adapt than floating point algorithms.

  • 32bit signed integers give a range of dates more than sufficient
    for all practical calculations for Julian and Gregorian calendars.

  • For typical computers these days 32bit signed integer arithmetic
    seems to be usually faster than 64bit floating point arithmetic.

The signed integer arithmetic algorithms do have a few more lines of code
than the floating point algorithms, but in my opinion that is a minor
disadvantage compared with the above major advantages.

So I consider that practical implementations of algorithms for working with
Julian and Gregorian calendars should use 32bit signed integer arithmetic
unless it can be shown that for the application in question
floating point arithmetic is likely to be faster.

Abbreviations and definitions

  • CJDN is Chronological Julian Day Number;
  • floor(x) means round x to nearest integer less than or equal to x,

    in other words round towards “minus infinity”;

  • floorDiv(a, b) means divide a by b; if result is not an exact integer
    round result to nearest integer less than or equal to result,

    in other words round result towards “minus infinity”;

  • floorMod(a, b) is the value such that: a == b * floorDiv(a, b) + floorMod(a, b);

    in other words: floorMod(a, b) = a – b * floorDiv(a, b);

  • truncateDiv(a, b) means divide a by b; if result is not an exact integer
    round result to nearest integer between “a” and zero,

    in other words round result towards zero.

Limitation of integer and floating point algorithms

  • In implementations of integer arithmetic algorithms we must take care
    to avoid integer overflows.

    Using 32bit signed integers we must
    be careful when (very roughly):

    • CJDN 530,000,000,
      where 530,000,000 is roughly (2**31) / 4.0

    • or year 1,400,000,
      where 1,400,000 is roughly 530,000,000 / 365.25.
  • In implementations of floating point algorithms there are likely to be
    limitations on the range of integers which can be stored accurately
    in floating point numbers. For example, Meeus’s floating point algorithm
    will be unreliable if any of the calculations in the algorithm result in
    numbers outside that range: using 64bit floating point (Java or C double)
    the number 9,007,199,254,740,995 cannot be held exactly as floating point,
    and we must be careful when (very roughly):

    • CJDN 2,200,000,000,000,000,
    • or year 6,100,000,000,000.

So whether we use integer or floating point arithmetic for an algorithm,
we should check the CJDN is in a range for which the arithmetic works.

But for all practical calculations with Julian and Gregorian calendars
it is very very unlikely that the CJDN will be outside that “working” range:

using 32bit signed integers we have a range of years of
(very roughly) -5,800,000 to 5,800,000,

the calculation being
(-2,147,483,648 / 365.25 – 4712) to (2,147,483,648 / 365.25 – 4712),

where 2,147,483,648 is 2 to the power 31, that is 2**31.

So if we work with CJDNs being stored in 32bit signed integers:

  • We have a range of years of about -5,800,000 to 5,800,000, more than
    sufficient for practical calculations with Julian and Gregorian calendars.

  • If we use integer algorithms, we should check for the (very very unlikely)
    possibility that integer overflow may occur. But such checks are fairly
    easy to make, and in my experience 32bit signed integer algorithms using
    such checks are still much faster than 64bit floating point algorithms.

Speeds of integer and floating point algorithms

Some benchmarking in Java suggests that 32bit signed integer arithmetic is
substantially faster than 64bit floating point arithmetic.
So for practical calculations with Julian and Gregorian calendars
it is likely to be faster to use 32bit signed integer arithmetic algorithms
than 64bit floating point (Java or C double) arithmetic algorithms.

The table below has some benchmarking which show the substantial
increases in speed by avoiding 64bit floating point arithmetic:

  • if we use ordinary 32bit signed integer arithmetic,
    the algorithms take about 50% of the time for floating point,
    a saving of time taken of about 50%;

  • if we also use some faster special algorithms,
    the algorithms take about 30% of the time for floating point,
    a saving of time taken of about 70%;

  • if we also use these faster algorithms, and furthermore
    use shortcuts for Gregorian dates in the range year 1901 to year 2099,

    the algorithms take about 20% of the time for floating point,
    a saving of time taken of about 80%.

I believe that for most current real world applications the most likely range
of dates is Gregorian years 1901 to 2099, so for most real world
applications if we adopt all of the above time-saving measures we can
save a large amount of processing time for Gregorian calendar calculations.

The overall savings in time will naturally be much less because typically
the calendar calculations will be a relatively small part of the overall
processing.

But, as I have argued above, using the integer arithmetic algorithms is faster
and – in my view – the code of the integer algorithms is more understandable.
So I consider that using integer algorithms is a “win-win”.

Table comparing the speeds of the algorithms

  • Benchmark tests were run for two computers, a Toshiba laptop and a Samsung netbook.

  • For each computer a benchmark test was 4 settle down runs, then 12 timed runs.

  • For each computer two benchmark tests were made: the timings were very similar,
    as can be seen in the first two columns for each computer.

  • The standard deviations of the averages of the times for each run in each benchmark test
    were all less than 3%, and most standard deviations were less than 1%.

  • The benchmark tests were for ranges of dates that are most likely to be used
    in date calculations, that is 1901 to 2099.

  • For CJDN to YMD runs the CJDN range was 2415386 to 2488069, representing a range
    in years of about 1901 to 2099,

    with the range repeated 10 times,
    for a total of about 726,000 iterations in each timed run.

  • For YMD to CJDN runs the year range was 1901 to 2099, the month range 1 to 12, and
    the day range 1 to 28, so that all the dates would be valid dates.

    These ranges were repeated 20 times, for a total of about 1,337,000 iterations in each timed run.

  • The “null” runs are so we can exclude the effect of benchmarking overheads,
    making it easier to compare the speeds of the algorithm implementations.

    The “null” runs do not make any calculations: they just do the iterations and return a constant value.

Columns in the table:

  • two runs: Timings for two benchmark runs – the figures are the average nano-seconds per iteration.

  • aa: The “adjusted average” timings are the average of the averages of the two runs
    minus the average of the relevant null runs.

  • aas: These are the “adjusted average” timings “scaled” so the longest timing is 100.

    • The Samsung netbook has a lower specification than the Toshiba laptop,
      and we expect it to take longer.

      So we “scale” to make it easier to compare the Toshiba laptop and Samsung netbook timings.
Laptop Netbook
Algorithm     two runs;   aa, aas;     two runs;   aa, aas;
cjdn2ymd_meeus_float 267, 268; 247, 100; 605, 605; 561, 100;
cjdn2ymd_meeus_float_shortcuts 188, 186; 166, 67; 422, 419; 376, 67;
cjdn2ymd_meeus_integer 131, 131; 110, 45; 343, 344; 299, 53;
cjdn2ymd_integer_etc 69, 69; 48, 19; 189, 190; 145, 26;
cjdn2ymd_integer_shortcuts_etc 51, 51; 30, 12; 135, 135; 91, 16;
cjdn2ymd_null 21, 21; 0, 0; 44, 45; 0, 0;
cjdn2ymd_null_null 3, 3; -18, -7; 10, 10; -35, -6;
ymd2cjdn_float 146, 146; 141, 100; 373, 374; 362, 100;
ymd2cjdn_float_shortcuts 74, 74; 69, 49; 186, 194; 178, 49; (*r)
ymd2cjdn_integer 58, 58; 53, 38; 238, 238; 226, 63; (*r)
ymd2cjdn_integer_etc; 29, 28; 24, 17; 102, 103; 91, 25;
ymd2cjdn_integer_shortcuts_etc 20, 20; 15, 11; 52, 52; 40, 11;
ymd2cjdn_null 5, 5; 0, 0; 12, 12; 0, 0;

(*r): Note that for these two timings the Samsung relative speeds are the reverse of the Toshiba.

These CJDN to YMD algorithms return a 3 element array of 32bit signed integers:

cjdn2ymd_meeus_float: Meeus’s floating point algorithm;
cjdn2ymd_meeus_float_shortcuts: Meeus’s floating point algorithm, but with a
shortcut for CJDNs of years 1901 to 2099;
cjdn2ymd_meeus_integer integer version of Meeus’s float algorithm;
cjdn2ymd_meeus_integer_etc: integer algorithm using some faster algorithms;
cjdn2ymd_integer_shortcuts_etc: as cjdn2ymd_meeus_integer_etc, but with a
shortcut for CJDNs of years 1901 to 2099;
cjdn2ymd_null: “null”.

This CJDN to YMD “null” algorithm returns a 32bit signed integer:

cjdn2ymd_null_null: “null” but returns a 32bit signed integer;

it seems faster to return a “packed” 64bit signed integer than an integer array,

so the (aa) and (aas) timings above exclude creating an array object.

But if anyone wishes to adjust the table to include creating an array object,

then a rough adjustment is to use cjdn2ymd_null_null as the “base” null.

These YMD to CJDN algorithms return a 32bit signed integer:

ymd2cjdn_float: a floating point algorithm;
ymd2cjdn_float_shortcuts: as ymd2cjdn_float, but with a shortcut for years 1901 to 2099;
ymd2cjdn_integer: integer algorithm;
ymd2cjdn_integer_etc: integer algorithm using some faster algorithms;
ymd2cjdn_integer_shortcuts_etc: as ymd2cjdn_integer_etc, but with a shortcut for years 1901 to 2099;
ymd2cjdn_null: “null”.

Explanation of code using Meeus’s floating point algorithm

For each step we first show code for the underlying integer arithmetic,
then (indented) we derive code which uses Meeus’s floating point algorithm.

The idea is use the variable “d” to hold days, starting with d as the CJDN.

  • Steps 1 to 4: from d calculate the year, and reduce d to exclude the year;
    later we may need to adjust the year by 1.

  • Steps 5 and 6: from new d calculate the month in year, and reduce d to exclude the month;
    later we may need to adjust the month, depending on the algorithm.

  • Now the new d is the day in the month, but we may need to adjust
    the day by 1, depending on the algorithm.

  • Step 7: we may need to adjust the month and year.
  • 1: d = CJDN of date;
    • if date is Gregorian change d to CJDN which in the Julian calendar
      has the same year, month, day as the Gregorian date;

      1867217 is CJDN of Gregorian 0400-03-01;

      • x = floorDiv((CJDN – 1867217) * 4 + 3, 146097);

        • or: x = floor((CJDN – 1867217 + 3/4.0) / 36524.25);

          so Meeus: x = floor((CJDN – 1867216.25) / 36524.25);
      • d = CJDN + 1 + x – floorDiv(x, 4);
        • or Meeus: d = CJDN + 1 + x – floor(x / 4.0);

  • 2: rebase from Julian -4712-01-01 to -4717-10-30 == 123 days before -4716-03-01;
    • d = d + 1524;

  • 3: calculate an offsetted year; (d – 123) rebases d to -4716-03-01;
    • y = floorDiv((d – 123) * 4 + 3, 1461);
      • or: y = floor((d – 123 + 3/4.0) / 365.25);

        in fact “+f” works for any f such that 3/4.0 <= f < 4/4.0;

      • so: y = floor((d – 122.25) / 365.25);

        in fact “-n” works for any n such that (just over 122) < n <= 122.25,

      • so Meeus: y = floor((d – 122.1) / 365.25);

        Using 122.1 gives a safety margin from both 122.0 and 122.25.

        (Or use 122.125 which can be held exactly as a floating point number.)

  • 4: change d so: 123 <= d < 489 == 123 + 366;
    • d = d – floorDiv(y * 1461, 4);
      • or: d = d – floor((y * 1461.0) / 4.0);

        so Meeus: d = d – floor(y * 365.25);

  • 5: calculate an offsetted month so: 4 <= m < 16 == 4 + 12;
    • m = floorDiv(d * 5 – 1, 153);
      • or: m = floor((d * 5 – 1) / 153.0);

        in fact “-f” works if: (just over 0) < f <= 1;

      • so: m = floor((d * 10 – 2) / 306.0);

        in fact “-f” works if: (just over 0) < f <= 2;

      • so: m = floor((d – 0.1) / 30.6);

        in fact “-f” works if: (just over 0) < f < (just less than 0.2);

      • so Meeus: m = floor(d / 30.6001),

        which is roughly: floor((d – small positive fraction) / 30.6);

  • 6: change d to the correct day in month so: 1 <= d <= 31;
    • d = d – floorDiv(m * 153 + 0, 5);
      • or: d = d – floor((m * 153 + 0) / 5.0);

        in fact “+f” works if: 0 <= f < (just less than 1);

      • so: d = d – floor((m * 306 + 0) / 10.0);

        in fact “+f” works if: 0 <= f < (just less than 2);

      • so: d = d – floor(m * 30.6 + 0.1);

        in fact “+f” works if: 0 <= f < (just less than 0.2);

      • so Meeus: d = d – floor(m * 30.6001);

        which is roughly: d – floor(m * 30.6 + small positive fraction);

  • 7: change m and y to the correct month and year;
    • if m >= 14 then
      • m = m – 13;
      • y = y – 4715;
    • else
      • m = m – 1;
      • y = y – 4716;

Now y, m, d are the year, month, day of the CJDN.

Note that for steps 2 to 7 we can use a different integer algorithm
which is a bit clearer than the integer version of Meeus.

But compared with that alternative, using the integer version of Meeus
gives an average saving of 14/12 integer additions per calculation.

So we use the integer version of Meeus, partly to make the saving,
but mostly to recognise that Meeus really thought about the algorithm,

and it would be a pity not to perpetuate that.

(Until I made a *proper* analysis I thought that the alternative
had on average slightly fewer integer additions per calculation.)

Implementations in Java of floating point and integer algorithms

Below are integer and floating point implementations of:

  • CJDN to year, month, day;
  • year, month, day to CJDN.

The Java implementations of CJDN to YMD return an integer array of year, month, day.

  • Faster Java code would return year, month, day “packed” into a 64bit integer.
  • Faster implementations in C would return year, month, day directly by using
    void functions which have year, month, day as “*arguments” to the function.

In this implementation of the CJDN to YMD integer algorithm the code is rearranged
to minimise the number of operations.

Faster algorithms and shortcuts

These algorithms avoid some divisions. This is because in Java (and in C?)
integer division is a relatively slow operation compared with addition, subtraction,
comparison, and – to a lesser extent – multiplication.

The divisions are avoided by:

  • in steps 1 and 4 use “bit shifts” to divide by 4;
  • in step 5 use “if statements” for binary search to get month in year from day in year;
  • in step 6 use an array to get days in year before a month in year;
  • in step 1 use a shortcut to convert the most likely Gregorian CJDNs to Julian CJDNs.

Minor purists may omit the shortcuts, but the code will usually be quite a bit slower.

Hardcore purists wanting minimal code will not use any of the faster algorithms or shortcuts,
but the minimal code will usually be much slower.

Leap Years in a mixed Julian and Gregorian calendar

I include this because I am aware of two projects which use this algorithm
for deciding on whether a date is in a leap year:

  • if the date is before the switch from the Julian to the Gregorian calendar:

    the year of the date is only a leap year if it is divisible by 4;

  • if the date is on or after the switch from the Julian to the Gregorian calendar:

    the year of the date is only a leap year if it is divisible by 4

    *and* (it is not divisible by 100 or it is divisible by 400).

In practice this works, but I believe it is theoretically wrong.

Suppose we are unwise and choose the last Julian date as 1900-10-10 with the first Gregorian date being 1900-10-24.

Then using the above algorithm the date 1900-09-09 is in a leap year,
but the date 1900-11-11 is not in a leap year.

It seems contradictory that whether a year is a leap year or not depends
on which date in the year we use.

Moreover, there is a 1900-02-29 in this mixed calendar, so it seems contradictory
to say 1900 is not a leap year.

So I suggest this algorithm which returns true for that “unwise” choice,
and which works for other awkward cases:

  • if the year of the date is not divisible by 4 it is not a leap year;
  • if the year of the date is divisible by 4
    and (is not divisible by 100 or is divisible by 400) it is a leap year;

  • this leaves the case of the year of the date is divisible by 100
    but is not divisible by 400:

    • if the calendar is an all Julian calendar the year of the date is a leap year;
    • if the calendar is an all Gregorian calendar the year of the date is not a leap year;
    • if the calendar is a mixed Julian and Gregorian calendar
      the year of the date is a leap year

      *unless* the year of the date has a February.28 and that February.28 is
      in the Gregorian part of the calendar,

      in which case the year of the date is *not* a leap year.

Here is pseudo-Java (and C) code for that algorithm (actual Java code is
here):


    public static double Infinity = 1.0 / 0.0;
    public static double MinusInfinity = -1.0 / 0.0;
    public static boolean isLeapYear(int yearv, double firstGregorianCJDNv) {
      if (isNotDivisibleBy4(yearv)) return false;
      if (1901 <= yearv && yearv = 0) {
        iy = yearv / 100;
      } else if (yearv > -100) {
        iy = 0;
      } else {
        iy = -((-(yearv + 100)) / 100) - 1; // avoid integer overflow
      }
      if (iy * 100 != yearv || isDivisibleBy4(iy)) return true;
      // if here yearv is divisible by 100 but not by 400
      if (firstGregorianCJDNv == Infinity) return true; // all Julian calendar
      if (firstGregorianCJDNv == MinusInfinity) return false; // all Gregorian calendar
      int iiy; // or instead of next lines use: int iiy = floorDivBy4(iy);
      if (iy >= 0) {
        iiy = iy / 4;
      } else if (iy > -4) {
        iiy = 0;
      } else {
        iiy = -((-(iy + 4)) / 4) - 1; // avoid integer overflow
      }
      // Next line says is leap year if firstGregorianCJDNv >= CJDN of yearv.03.01.
      // In next line cast to 64bit integer (long) to avoid integer overflow.
      // We can combine (1721120 + iiy) as 32bit integer which will not overflow.
      return firstGregorianCJDNv >= (1721120 + iiy) + ((long) iy) * 36524;
    }

Java code for algorithms

Caveat: while I have tried to avoid errors and have tested the code (see the
code examples below) I make no warranties that code is correct, and anyone who
wishes to use any of the code below (unmodified or modified) should make their
own judgment on the correctness and suitability of the code.

Although this code is in Java it is intended to be easily convertible to C.
Code for Java only would be a little more concise, and possibly a bit faster.

The code below tests the algorithms – the contents of the code are these methods:

  • main method, which has several calls to the method testMeeusEtc;
  • testMeeusEtc, which tests algorithms CJDN to YMD, and YMD to CJDN;
  • the basic integer and floating point algorithms:
    • cjdnToYMDuseMeeusInteger;
    • cjdnToYMDuseMeeusFloat;
    • ymdToCJDNinteger;
    • ymdToCJDNfloat;
  • faster integer and floating point algorithms:
    • cjdnToYMDintegerEtc;
    • cjdnToYMDintegerShortcutsEtc;
    • cjdnToYMDuseMeeusFloatShortcuts;
    • ymdToCJDNintegerEtc;
    • ymdToCJDNintegerShortcutsEtc;
    • ymdToCJDNfloatShortcuts.
  • utility integer date array and methods:
    • DaysInYearBeforeMonthBaseMonth3ArrayBeware – array;
    • dayInYearFirstDay0Base0301ToMonthBeware – faster way of getting month in year;
    • isLeapYear – checks whether a year is a leap year or not;
  • utility integer methods:
    • floorDivBeware – read the comments for why “Beware”;
    • floorModUseFloorDivResult – a quick way of calculating floorMod;
    • isDivisibleBy4(a) – a quick way of checking if a is divisible by 4;
    • isNotDivisibleBy4(a) – a quick way of checking if a is not divisible by 4;
    • floorDivBy4(a) – a faster floorDiv(a, 4);
    • floorModBy4(a) – a faster floorMod(a, 4);
    • floorModBy4UseFloorDivResult(a) – also a faster floorDiv(a, 4);

// package net.suliram.java.datejg;  // packaging may be useful
public class TestDateCalcs {

    // This has code which tests integer and floating point implementations of
    // algorithms for working with dates in Julian and Gregorian calendars.

    public static void main(String[] args) {
        double firstGregorianAllJulian    =  1.0 / 0.0;
        double firstGregorianAllGregorian = -1.0 / 0.0;
        testMeeusEtc(firstGregorianAllJulian, -146100 * 5, 146100 * 5);
        testMeeusEtc(firstGregorianAllJulian,
                              2415080 - 146100 * 5, 2488143 + 146100 * 5);
        testMeeusEtc(firstGregorianAllJulian,
                              500100100 - 146100 * 5, 500100100 + 146100 * 5);
        testMeeusEtc(firstGregorianAllGregorian, -146100 * 5, 146100 * 5);
        testMeeusEtc(firstGregorianAllGregorian,
                              2415080 - 146100 * 5, 2488143 + 146100 * 5);
        testMeeusEtc(firstGregorianAllGregorian,
                              500100100 - 146100 * 5, 500100100 + 146100 * 5);
    }

    public static void testMeeusEtc(double firstGregorianCJDNv,
                                    int minCJDNv, int maxCJDNv) {
        System.out.println("");
        System.out.println("*firstGregorian= " + firstGregorianCJDNv +
                           ", CJDN " + minCJDNv + " to " + maxCJDNv);
        int exCJDNv = minCJDNv;  // avoid overflow in for;
        int numCJDN = 0;
        int numNotOKtoYMD = 0;
        int numNotOKtoCJDNf = 0;
        int numNotOKtoCJDNi = 0;
        int numNotOKtoYMDsc = 0;
        int numNotOKtoCJDNfsc = 0;
        int numNotOKtoCJDNisc = 0;
        for (int cjdnv = minCJDNv; cjdnv = exCJDNv; cjdnv++) {
            numCJDN = numCJDN + 1;
            int[] ymdi = cjdnToYMDuseMeeusInteger(cjdnv, firstGregorianCJDNv);
            int[] ymdf = cjdnToYMDuseMeeusFloat(cjdnv, firstGregorianCJDNv);
            if (ymdi[0] != ymdf[0] || ymdi[1] != ymdf[1] || ymdi[2] != ymdf[2]) {
                numNotOKtoYMD = numNotOKtoYMD + 1;
            }
            int cjdni = ymdToCJDNinteger(ymdi[0], ymdi[1], ymdi[2],
                                         firstGregorianCJDNv);
            if (cjdni != cjdnv) {
                numNotOKtoCJDNi = numNotOKtoCJDNi + 1;
            }
            int cjdnf = ymdToCJDNfloat(ymdf[0], ymdf[1], ymdf[2],
                                       firstGregorianCJDNv);
            if (cjdnf != cjdnv) {
                numNotOKtoCJDNf = numNotOKtoCJDNf + 1;
            }
            int[] ymdisc = cjdnToYMDintegerShortcutsEtc(cjdnv, firstGregorianCJDNv);
            int[] ymdfsc = cjdnToYMDuseMeeusFloatShortcuts(cjdnv, firstGregorianCJDNv);
            if (ymdisc[0] != ymdi[0] || ymdisc[1] != ymdi[1] || ymdisc[2] != ymdi[2] ||
                ymdfsc[0] != ymdi[0] || ymdfsc[1] != ymdi[1] || ymdfsc[2] != ymdi[2]) {
                numNotOKtoYMDsc = numNotOKtoYMDsc + 1;
            }
            int cjdnisc = ymdToCJDNintegerShortcutsEtc(ymdi[0], ymdi[1], ymdi[2],
                                                       firstGregorianCJDNv);
            if (cjdnisc != cjdnv) {
                numNotOKtoCJDNisc = numNotOKtoCJDNisc + 1;
            }
            int cjdnfsc = ymdToCJDNfloatShortcuts(ymdf[0], ymdf[1], ymdf[2],
                                                  firstGregorianCJDNv);
            if (cjdnfsc != cjdnv) {
                numNotOKtoCJDNfsc = numNotOKtoCJDNfsc + 1;
            }
        }
        System.out.println("**count=" + numCJDN + "; NotOK:" +
                           " toYMD=" + numNotOKtoYMD +
                           ", toCJDNi=" + numNotOKtoCJDNi +
                           ", toCJDNf=" + numNotOKtoCJDNf +
                           "; shortcuts: toYMD=" + numNotOKtoYMDsc +
                           ", toCJDNi=" + numNotOKtoCJDNisc +
                           ", toCJDNf=" + numNotOKtoCJDNfsc + ";");
    }

    public static int[] cjdnToYMDuseMeeusInteger(int cjdnv, double firstGregorianCJDNv) {
        // *** uses integer Meeus ***
        int y, m, d;
        y = -4716;
        d = cjdnv;
        // This code avoids integer overflow in (m * 1461) and (d * 4),
        // and ensures that (d * 4 - 489) >= 0.
        // The constants 500100100 and 5000 have margins to allow small changes
        // to this algorithm without causing integer overflow problems.
        if (d  500100100 || d  500100100 || d < 5000) {
                m = floorDivBeware(d, 146097); // here "m" is a working variable
                y = y + m * 400;
                d = floorModUseFloorDivResult(d, 146097, m);
            }
            m = floorDivBeware(d * 4 - 7468865, 146097);
            d = d + 1525 + m - floorDivBeware(m, 4);  // +1525 == +1524 +1
        }
        m = (d * 4 - 489) / 1461; // here "m" is a working variable
        y = y + m;
        d = d - (m * 1461) / 4;
        m = (d * 5 - 1) / 153;
        d = d - (m * 153) / 5;
        if (m < 14) {
            m = m - 1;
        } else {
            m = m - 13;
            y = y + 1;
        }
        int[] ymdv = new int[3];
        ymdv[0] = y;
        ymdv[1] = m;
        ymdv[2] = d;
        return ymdv;
    }

    public static int[] cjdnToYMDuseMeeusFloat(int cjdnv, double firstGregorianCJDNv) {
        // *** uses floating point Meeus ***
        double yy, mm, dd;  // use fewer variables than usual for this algorithm
        if (cjdnv = 14) {
            mm = mm - 13;
            yy = yy - 4715;
        } else {
            mm = mm - 1;
            yy = yy - 4716;
        }
        int[] ymdv = new int[3];
        ymdv[0] = (int) yy;
        ymdv[1] = (int) mm;
        ymdv[2] = (int) dd;
        return ymdv;
    }

    public static int ymdToCJDNinteger(int y, int m, int d, double firstGregorianCJDNv) {
        // note this does not check for integer overflow;
        int b, cjdnv;
        if (m < 3) {
            y = y - 1;
            m = m + 12;
        }
        b = floorDivBeware(y, 100);
        b = 2 - b + floorDivBeware(b, 4);
        y = y + 4716;
        cjdnv = 365 * y + floorDivBeware(y, 4) + (306 * (m + 1)) / 10 + d + b - 1524;
        if (cjdnv < firstGregorianCJDNv) {
            cjdnv = cjdnv - b;
        }
        return cjdnv;
    }

    public static int ymdToCJDNfloat(int y, int m, int d, double firstGregorianCJDNv) {
        // note this does not check for integer overflow in the return value;
        double bb, cjdnvv;
        if (m < 3) {
            y = y - 1;
            m = m + 12;
        }
        bb = Math.floor(y / 100.0);
        bb = 2 - bb + Math.floor(bb/ 4.0);
        cjdnvv = Math.floor(365.25 * (y + 4716)) + Math.floor(30.6001 * (m + 1)) +
                 d + bb - 1524;
        if (cjdnvv < firstGregorianCJDNv) {
            cjdnvv = cjdnvv - bb;
        }
        return (int) cjdnvv;
    }

    public static int[] cjdnToYMDintegerEtc(int cjdnv, double firstGregorianCJDNv) {
        int y, m, d;
        y = -4712;
        // The constants 500100100 and 5000 have margins to allow small changes
        // to this algorithm without causing integer overflow problems.
        if (cjdnv  500100100 || d  500100100 || d < 5000) {
                m = floorDivBeware(d, 146097); // number of 400 year Gregorian cycles
                y = y + m * 400 - 400;
                d = floorModUseFloorDivResult(d, 146097, m) + 146097;
            }
            m = floorDivBeware(d * 4 - 7468865, 146097);
            d = d - 59 + m - floorDivBy4(m);  // -59 == -60 +1
        }
        m = (d * 4 + 3) / 1461; // "m" is a working variable
        y = y + m;
        d = d - floorDivBy4(m * 1461);                             // now 0 <= d < 366
        m = dayInYearFirstDay0Base0301ToMonthBeware(d);            // now 1 <= m <= 12
        d = d + 1 - DaysInYearBeforeMonthBaseMonth3ArrayBeware[m]; // now 1 <= d <= 31
        if (m < 3) {
            y = y + 1;
        }
        int[] ymdv = new int[3];
        ymdv[0] = y;
        ymdv[1] = m;
        ymdv[2] = d;
        return ymdv;
    }

    public static int[] cjdnToYMDintegerShortcutsEtc(int cjdnv, double firstGregorianCJDNv) {
        int y, m, d;
        y = -4712;
        // The constants 500100100 and 5000 have margins to allow small changes
        // to this algorithm without causing integer overflow problems.
        if (cjdnv  500100100 || d < 5000) {
                m = floorDivBeware(d, 1461); // here "m" is a working variable
                y = y + m * 4 - 4;
                d = floorModUseFloorDivResult(d, 1461, m) + 1461;
            }
            d = d - 60; // rebase from -4712-01-01 CJDN 0 to -4712-03-01 CJDN 60;
        } else if (2415080 <= cjdnv && cjdnv < 2488129) {
            // shortcut if Gregorian 1900-03-01 <= cjdnv  500100100 || d < 5000) {
                m = floorDivBeware(d, 146097); // number of 400 year Gregorian cycles
                y = y + m * 400 - 400;
                d = floorModUseFloorDivResult(d, 146097, m) + 146097;
            }
            m = floorDivBeware(d * 4 - 7468865, 146097);
            d = d - 59 + m - floorDivBy4(m);  // -59 == -60 +1
        }
        m = (d * 4 + 3) / 1461; // "m" is a working variable
        y = y + m;
        d = d - floorDivBy4(m * 1461);                             // now 0 <= d < 366
        m = dayInYearFirstDay0Base0301ToMonthBeware(d);            // now 1 <= m <= 12
        d = d + 1 - DaysInYearBeforeMonthBaseMonth3ArrayBeware[m]; // now 1 <= d <= 31
        if (m < 3) {
            y = y + 1;
        }
        int[] ymdv = new int[3];
        ymdv[0] = y;
        ymdv[1] = m;
        ymdv[2] = d;
        return ymdv;
    }

    public static int[] cjdnToYMDuseMeeusFloatShortcuts(int cjdnv, double firstGregorianCJDNv) {
        // *** uses floating point Meeus ***
        double yy, mm, dd;  // use fewer variables than usual for this algorithm
        if (cjdnv < firstGregorianCJDNv) {
            dd = cjdnv;
        } else if (2415080 <= cjdnv && cjdnv < 2488129) {
            // shortcut if Gregorian 1900-03-01 <= cjdnv = 14) {
            mm = mm - 13;
            yy = yy - 4715;
        } else {
            mm = mm - 1;
            yy = yy - 4716;
        }
        int[] ymdv = new int[3];
        ymdv[0] = (int) yy;
        ymdv[1] = (int) mm;
        ymdv[2] = (int) dd;
        return ymdv;
    }


    public static int ymdToCJDNintegerEtc(int y, int m, int d,
                                          double firstGregorianCJDNv) {
        // note this does not check for integer overflow;
        int b, cjdnv;
        if (m < 3) {
            y = y - 1;
            m = m + 12;
        }
        b = floorDivBeware(y, 100);
        b = 2 - b + floorDivBy4(b);
        y = y + 4712;
        // in next line "+59" is "+60 -1", where 60 is CJDN of Julian -4712-03-01;
        cjdnv = 365 * y + floorDivBy4(y) + d + b + 59 +
                DaysInYearBeforeMonthBaseMonth3ArrayBeware[m];
        if (cjdnv < firstGregorianCJDNv) {
            cjdnv = cjdnv - b;
        }
        return cjdnv;
    }

    public static int ymdToCJDNintegerShortcutsEtc(int y, int m, int d,
                                                   double firstGregorianCJDNv) {
        // note this does not check for integer overflow;
        int b, cjdnv;
        if (m < 3) {
            y = y - 1;
            m = m + 12;
        }
        if (1900 < y && y < 2099) {
            b = -13;
        } else {
            b = floorDivBeware(y, 100);
            b = 2 - b + floorDivBy4(b);
        }
        y = y + 4712;
        // in next line "+59" is "+60 -1", where 60 is CJDN of Julian -4712-03-01;
        cjdnv = 365 * y + floorDivBy4(y) + d + b + 59 +
                DaysInYearBeforeMonthBaseMonth3ArrayBeware[m];
        if (cjdnv < firstGregorianCJDNv) {
            cjdnv = cjdnv - b;
        }
        return cjdnv;
    }

    public static int ymdToCJDNfloatShortcuts(int y, int m, int d,
                                              double firstGregorianCJDNv) {
        // note this does not check for integer overflow in the return value;
        double bb, cjdnvv;
        if (m < 3) {
            y = y - 1;
            m = m + 12;
        }
        if (1900 < y && y < 2099) {
            bb = -13;
        } else {
            bb = Math.floor(y / 100.0);
            bb = 2 - bb + Math.floor(bb/ 4.0);
        }
        cjdnvv = Math.floor(365.25 * (y + 4716)) + Math.floor(30.6001 * (m + 1)) +
                 d + bb - 1524;
        if (cjdnvv < firstGregorianCJDNv) {
            cjdnvv = cjdnvv - bb;
        }
        return (int) cjdnvv;
    }

    // days in a year before the month with base March(month 3).01;
    // allows m = 0, 1 to 12, 13, 14; the 0 & 13 & 14 make ymd_to_cjdn a bit faster;
    // "Beware" is because for month 0, 1, 2 it returns the values for 12, 13, 14;
    // initial capital letter is used in name to make conversion to Ruby code easier;
    public static final int[] DaysInYearBeforeMonthBaseMonth3ArrayBeware =
                                  { 275, 306, 337,
                                           0,  31,  61,  92, 122,
                                         153, 184, 214, 245, 275,
                                         306, 337 };

    public static int dayInYearFirstDay0Base0301ToMonthBeware(int d) {
        // This is for the case 0 <= d = 306)
        // this returns 1 or 2 instead of 13 or 14.
        if (d < 122) {
            if (d < 61) {
                if (d < 31) {
                    return 3;
                }
                return 4;
            } else if (d < 92) {
                return 5;
            }
            return 6;
        }
        if (d < 245) {
            if (d < 184) {
                if (d < 153) {
                    return 7;
                }
                return 8;
            } else if (d < 214) {
                return 9;
            }
            return 10;
        }
       if (d < 306) {
            if (d < 275) {
                return 11;
            }
            return 12;
        }
        if (d < 337) {
            return 1;
        }
        return 2;
    }

    public static double Infinity = 1.0 / 0.0;
    public static double MinusInfinity = -1.0 / 0.0;

    public static boolean isLeapYear(int yearv, double firstGregorianCJDNv) {
      if (isNotDivisibleBy4(yearv)) return false;
      if (1901 <= yearv && yearv = 0) {
        iy = yearv / 100;
      } else if (yearv > -100) {
        iy = 0;
      } else {
        iy = -((-(yearv + 100)) / 100) - 1; // avoid integer overflow
      }
      if (iy * 100 != yearv || isDivisibleBy4(iy)) return true;
      // if here yearv is divisible by 100 but not by 400
      if (firstGregorianCJDNv == Infinity) return true; // all Julian calendar
      if (firstGregorianCJDNv == MinusInfinity) return false; // all Gregorian calendar
      // If here yearv is a leap year if firstGregorianCJDNv >= CJDN of yearv.03.01.
      // In next line cast to 64bit integer (long) to avoid integer overflows.
      // We can combine (1721120 + ...) as 32bit integer which will not overflow.
      return firstGregorianCJDNv >= (1721120 + floorDivBy4(iy)) + ((long) iy) * 36524;
    }

    public static int INT32_MIN_VALUE = -2147483648;  // -(2 to power 31)

    public static int floorDivBeware(int a, int b) {
        // This is useful in general calendar calculations.
        // It returns result of dividing a by b: if b does not exactly divide a
        // the result is rounded to the nearest integer  0 then we need not check the result
        // for overflow or divide by zero.
        //
        // This is for 32bit signed integers using the standard method of
        // representing negative integers as "complements" of positive integers.
        // It can be easily adapted for other sizes of signed integers which
        // use the standard method of representing negative integers
        // as "complements" of positive integers.
        //
        // This is written in Java, but it can also be used for C because
        // it assumes that the only well defined results of (a / b) are
        // for (a >= 0) and (b > 0). Where (a / b) is well defined for other
        // values of a and b (for example in Java) the code could be shortened.
        //
        if (b > 0) {           // most likely case, so check for b > 0 first;
            if (a >= 0) {
                return a / b;
            }
            return -((-(a + 1)) / b) - 1;
        }
        // Now deal with b  0) {
                return -1;
            }
            if (a != INT32_MIN_VALUE) {
                return 0;  // INT32_MIN_VALUE == b < a <= 0;
            }
            return 1;  // INT32_MIN_VALUE == b == a;
        }
        // Now deal with INT32_MIN_VALUE < b <= 0.
        if (b  0) {
                return -((a - 1) / (-b)) - 1;
            }
            if (a != INT32_MIN_VALUE) {
                return (-a) / (-b);
            }
            return (-(a - b)) / (-b) + 1;  // INT32_MIN_VALUE == a < b  0) {           // most likely case, so check for b > 0 first;
            // if here then b > 0 and we should have floorMod(a, b) >= 0;
            if (a >= 0) {
                return a - b * floorDivResult;  // a >= 0 and floorDivResult >= 0
            }
        } else if (a < 0) {
            // if here then b < 0 and we should have floorMod(a, b) <= 0;
            return a - b * floorDivResult;  // a < 0 && b < 0
        }
        // if here: a  0 and floorDivResult  0 and b < 0 and floorDivResult = 0) {
            return a >> 2;
        }
        return -((-(a + 1)) >> 2) - 1;
    }

    public static int floorModBy4(int a) {
        return a & 3;
    }

    public static int floorModBy4UseFloorDivResult(int a, int floorDivResult) {
        // Note that we do not need floorDivResult: this method is here
        // so if someone guesses that the divide by 4 equivalent of
        // floorModUseFloorDivResult is floorModBy4UseFloorDivResult
        // then their guess will not be wrong.
        // And for some languages we may not be able to use (a & 3), and must do:
        //     if (a >= 0) {
        //         return a - 4 * floorDivResult;
        //     }
        //     return (a - 4 * (floorDivResult + 1)) + 4;
        return a & 3;
    }

}
Advertisements

About Colin Bartlett

I'm interested in arts, mathematics, science. Suliram is a partial conflation of the names of three good actors: Ira Aldridge, Anna May Wong, and another. My intention is to use a personal experience of arts to make some points, but without being too "me me me" about it. And to follow Strunk's Elements of Style. Except that I won't always "be definite": I prefer Niels Bohr's precept that you shouldn't write clearer than you think.
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s