diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py
index 9a74b0cb..7c4650b3 100644
--- a/demo/python/astronomy.py
+++ b/demo/python/astronomy.py
@@ -577,7 +577,7 @@ def _UniversalTime(tt):
             return ut
         dt += err
 
-_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
+_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
 
 class Time:
     """Represents a date and time used for performing astronomy calculations.
@@ -705,7 +705,7 @@ class Time:
         Parameters
         ----------
         year : int
-            The UTC 4-digit year value, e.g. 2019.
+            The UTC year value, e.g. 2019.
         month : int
             The UTC month in the range 1..12.
         day : int
@@ -721,10 +721,17 @@ class Time:
         -------
         Time
         """
-        micro = round(math.fmod(second, 1.0) * 1000000)
-        second = math.floor(second - micro/1000000)
-        d = datetime.datetime(year, month, day, hour, minute, second, micro)
-        ut = (d - _EPOCH).total_seconds() / 86400
+        # This formula is adapted from NOVAS C 3.1 function julian_date().
+        y = int(year)
+        m = int(month)
+        d = int(day)
+        y2000 = float(
+            (d - 2483620)
+            + 1461 * (y + 4800 - (14 - m) // 12) // 4
+            + 367 * (m - 2 + (14 - m) // 12 * 12) // 12
+            - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4
+        )
+        ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0)
         return Time(ut)
 
     @staticmethod
@@ -774,9 +781,38 @@ class Time:
         return 'Time(\'' + str(self) + '\')'
 
     def __str__(self):
-        millis = round(self.ut * 86400000.0)
-        n = _EPOCH + datetime.timedelta(milliseconds=millis)
-        return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000))
+        # Adapted from the NOVAS C 3.1 function cal_date().
+        djd = self.ut + 2451545.5
+        jd = int(djd)
+        x = 24.0 * math.fmod(djd, 1.0)
+        hour = int(x)
+        x = 60.0 * math.fmod(x, 1.0)
+        minute = int(x)
+        second = round(60.0 * math.fmod(x, 1.0), 3)
+        if second >= 60.0:
+            second -= 60.0
+            minute += 1
+            if minute >= 60:
+                minute -= 60
+                hour += 1
+                if hour >= 24:
+                    hour -= 24
+                    jd += 1
+        k = jd + 68569
+        n = 4 * k // 146097
+        k = k - (146097 * n + 3) // 4
+        m = 4000 * (k + 1) // 1461001
+        k = k - 1461 * m // 4 + 31
+        month = (80 * k) // 2447
+        day = k - 2447 * month // 80
+        k = month // 11
+        month = month + 2 - 12 * k
+        year = 100 * (n - 49) + m + k
+        millis = max(0, min(59999, round(1000.0 * second)))
+        text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000)
+        if year > 9999:
+            text = '+' + text
+        return text
 
     def Utc(self):
         """Returns the UTC date and time as a `datetime` object.
diff --git a/demo/python/correct/lunar_angles.txt b/demo/python/correct/lunar_angles.txt
index 865eabf5..80b97880 100644
--- a/demo/python/correct/lunar_angles.txt
+++ b/demo/python/correct/lunar_angles.txt
@@ -1,4 +1,4 @@
-2021-05-15T01:45:15.503Z Jupiter  120
+2021-05-15T01:45:15.502Z Jupiter  120
 2021-05-15T17:53:54.557Z Venus     30
 2021-05-16T04:23:24.323Z Saturn   150
 2021-05-16T05:06:25.856Z Mars       0
@@ -15,7 +15,7 @@
 2021-05-20T23:01:26.198Z Venus     90
 2021-05-21T02:57:43.803Z Mars      60
 2021-05-21T11:47:54.388Z Mercury   90
-2021-05-22T02:58:57.124Z Jupiter  210
+2021-05-22T02:58:57.125Z Jupiter  210
 2021-05-22T03:46:48.175Z Sun      120
 2021-05-23T00:10:42.249Z Saturn   240
 2021-05-23T06:36:37.117Z Venus    120
@@ -53,7 +53,7 @@
 2021-06-03T21:09:24.870Z Venus    270
 2021-06-03T21:22:13.482Z Jupiter   30
 2021-06-04T20:23:38.951Z Saturn    60
-2021-06-04T22:38:31.596Z Sun      300
+2021-06-04T22:38:31.597Z Sun      300
 2021-06-05T15:57:17.165Z Mercury  300
 2021-06-05T22:47:34.963Z Mars     270
 2021-06-06T09:32:04.351Z Jupiter   60
diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c
index e79537f0..35e0fcf7 100644
--- a/generate/template/astronomy.c
+++ b/generate/template/astronomy.c
@@ -910,19 +910,19 @@ astro_time_t Astronomy_CurrentTime(void)
 astro_time_t Astronomy_MakeTime(int year, int month, int day, int hour, int minute, double second)
 {
     astro_time_t time;
-    long int jd12h;
-    long int y2000;
+    long y = (long)year;
+    long m = (long)month;
+    long d = (long)day;
 
     /* This formula is adapted from NOVAS C 3.1 function julian_date() */
-    jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L
-        + ((long) month - 14L) / 12L) / 4L
-        + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L)
-        / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L)
-        / 100L) / 4L;
+    long y2000 = (
+        (d - 2483620L)
+        + 1461L*(y + 4800L - (14L - m) / 12L)/4L
+        + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
+        - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+    );
 
-    y2000 = jd12h - 2451545L;
-
-    time.ut = (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+    time.ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
     time.tt = TerrestrialTime(time.ut);
     time.psi = time.eps = time.st = NAN;
 
diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs
index ef971fae..80881efe 100644
--- a/generate/template/astronomy.cs
+++ b/generate/template/astronomy.cs
@@ -321,15 +321,19 @@ namespace CosineKitty
             // It given a Gregorian calendar date/time, it calculates the fractional
             // number of days since the J2000 epoch.
 
-            long jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L
-                + ((long) month - 14L) / 12L) / 4L
-                + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L)
-                / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L)
-                / 100L) / 4L;
+            long y = (long)year;
+            long m = (long)month;
+            long d = (long)day;
 
-            long y2000 = jd12h - 2451545L;
+            long y2000 = (
+                (d - 2483620L)
+                + 1461L*(y + 4800L - (14L - m) / 12L)/4L
+                + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
+                - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+            );
 
-            return (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+            double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+            return ut;
         }
     }
 
diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index b304fba9..99b0a43a 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -577,7 +577,7 @@ def _UniversalTime(tt):
             return ut
         dt += err
 
-_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
+_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
 
 class Time:
     """Represents a date and time used for performing astronomy calculations.
