Faster `Date::from_julian_day` using Ben Joffe's algorithms?

[madsmtm]

I stumbled upon this article on Hack News today.

In it, Ben Joffe discusses a lot of tricks to speed up int -> date parsing, specifically the calculation in Date::from_julian_day_unchecked that time currently implements using Peter Baum's algorithms IIUC. Ben Joffe also has two earlier articles detailing tricks for operating in the 32-bit range (which are probably more relevant in time's case).

Idk., I thought you might be interested in it, especially since operations like Date::checked_add are implemented in terms of this conversion - but feel free to close this issue if you don't think it's worth it.

[jhpratt]

Yeah, the Julian day (technically Julian day number) is used all over the place internally, so this is something worth looking in to.

Off the top of my head, I believe that Peter Baum's algorithm is what is used, though I modified it slightly to work with year-ordinal instead of year-month-day. This is now reflected in his paper, as I reached out to him at the time.

When I get a chance, I'll read through all three articles. Just by reading the introduction from the third, I can't immediately dismiss it for the reason I've dismissed every other algorithm. Many only work with positive years (or even after 1970), but I need negative as well. Given the claim of ±1.89 trillion years, that's certainly more than sufficient. In any situation, I likely have a script lying around to validate the full range such that I don't have to manually check things.

[benjoffe]

Hi, I'm the author of those articles and have stumbled upon this thread.

As you know, the algorithms in those blog posts are optimised for the YMD case, whereas this library appears to primarily use Year/Ordinal/Leap. A direct translation of those techniques will likely speed things up a little, but won't be optimal

The next blog post in my series will be focused mainly on 32-bit ordinal date algorithms, and will introduce an entirely new approach which should be much faster for this case. I will be happy to bump this thread once the algorithms are finalised and the blog post is ready, which I expect to be in the next week or two.

Once that is done, if you like, and time permitting, I will also be happy to create a pull request with the algorithm implementation.

I have reviewed the existing implementation of from_julian_day_unchecked and it looks like Peter Baum's algorithm was swapped out for Neri-Schneider around 10 months ago in
commit ec327a2 - although the comment and documentation were not updated.

On a related note, would you be aware of any other date libraries that use Year/Ordinal as the primary date implementation? It is a very interesting alternative to the more common YMD implementation, and comes with a unique set of advantages. Reviewing these will help make my comparison benchmarks more broad, and might surface ideas I haven’t considered yet.

Thank you for your time.

[jhpratt]

Looks like I did swap out the implementation; good catch. I'll keep an eye out for your next post, though I'm certainly not going to stop you from posting it here if you'd like.

On a related note, would you be aware of any other date libraries that use Year/Ordinal as the primary date implementation? It is a very interesting alternative to the more common YMD implementation, and comes with a unique set of advantages. Reviewing these will help make my comparison benchmarks more broad, and might surface ideas I haven’t considered yet.

It's been quite a while since I looked into other datetime libraries across various languages, so the notes I took are long gone. I do recall running benchmarks when I chose this layout that showed an improvement for most functions (circa 10% iirc) with only minimal harm to year-month-day calculations, so I deemed it to be worth it. Without anything concrete in front of me, it makes quite a few calculations simpler as you ultimately only care about the number of days and don't have to account for the varying lengths of months.

[benjoffe]

Hi,

I have committed a draft of the new algorithm here:
https://github.com/benjoffe/fast-date-benchmarks/blob/main/algorithms_ordinal/ordinal_benjoffe_fast32.hpp

There is not much commentary yet, so I will note a couple of things:

• It is accurate over more than 40% of the 32-bit range, which exceeds the range of the existing implementation, which is around 25%
• My preliminary benchmarks indicate about 20% speedup on some computers, and 50% speedup on Apple Silicon (i.e. up to 2x speed). This is in C++, things may differ in Rust.

This may not be the final version, but I think it's pretty good. I was side-tracked creating a new blog post about a faster full-range 32-bit isLeapYear function, which is now online. This may be of interest as well:
https://www.benjoffe.com/fast-leap-year

Since I haven't yet found any other date libraries that use year/ordinal/leap as the primary date representation, the blog post about this algorithm, when it's finished, will be primarily referencing this repo. As such, I think it will make for a much more positive tone if it is framed as "How Rust Time is now X% faster in some places", rather than just showing that Rust Time is slow compared to this one. That is to say, I am happy to delay the blog post until this is implemented in this library, if you like.

