[ruby/date] Optimize Gregorian date conversions with Neri-Schneider algorithm

Replace floating-point arithmetic and iterative loops with pure integer
operations for ~40% faster Date operations. Date.ordinal and Date.commercial
are ~2x faster due to O(1) first-day-of-year calculation.

Reference: https://arxiv.org/abs/2102.06959

Author: mensfeld

ext/date/date_core.c

@@ -7,6 +7,7 @@
 #include "ruby/util.h"
 #include <math.h>
 #include <time.h>
+#include <stdint.h>  /* For uint64_t in Neri-Schneider algorithm */
 #if defined(HAVE_SYS_TIME_H)
 #include <sys/time.h>
 #endif
@@ -452,11 +453,27 @@ do {\
 static int c_valid_civil_p(int, int, int, double,
 			   int *, int *, int *, int *);
 
+/* Forward declarations for Neri-Schneider optimized functions */
+static int c_gregorian_civil_to_jd(int y, int m, int d);
+static void c_gregorian_jd_to_civil(int jd, int *ry, int *rm, int *rd);
+inline static int c_gregorian_fdoy(int y);
+inline static int c_gregorian_ldoy(int y);
+inline static int c_gregorian_ldom_jd(int y, int m);
+inline static int ns_jd_in_range(int jd);
+
 static int
 c_find_fdoy(int y, double sg, int *rjd, int *ns)
 {
     int d, rm, rd;
 
+    /* Fast path: pure Gregorian calendar */
+    if (isinf(sg) && sg < 0) {
+	*rjd = c_gregorian_fdoy(y);
+	*ns = 1;
+	return 1;
+    }
+
+    /* Keep existing loop for Julian/reform period */
     for (d = 1; d < 31; d++)
 	if (c_valid_civil_p(y, 1, d, sg, &rm, &rd, rjd, ns))
 	    return 1;
@@ -468,6 +485,14 @@ c_find_ldoy(int y, double sg, int *rjd, int *ns)
 {
     int i, rm, rd;
 
+    /* Fast path: pure Gregorian calendar */
+    if (isinf(sg) && sg < 0) {
+	*rjd = c_gregorian_ldoy(y);
+	*ns = 1;
+	return 1;
+    }
+
+    /* Keep existing loop for Julian/reform period */
     for (i = 0; i < 30; i++)
 	if (c_valid_civil_p(y, 12, 31 - i, sg, &rm, &rd, rjd, ns))
 	    return 1;
@@ -493,6 +518,14 @@ c_find_ldom(int y, int m, double sg, int *rjd, int *ns)
 {
     int i, rm, rd;
 
+    /* Fast path: pure Gregorian calendar */
+    if (isinf(sg) && sg < 0) {
+	*rjd = c_gregorian_ldom_jd(y, m);
+	*ns = 1;
+	return 1;
+    }
+
+    /* Keep existing loop for Julian/reform period */
     for (i = 0; i < 30; i++)
 	if (c_valid_civil_p(y, m, 31 - i, sg, &rm, &rd, rjd, ns))
 	    return 1;
@@ -502,55 +535,74 @@ c_find_ldom(int y, int m, double sg, int *rjd, int *ns)
 static void
 c_civil_to_jd(int y, int m, int d, double sg, int *rjd, int *ns)
 {
-    double a, b, jd;
+    int jd;
 
-    if (m <= 2) {
-	y -= 1;
-	m += 12;
+    /* Fast path: pure Gregorian calendar (sg == -infinity) */
+    if (isinf(sg) && sg < 0) {
+	*rjd = c_gregorian_civil_to_jd(y, m, d);
+	*ns = 1;
+	return;
     }
-    a = floor(y / 100.0);
-    b = 2 - a + floor(a / 4.0);
-    jd = floor(365.25 * (y + 4716)) +
-	floor(30.6001 * (m + 1)) +
-	d + b - 1524;
+
+    /* Calculate Gregorian JD using optimized algorithm */
+    jd = c_gregorian_civil_to_jd(y, m, d);
+
     if (jd < sg) {
-	jd -= b;
+	/* Before Gregorian switchover - use Julian calendar */
+	int y2 = y, m2 = m;
+	if (m2 <= 2) {
+	    y2 -= 1;
+	    m2 += 12;
+	}
+	jd = (int)(floor(365.25 * (y2 + 4716)) +
+		   floor(30.6001 * (m2 + 1)) +
+		   d - 1524);
 	*ns = 0;
     }
-    else
+    else {
 	*ns = 1;
+    }
 
-    *rjd = (int)jd;
+    *rjd = jd;
 }
 
 static void
 c_jd_to_civil(int jd, double sg, int *ry, int *rm, int *rdom)
 {
-    double x, a, b, c, d, e, y, m, dom;
-
-    if (jd < sg)
-	a = jd;
-    else {
-	x = floor((jd - 1867216.25) / 36524.25);
-	a = jd + 1 + x - floor(x / 4.0);
-    }
-    b = a + 1524;
-    c = floor((b - 122.1) / 365.25);
-    d = floor(365.25 * c);
-    e = floor((b - d) / 30.6001);
-    dom = b - d - floor(30.6001 * e);
-    if (e <= 13) {
-	m = e - 1;
-	y = c - 4716;
-    }
-    else {
-	m = e - 13;
-	y = c - 4715;
+    /* Fast path: pure Gregorian or date after switchover, within safe range */
+    if (((isinf(sg) && sg < 0) || jd >= sg) && ns_jd_in_range(jd)) {
+	c_gregorian_jd_to_civil(jd, ry, rm, rdom);
+	return;
     }
 
-    *ry = (int)y;
-    *rm = (int)m;
-    *rdom = (int)dom;
+    /* Original algorithm for Julian calendar or extreme dates */
+    {
+	double x, a, b, c, d, e, y, m, dom;
+
+	if (jd < sg)
+	    a = jd;
+	else {
+	    x = floor((jd - 1867216.25) / 36524.25);
+	    a = jd + 1 + x - floor(x / 4.0);
+	}
+	b = a + 1524;
+	c = floor((b - 122.1) / 365.25);
+	d = floor(365.25 * c);
+	e = floor((b - d) / 30.6001);
+	dom = b - d - floor(30.6001 * e);
+	if (e <= 13) {
+	    m = e - 1;
+	    y = c - 4716;
+	}
+	else {
+	    m = e - 13;
+	    y = c - 4715;
+	}
+
+	*ry = (int)y;
+	*rm = (int)m;
+	*rdom = (int)dom;
+    }
 }
 
 static void
@@ -725,6 +777,113 @@ c_gregorian_last_day_of_month(int y, int m)
     return monthtab[c_gregorian_leap_p(y) ? 1 : 0][m];
 }
 
+/*
+ * Neri-Schneider algorithm for optimized Gregorian date conversion.
+ * Reference: Neri & Schneider, "Euclidean Affine Functions and Applications
+ *            to Calendar Algorithms", Software: Practice and Experience, 2023.
+ *            https://arxiv.org/abs/2102.06959
+ *
+ * This algorithm provides ~2-3x speedup over traditional floating-point
+ * implementations by using pure integer arithmetic with multiplication
+ * and bit-shifts instead of expensive division operations.
+ */
+
+/* JDN of March 1, Year 0 in proleptic Gregorian calendar */
+#define NS_GREGORIAN_EPOCH 1721120
+
+/*
+ * Safe bounds for Neri-Schneider algorithm to avoid integer overflow.
+ * These correspond to approximately years -1,000,000 to +1,000,000.
+ */
+#define NS_JD_MIN -364000000
+#define NS_JD_MAX  538000000
+
+inline static int
+ns_jd_in_range(int jd)
+{
+    return jd >= NS_JD_MIN && jd <= NS_JD_MAX;
+}
+
+/* Optimized: Gregorian date -> Julian Day Number */
+static int
+c_gregorian_civil_to_jd(int y, int m, int d)
+{
+    /* Shift epoch to March 1 of year 0 (Jan/Feb belong to previous year) */
+    int j  = (m < 3) ? 1 : 0;
+    int y0 = y - j;
+    int m0 = j ? m + 12 : m;
+    int d0 = d - 1;
+
+    /* Calculate year contribution with leap year correction */
+    int q1 = DIV(y0, 100);
+    int yc = DIV(1461 * y0, 4) - q1 + DIV(q1, 4);
+
+    /* Calculate month contribution using integer arithmetic */
+    int mc = (979 * m0 - 2919) / 32;
+
+    /* Combine and add epoch offset to get JDN */
+    return yc + mc + d0 + NS_GREGORIAN_EPOCH;
+}
+
+/* Optimized: Julian Day Number -> Gregorian date */
+static void
+c_gregorian_jd_to_civil(int jd, int *ry, int *rm, int *rd)
+{
+    int r0, n1, q1, r1, n2, q2, r2, n3, q3, r3, y0, j;
+    uint64_t u2;
+
+    /* Convert JDN to rata die (March 1, Year 0 epoch) */
+    r0 = jd - NS_GREGORIAN_EPOCH;
+
+    /* Extract century and day within 400-year cycle */
+    /* Use Euclidean (floor) division for negative values */
+    n1 = 4 * r0 + 3;
+    q1 = DIV(n1, 146097);                   /* Century */
+    r1 = MOD(n1, 146097) / 4;               /* Day within 400-year cycle */
+
+    /* Use 64-bit arithmetic for year calculation within century */
+    n2 = 4 * r1 + 3;
+    u2 = (uint64_t)2939745 * (uint64_t)n2;
+    q2 = (int)(u2 >> 32);                   /* Year within century */
+    r2 = (int)((uint32_t)u2 / 2939745 / 4); /* Day of year */
+
+    /* Calculate month and day using integer arithmetic */
+    n3 = 2141 * r2 + 197913;
+    q3 = n3 >> 16;                          /* Month (3-14) */
+    r3 = (n3 & 0xFFFF) / 2141;              /* Day of month (0-based) */
+
+    /* Combine century and year */
+    y0 = 100 * q1 + q2;
+
+    /* Adjust for January/February (shift from fiscal year) */
+    j = (r2 >= 306) ? 1 : 0;
+
+    *ry = y0 + j;
+    *rm = j ? q3 - 12 : q3;
+    *rd = r3 + 1;
+}
+
+/* O(1) first day of year for Gregorian calendar */
+inline static int
+c_gregorian_fdoy(int y)
+{
+    return c_gregorian_civil_to_jd(y, 1, 1);
+}
+
+/* O(1) last day of year for Gregorian calendar */
+inline static int
+c_gregorian_ldoy(int y)
+{
+    return c_gregorian_civil_to_jd(y, 12, 31);
+}
+
+/* O(1) last day of month (JDN) for Gregorian calendar */
+inline static int
+c_gregorian_ldom_jd(int y, int m)
+{
+    return c_gregorian_civil_to_jd(y, m, c_gregorian_last_day_of_month(y, m));
+}
+
 static int
 c_valid_julian_p(int y, int m, int d, int *rm, int *rd)
 {