@@ -705,7 +705,7 @@ class Time:
         Parameters
         ----------
         year : int
-            The UTC 4-digit year value, e.g. 2019.
+            The UTC year value, e.g. 2019.
         month : int
             The UTC month in the range 1..12.
         day : int
@@ -721,10 +721,17 @@ class Time:
         -------
         Time
         """
-        micro = round(math.fmod(second, 1.0) * 1000000)
-        second = math.floor(second - micro/1000000)
-        d = datetime.datetime(year, month, day, hour, minute, second, micro)
-        ut = (d - _EPOCH).total_seconds() / 86400
+        # This formula is adapted from NOVAS C 3.1 function julian_date().
+        y = int(year)
+        m = int(month)
+        d = int(day)
+        y2000 = float(
+            (d - 2483620)
+            + 1461 * (y + 4800 - (14 - m) // 12) // 4
+            + 367 * (m - 2 + (14 - m) // 12 * 12) // 12
+            - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4
+        )
+        ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0)
         return Time(ut)
 
     @staticmethod
@@ -774,9 +781,38 @@ class Time:
         return 'Time(\'' + str(self) + '\')'
 
     def __str__(self):
-        millis = round(self.ut * 86400000.0)
-        n = _EPOCH + datetime.timedelta(milliseconds=millis)
-        return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000))
+        # Adapted from the NOVAS C 3.1 function cal_date().
+        djd = self.ut + 2451545.5
+        jd = int(djd)
+        x = 24.0 * math.fmod(djd, 1.0)
+        hour = int(x)
+        x = 60.0 * math.fmod(x, 1.0)
+        minute = int(x)
+        second = round(60.0 * math.fmod(x, 1.0), 3)
+        if second >= 60.0:
+            second -= 60.0
+            minute += 1
+            if minute >= 60:
+                minute -= 60
+                hour += 1
+                if hour >= 24:
+                    hour -= 24
+                    jd += 1
+        k = jd + 68569
+        n = 4 * k // 146097
+        k = k - (146097 * n + 3) // 4
+        m = 4000 * (k + 1) // 1461001
+        k = k - 1461 * m // 4 + 31
+        month = (80 * k) // 2447
+        day = k - 2447 * month // 80
+        k = month // 11
+        month = month + 2 - 12 * k
+        year = 100 * (n - 49) + m + k
+        millis = max(0, min(59999, round(1000.0 * second)))
+        text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000)
+        if year > 9999:
+            text = '+' + text
+        return text
 
     def Utc(self):
         """Returns the UTC date and time as a `datetime` object.
diff --git a/generate/test.py b/generate/test.py
index 25babef2..89f6cebd 100755
--- a/generate/test.py
+++ b/generate/test.py
@@ -63,25 +63,25 @@ def AstroTime():
     diff = time.ut - expected_ut
     if vabs(diff) > 1.0e-12:
         print('PY AstroTime: excessive UT error {}'.format(diff))