From a cursory look earlier, I think there may be some easy speed gains by ensuring the other code to calculate day-of-month and month-of-year reuse the 'leap' variable, rather than recalculate it (although I haven't looked deep to see if that actually is possible).

While I am happy to help with a pull request, it might be several weeks before I can have time, so feel free to try implementing it yourself if you are interested. I'll be happy to answer any questions about it.

[jhpratt]

Checking the leap year algorithm, it is a statistically significant improvement. However, I need to refine my benchmark before I assign a concrete number and switch the implementation.

For the ordinal date, I'll give it a shot soon!

[jhpratt]

I'll be switching the leap year implementation. After thorough benchmarking, it is 0.6% faster when i32 is used as a range. That's compared to % 25, where I was still using the older % 100 (that gave an additional 3.4% improvement).

A quick attempt at porting the ordinal date from Julian day number to Rust resulted in a number of failing tests, and nothing immediately stands out to me as to why.

[benjoffe]

Hi,
Thanks for trying out the algorithm.
If you want to share a copy of the initial implementation I'll be happy to debug it.

[jhpratt]

Swapping out the current implementation in date.rs for

pub(crate) const unsafe fn from_julian_day_unchecked(julian_day: i32) -> Self {
    debug_assert!(julian_day >= Self::MIN.to_julian_day());
    debug_assert!(julian_day <= Self::MAX.to_julian_day());

    const ERAS: u32 = 5949;
    // Rata Die shift:
    const D_SHIFT: u32 = 146097 * ERAS + 719162 + 366;
    // Year shift:
    const Y_SHIFT: u32 = 400 * ERAS;

    const CEN_MUL: u32 = ((4u64 << 47) / 146097) as u32;
    const JUL_MUL: u32 = ((4u64 << 40) / 1461 + 1) as u32;
    const CEN_CUT: u32 = ((365u64 << 32) / 36525) as u32;

    let day = julian_day.wrapping_add_unsigned(D_SHIFT) as u32; // Epoch: -XX00-01-01
    let c_n = (day as u64 * CEN_MUL as u64) >> 15; // Divide 36524.25
    let cen = (c_n >> 32) as u32; // Century
    let cpt = c_n as u32; // Century-part
    let ijy = (cen % 4 == 0) || (cpt > CEN_CUT); // "Is Julian Year"
    let jul = day.wrapping_sub(cen / 4).wrapping_add(cen); // Julian Map
    let y_n = (jul as u64 * JUL_MUL as u64) >> 8; // Divide 365.25
    let yrs = (y_n >> 32) as u32; // Year
    let ypt = y_n as u32; // Year-part

    let year = yrs.wrapping_sub(Y_SHIFT) as i32;
    let ordinal = ((ypt as u64 * 1461) >> 34) as u32 + ijy as u32;
    let leap = (yrs % 4 == 0) & ijy;

    unsafe { Self::from_parts(year, leap, ordinal as u16) }
}

and running the tests

cargo test -p time --all-features --tests -q

gives 101 failures. That number isn't as bad as it seems — tons of things rely on this function.

As the tests are run in debug mode, the assertions here and in from_parts (the final, unsafe call) are executed; none are triggered so that's not the issue.

[benjoffe]

Thanks for sharing.

It looks like the epoch is the issue.

I subtracted 2_440_588 from D_SHIFT and it then passed all tests:

const D_SHIFT: u32 = 146097 * ERAS + 719162 + 366 - 2_440_588;

[jhpratt]

🤦 I do an identical adjustment in the current implementation. It's needed to account for the difference in epoch, which for some reason I didn't think twice about the comment in your code.

[benjoffe]

Yep, I grabbed that number from the previous code.

I'll explain the algorithm a little since it's probably not nice to include code that is mystical.

First, I'll define a Gregorian Year as a year divisible by 100 but not divisible by 400. Eg. 100,200,300,500... etc.
That is, years that are not "Julian Years" - following logic of the Julian Calendar (all other years).

• day is a remapped rata die counting from an epoch some large multiple of 400 years in the past. The value of ERAS can be tweaked to adjust which particular year it counts up from, so that the range of supported years cover the range required by the function. I chose 5949 as it seemed to give a balanced range centred around the year zero, but it can be tweaked to give other ranges (such as those balanced around the year 1970 etc.). I assume this library won't need to tweak ERAS since the covered range is already in millions of years, which I believe exceeds the requirements.
• c_n is calculated from a fixed-point multiplication that places the number of elapsed centuries in the upper 32-bits of the 64-bit result. It bit-shifts by 15 because CEN_MUL uses a bit-shift of 47, which is 15 more than 32. CEN_MUL uses a bit-shift of 47 because that gives the largest possible fixed-point multiplier that fits in a 32-bit number.
• cpt is the lower 32 bits, and represents the amount of the century that has elapsed, where zero = start of century, and 2^32-1 is at the end of the century.
• ijy is a boolean that tells us whether the calendar year matches the rules of the Julian Calendar (i.e not the Gregorian rules). This is true when we are past the first year of the century (CEN_CUT is the cutoff point for 1-year out of 100 years in 32-bit space, as used by cpt) - but is ALSO true when the century is divisible by 4. The only times this is NOT true is in years divisible by 100 that are not divisible by 400. This is the trick that allows this new algorithm to avoid calling isLeap(year).
• jul is the current rata die remapped to Julian year logic (as outlined in my first article) - with one exception: it is inserting a fake "Jan 0" instead of a fake "Feb 29". This will be corrected for the in calculation of ordinal.
• y_n is calculated from a fixed-point multiplication the same way c_n is performed, except done for years instead of centuries. This bit-shift by 8 is for the same reason as before, but this time because JUL_MUL has a bit-shift of 40, which exceeds 32 by 8.
• yrs is the upper 32-bits (number of years elapsed)
• ypt is the lower 32-bits, where 0 is start of year and 2^32-1 is end of year
• ordinal is mapping ypt by doing 365.25 / 2^32, but multiplying numerator and denominator by 4 so we can avoid fractional multiplication. We usually want to add 1 since this is 1-indexed and not 0-indexed, but we also want to subtract by 1 in the Gregorian special cases of 100,200,300,500... That is why we add ijy - it achieves those two requirements. The reason Gregorian special cases should have subtraction by 1 is because as noted previously in the calculation of jul, the timeline is being padded at the start of the year in Gregorian years already, which means we don't need to pad by 1 in those years.
• leap is the last part of the trick to avoid directly computing whether it's a leap year. i.e. years divisible by 4 that are not those Gregorian special cases of 100,200,300,500...

I note that the fake padding at the start of the year is the reason that this algorithm doesn't need to use 1 March as the starting point, and not benefit from calculating dates backwards (as was done in my 64-bit algorithm). This is uniquely possible due to using ordinal dates, rather than y/m/d.

I note that this algorithm is MOSTLY 32-bit friendly, but will involve 4 extra cycles on 32-bit computers, which is pretty minor.
The reason for the slight penalty on 32-bit is the bit-shifts in the calculation of c_n and y_n involve shifting a 64-bit result (shift by 15 and 8 respectively). Compiler will typically not perform that bit-shift there, and merge the bit-shifts with the subsequent attempts to get the upper and lower 32-bit results. The upper 32-bit result with thus be a single-cycle bit-shift of 15 and 8 respectively of the upper half result (just as cheap), but the lower-32 bit results will involve the cpu taking the upper result and shifting it LEFT, and taking the lower result and shifting it RIGHT, and adding the results. The above is thus 2 extra cycles for both of those steps (4 in total). I suspect that even with this slight penalty on 32-bit computers the result is still faster overall, as the speed gain should be much more than just 4 cycles.

I will be interested to see if this algorithm benchmarks much faster for you.

[benjoffe]

Hi,

If you are interested, I have worked on a faster version of to_calendar_date, which can also be used for the standalone month() and day(). There seems to be a notable speedup for each, although day() may not change much. I presume that most performance critical code will use the to_calendar_date to print the entire date in one go?

This is the performance improvement to to_calendar_date measured on M4 Mac Pro:

Date: to_calendar_date  time:   [22.044 ns 22.066 ns 22.088 ns]
                        change: [-15.501% -15.362% -15.221%] (p = 0.00 < 0.00)
                        Performance has improved.

It builds upon your existing algorithm, which itself is an impressive algorithm as it appears to have independently discovered some neat tricks. The changes are:

1. it uses a larger scale factor so the modulus can be 2^16 on x86, which is "free".
2. it pre-computes some constants so that there are fewer additions in relation to leap years and so on. Additionally, it is structured so that the first multiplication can initiate immediately on super-scalar processors, and allow the leap year checks and so-on to be done in parallel (processor permitting)
3. it uses the high and low bits from the multiplication, which is the technique pioneered by Neri-Schneider

The algorithm in Rust is as follows (Note: I haven't implemented the smaller SCALE on ARM as I'm not sure what the correct way to check for ARM would be in this codebase, I also don't know if you are interested in such platform specific optimisations):

        let (year, ordinal) = self.to_ordinal_date();
        let ordinal = ordinal as u32;
        let is_leap = self.is_in_leap_year();
        let jan_feb_len = 59 + is_leap as u32;

        const SCALE: u32 = 2; // TODO: make this `1` for ARM devices
        const STEP: u32 = 1071 * SCALE;
        const DIVISOR: u32 = SCALE << 15;
        const SHIFT_0: u32 = DIVISOR - 439 * SCALE;
        const SHIFT_1: u32 = SHIFT_0 + STEP;
        const SHIFT_2: u32 = SHIFT_1 + STEP;

        let shift = if ordinal <= jan_feb_len {
            SHIFT_0
        }
        else if is_leap {
            SHIFT_1
        } else {
            SHIFT_2
        };
        
        let num = ordinal * STEP + shift;
        let month = num / DIVISOR;
        let day = (num % DIVISOR) / STEP + 1;

A c++ implementation is also available to review here:
https://github.com/benjoffe/fast-date-benchmarks/blob/main/algorithms/benjoffe_ordinal_alternative.hpp

I have also put a draft blog page up that documents the two algorithms I have shared in this thread:
https://www.benjoffe.com/fast-ordinal-date

If you feel inclined to share a URL in either the documentation or codebase for any of these algorithms, I would ask that you please use that URL, as that will be the canonical representation of these algorithms and the article will be iterated on and properly released if/once this library uses them.