-        sys.exit(1)
+        return 1
     diff = time.tt - expected_tt
     if vabs(diff) > 1.0e-12:
         print('PY AstroTime: excessive TT error {}'.format(diff))
-        sys.exit(1)
+        return 1
     s = str(time.Utc())
     if s != '2018-12-02 18:30:12.543000':
         print('PY AstroTime: Utc() returned incorrect string "{}"'.format(s))
-        sys.exit(1)
+        return 1
     time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9994)
     s = str(time)
     if s != '2018-12-31T23:59:59.999Z':
         print('PY AstroTime: expected 2018-12-31T23:59:59.999Z but found {}'.format(s))
-        sys.exit(1)
+        return 1
     time = astronomy.Time.Make(2018, 12, 31, 23, 59, 59.9995)
     s = str(time)
     if s != '2019-01-01T00:00:00.000Z':
         print('PY AstroTime: expected 2019-01-01T00:00:00.000Z but found {}'.format(s))
-        sys.exit(1)
+        return 1
     print('PY Current time =', astronomy.Time.Now())
     AssertGoodTime('2015-12-31', '2015-12-31T00:00:00.000Z')
     AssertGoodTime('2015-12-31T23:45Z', '2015-12-31T23:45:00.000Z')
@@ -94,6 +94,9 @@ def AstroTime():
     AssertBadTime('1971-12-31T23:60:00Z')
     AssertBadTime('1971-12-31T23:00:60Z')
     AssertBadTime('1971-03-17T03:30:55.976')
+    # Extreme year values...
+    AssertGoodTime('-2300-12-19T16:22:26.325Z', '-2300-12-19T16:22:26.325Z')
+    AssertGoodTime('+12345-12-11T13:30:10.041Z', '+12345-12-11T13:30:10.041Z')
     return 0
 
 #-----------------------------------------------------------------------------------------------------------
@@ -2851,11 +2854,32 @@ def ArcminVelError(correct, calc):
 
 #-----------------------------------------------------------------------------------------------------------
 
+def CheckDecemberSolstice(year, expected):
+    si = astronomy.Seasons(year)
+    actual = str(si.dec_solstice)
+    if actual != expected:
+        print('PY DatesIssue250: FAIL: year {}, expected [{}], actual [{}]'.format(year, expected, actual))
+        return 1
+    return 0
+
+def DatesIssue250():
+    # Make sure we can handle dates outside the range supported by System.DateTime.
+    # https://github.com/cosinekitty/astronomy/issues/250
+    return (
+        CheckDecemberSolstice( 2022, "2022-12-21T21:47:58.189Z") or
+        CheckDecemberSolstice(-2300, "-2300-12-19T16:22:26.325Z") or
+        CheckDecemberSolstice(12345, "+12345-12-11T13:30:10.041Z") or
+        Pass('DatesIssue250')
+    )
+
+#-----------------------------------------------------------------------------------------------------------
+
 UnitTests = {
     'aberration':               Aberration,
     'axis':                     Axis,
     'barystate':                BaryState,
     'constellation':            Constellation,
+    'dates250':                 DatesIssue250,
     'elongation':               Elongation,
     'geoid':                    Geoid,
     'global_solar_eclipse':     GlobalSolarEclipse,
diff --git a/source/c/astronomy.c b/source/c/astronomy.c
index a547f5af..4fec4cbc 100644
--- a/source/c/astronomy.c
+++ b/source/c/astronomy.c
@@ -916,19 +916,19 @@ astro_time_t Astronomy_CurrentTime(void)
 astro_time_t Astronomy_MakeTime(int year, int month, int day, int hour, int minute, double second)
 {
     astro_time_t time;
-    long int jd12h;
-    long int y2000;
+    long y = (long)year;
+    long m = (long)month;
+    long d = (long)day;
 
     /* This formula is adapted from NOVAS C 3.1 function julian_date() */
-    jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L
-        + ((long) month - 14L) / 12L) / 4L
-        + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L)
-        / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L)
-        / 100L) / 4L;
+    long y2000 = (
+        (d - 2483620L)
+        + 1461L*(y + 4800L - (14L - m) / 12L)/4L
+        + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
+        - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+    );
 
-    y2000 = jd12h - 2451545L;
-
-    time.ut = (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+    time.ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
     time.tt = TerrestrialTime(time.ut);
     time.psi = time.eps = time.st = NAN;
 
diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs
index d5222c21..32f42faa 100644
--- a/source/csharp/astronomy.cs
+++ b/source/csharp/astronomy.cs
@@ -321,15 +321,19 @@ namespace CosineKitty
             // It given a Gregorian calendar date/time, it calculates the fractional
             // number of days since the J2000 epoch.
 
-            long jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L
-                + ((long) month - 14L) / 12L) / 4L
-                + 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L)
-                / 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L)
-                / 100L) / 4L;
+            long y = (long)year;
+            long m = (long)month;
+            long d = (long)day;
 
-            long y2000 = jd12h - 2451545L;
+            long y2000 = (
+                (d - 2483620L)
+                + 1461L*(y + 4800L - (14L - m) / 12L)/4L
+                + 367L*(m - 2L + (14L - m) / 12L * 12L)/12L
+                - 3L*((y + 4900L - (14L - m) / 12L) / 100L)/4L
+            );
 
-            return (double)y2000 - 0.5 + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+            double ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
+            return ut;
         }
     }
 
diff --git a/source/python/README.md b/source/python/README.md
index ea6477f2..eecb0ac8 100644
--- a/source/python/README.md
+++ b/source/python/README.md
@@ -917,7 +917,7 @@ The value of the calling object is not modified. This function creates a brand n
 
 | Type | Parameter | Description |
 | --- | --- | --- |
-| `int` | `year` | The UTC 4-digit year value, e.g. 2019. |
+| `int` | `year` | The UTC year value, e.g. 2019. |
 | `int` | `month` | The UTC month in the range 1..12. |
 | `int` | `day` | The UTC day of the month, in the range 1..31. |
 | `int` | `hour` | The UTC hour, in the range 0..23. |
diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py
index 9a74b0cb..7c4650b3 100644
--- a/source/python/astronomy/astronomy.py
+++ b/source/python/astronomy/astronomy.py
@@ -577,7 +577,7 @@ def _UniversalTime(tt):
             return ut
         dt += err
 
-_TimeRegex = re.compile(r'^([0-9]{1,4})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
+_TimeRegex = re.compile(r'^([\+\-]?[0-9]{1,5})-([0-9]{2})-([0-9]{2})(T([0-9]{2}):([0-9]{2})(:([0-9]{2}(\.[0-9]+)?))?Z)?$')
 
 class Time:
     """Represents a date and time used for performing astronomy calculations.
@@ -705,7 +705,7 @@ class Time:
         Parameters
         ----------
         year : int
-            The UTC 4-digit year value, e.g. 2019.
+            The UTC year value, e.g. 2019.
         month : int
             The UTC month in the range 1..12.
         day : int
@@ -721,10 +721,17 @@ class Time:
         -------
         Time
         """
-        micro = round(math.fmod(second, 1.0) * 1000000)
-        second = math.floor(second - micro/1000000)
-        d = datetime.datetime(year, month, day, hour, minute, second, micro)
-        ut = (d - _EPOCH).total_seconds() / 86400
+        # This formula is adapted from NOVAS C 3.1 function julian_date().
+        y = int(year)
+        m = int(month)
+        d = int(day)
+        y2000 = float(
+            (d - 2483620)
+            + 1461 * (y + 4800 - (14 - m) // 12) // 4
+            + 367 * (m - 2 + (14 - m) // 12 * 12) // 12
+            - 3 * ((y + 4900 - (14 - m) // 12) // 100) // 4
+        )
+        ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0)
         return Time(ut)
 
     @staticmethod
@@ -774,9 +781,38 @@ class Time:
         return 'Time(\'' + str(self) + '\')'
 
     def __str__(self):
-        millis = round(self.ut * 86400000.0)
-        n = _EPOCH + datetime.timedelta(milliseconds=millis)
-        return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000))
+        # Adapted from the NOVAS C 3.1 function cal_date().
+        djd = self.ut + 2451545.5
+        jd = int(djd)
+        x = 24.0 * math.fmod(djd, 1.0)
+        hour = int(x)
+        x = 60.0 * math.fmod(x, 1.0)
+        minute = int(x)
+        second = round(60.0 * math.fmod(x, 1.0), 3)
+        if second >= 60.0:
+            second -= 60.0
+            minute += 1
+            if minute >= 60:
+                minute -= 60
+                hour += 1
+                if hour >= 24:
+                    hour -= 24
+                    jd += 1
+        k = jd + 68569
+        n = 4 * k // 146097
+        k = k - (146097 * n + 3) // 4
+        m = 4000 * (k + 1) // 1461001
+        k = k - 1461 * m // 4 + 31
+        month = (80 * k) // 2447
+        day = k - 2447 * month // 80
+        k = month // 11
+        month = month + 2 - 12 * k
+        year = 100 * (n - 49) + m + k
+        millis = max(0, min(59999, round(1000.0 * second)))
+        text = '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(year, month, day, hour, minute, millis // 1000, millis % 1000)
+        if year > 9999:
+            text = '+' + text
+        return text
 
     def Utc(self):
         """Returns the UTC date and time as a `datetime` object.