ENH: datetime: Got datetime interunit comparisons working

6 files changed, 581 insertions, 244 deletions

diff --git a/numpy/core/include/numpy/ndarraytypes.h b/numpy/core/include/numpy/ndarraytypes.h
index 0b07ad762..3e6af3c08 100644
--- a/numpy/core/include/numpy/ndarraytypes.h
+++ b/numpy/core/include/numpy/ndarraytypes.h
@@ -213,6 +213,9 @@ typedef enum {
         NPY_RAISE=2
 } NPY_CLIPMODE;
 
+/* The special not-a-time (NaT) value */
+#define NPY_DATETIME_NAT NPY_MIN_INT64
+
 typedef enum {
         NPY_FR_Y, /* Years */
         NPY_FR_M, /* Months */
diff --git a/numpy/core/src/multiarray/_datetime.h b/numpy/core/src/multiarray/_datetime.h
index ed80ba2de..3a7e3426d 100644
--- a/numpy/core/src/multiarray/_datetime.h
+++ b/numpy/core/src/multiarray/_datetime.h
@@ -89,7 +89,24 @@ convert_datetime_divisor_to_multiple(PyArray_DatetimeMetaData *meta,
 NPY_NO_EXPORT PyObject *
 compute_datetime_metadata_greatest_common_divisor(
                         PyArray_Descr *type1,
-                        PyArray_Descr *type2);
+                        PyArray_Descr *type2,
+                        int strict_with_nonlinear_units);
+
+/*
+ * Computes the conversion factor to convert data with 'src_meta' metadata
+ * into data with 'dst_meta' metadata, not taking into account the events.
+ *
+ * To convert a npy_datetime or npy_timedelta, first the event number needs
+ * to be divided away, then it needs to be scaled by num/denom, and
+ * finally the event number can be added back in.
+ *
+ * If overflow occurs, both out_num and out_denom are set to 0, but
+ * no error is set.
+ */
+NPY_NO_EXPORT void
+get_datetime_conversion_factor(PyArray_DatetimeMetaData *src_meta,
+                                PyArray_DatetimeMetaData *dst_meta,
+                                npy_int64 *out_num, npy_int64 *out_denom);
 
 /*
  * Given an the capsule datetime metadata object,
@@ -160,4 +177,25 @@ convert_pyobject_to_datetime(PyArray_DatetimeMetaData *meta, PyObject *obj,
 NPY_NO_EXPORT PyObject *
 convert_datetime_to_pyobject(npy_datetime dt, PyArray_DatetimeMetaData *meta);
 
+/*
+ * Converts a datetime based on the given metadata into a datetimestruct
+ */
+NPY_NO_EXPORT int
+convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
+                                    npy_datetime dt,
+                                    npy_datetimestruct *out);
+
+/*
+ * Converts a datetime from a datetimestruct to a datetime based
+ * on some metadata. The date is assumed to be valid.
+ *
+ * TODO: If meta->num is really big, there could be overflow
+ *
+ * Returns 0 on success, -1 on failure.
+ */
+NPY_NO_EXPORT int
+convert_datetimestruct_to_datetime(PyArray_DatetimeMetaData *meta,
+                                    const npy_datetimestruct *dts,
+                                    npy_datetime *out);
+
 #endif
diff --git a/numpy/core/src/multiarray/datetime.c b/numpy/core/src/multiarray/datetime.c
index 7adb24d1e..f3fa0534c 100644
--- a/numpy/core/src/multiarray/datetime.c
+++ b/numpy/core/src/multiarray/datetime.c
@@ -23,6 +23,9 @@ numpy_pydatetime_import()
     PyDateTime_IMPORT;
 }
 
+static int
+is_leapyear(npy_int64 year);
+
 
 /* For defaults and errors */
 #define NPY_FR_ERR  -1
@@ -33,10 +36,6 @@ numpy_pydatetime_import()
 
 /* Calendar Structure for Parsing Long -> Date */
 typedef struct {
-    int year, month, day;
-} ymdstruct;
-
-typedef struct {
     int hour, min, sec;
 } hmsstruct;
 
@@ -78,26 +77,12 @@ NPY_NO_EXPORT char *_datetime_strings[] = {
         } \
     }
 
-/* Table with day offsets for each month (0-based, without and with leap) */
-static int month_offset[2][13] = {
-    { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
-    { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
-};
-
 /* Table of number of days in a month (0-based, without and with leap) */
 static int days_in_month[2][12] = {
     { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 },
     { 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }
 };
 
-/* Return 1/0 iff year points to a leap year in calendar. */
-static int
-is_leapyear(npy_int64 year)
-{
-    return (year % 4 == 0) && ((year % 100 != 0) || (year % 400 == 0));
-}
-
-
 /*
  * Return the day of the week for the given absolute date.
  * Monday is 0 and Sunday is 6
@@ -116,186 +101,180 @@ day_of_week(npy_longlong absdate)
     }
 }
 
-/*
- * Return the year offset, that is the absolute date of the day
- * 31.12.(year-1) since 31.12.1969 in the proleptic Gregorian calendar.
+/* Returns absolute seconds from an hour, minute, and second
  */
-static npy_longlong
-year_offset(npy_longlong year)
-{
-    /* Note that 477 == 1969/4 - 1969/100 + 1969/400 */
-    year--;
-    if (year >= 0 || -1/4 == -1)
-        return (year-1969)*365 + year/4 - year/100 + year/400 - 477;
-    else
-        return (year-1969)*365 + (year-3)/4 - (year-99)/100 + (year-399)/400 - 477;
-}
+#define secs_from_hms(hour, min, sec, multiplier) (\
+  ((hour)*3600 + (min)*60 + (sec)) * (npy_int64)(multiplier)\
+)
 
 /*
- * Modified version of mxDateTime function
- * Returns absolute number of days since Jan 1, 1970
- * assuming a proleptic Gregorian Calendar
- * Raises a ValueError if out of range month or day
- * day -1 is Dec 31, 1969, day 0 is Jan 1, 1970, day 1 is Jan 2, 1970
+ * Converts an integer number of seconds in a day to hours minutes seconds.
+ * It assumes seconds is between 0 and 86399.
  */
-static npy_longlong
-days_from_ymd(int year, int month, int day)
-{
-
-    /* Calculate the absolute date */
-    int leap;
-    npy_longlong yearoffset, absdate;
-
-    /* Is it a leap year ? */
-    leap = is_leapyear(year);
 
-    /* Negative month values indicate months relative to the years end */
-    if (month < 0) month += 13;
-    Py_AssertWithArg(month >= 1 && month <= 12,
-                     PyExc_ValueError,
-                     "month out of range (1-12): %i",
-                     month);
-
-    /* Negative values indicate days relative to the months end */
-    if (day < 0) day += days_in_month[leap][month - 1] + 1;
-    Py_AssertWithArg(day >= 1 && day <= days_in_month[leap][month - 1],
-                     PyExc_ValueError,
-                     "day out of range: %i",
-                     day);
+static hmsstruct
+seconds_to_hmsstruct(npy_longlong dlong)
+{
+    int hour, minute, second;
+    hmsstruct hms;
 
-    /*
-     * Number of days between Dec 31, (year - 1) and Dec 31, 1969
-     *    (can be negative).
-     */
-    yearoffset = year_offset(year);
+    hour   = dlong / 3600;
+    minute = (dlong % 3600) / 60;
+    second = dlong - (hour*3600 + minute*60);
 
-    if (PyErr_Occurred()) goto onError;
+    hms.hour   = hour;
+    hms.min = minute;
+    hms.sec = second;
 
-    /*
-     * Calculate the number of days using yearoffset
-     * Jan 1, 1970 is day 0 and thus Dec. 31, 1969 is day -1
-     */
-    absdate = day-1 + month_offset[leap][month - 1] + yearoffset;
+    return hms;
+}
 
-    return absdate;
+/*
+  ====================================================
+  == End of section adapted from mx.DateTime       ==
+  ====================================================
+*/
 
- onError:
-    return 0;
 
+static int
+is_leapyear(npy_int64 year)
+{
+    return (year & 0x3) == 0 && /* year % 4 == 0 */
+           ((year % 100) != 0 ||
+            (year % 400) == 0);
 }
 
-/* Returns absolute seconds from an hour, minute, and second
- */
-#define secs_from_hms(hour, min, sec, multiplier) (\
-  ((hour)*3600 + (min)*60 + (sec)) * (npy_int64)(multiplier)\
-)
-
 /*
- * Takes a number of days since Jan 1, 1970 (positive or negative)
- * and returns the year. month, and day in the proleptic
- * Gregorian calendar
- *
- * Examples:
- *
- * -1 returns 1969, 12, 31
- * 0  returns 1970, 1, 1
- * 1  returns 1970, 1, 2
+ * Calculates the days offset from the 1970 epoch.
  */
-
-static ymdstruct
-days_to_ymdstruct(npy_datetime dlong)
+static npy_int64
+get_datetimestruct_days(const npy_datetimestruct *dts)
 {
-    ymdstruct ymd;
-    long year;
-    npy_longlong yearoffset;
-    int leap, dayoffset;
-    int month = 1, day = 1;
-    int *monthoffset;
+    int i, month;
+    npy_int64 year, days = 0;
+    int *month_lengths;
 
-    dlong += 1;
-
-    /* Approximate year */
-    year = 1970 + dlong / 365.2425;
-
-    /* Apply corrections to reach the correct year */
-    while (1) {
-        /* Calculate the year offset */
-        yearoffset = year_offset(year);
+    year = dts->year - 1970;
+    days = year * 365;
 
+    /* Adjust for leap years */
+    if (days >= 0) {
         /*
-         * Backward correction: absdate must be greater than the
-         * yearoffset
+         * 1968 is the closest leap year before 1970.
+         * Exclude the current year, so add 1.
          */
-        if (yearoffset >= dlong) {
-            year--;
-            continue;
-        }
-
-        dayoffset = dlong - yearoffset;
-        leap = is_leapyear(year);
-
-        /* Forward correction: non leap years only have 365 days */
-        if (dayoffset > 365 && !leap) {
-            year++;
-            continue;
-        }
-        break;
+        year += 1;
+        /* Add one day for each 4 years */
+        days += year / 4;
+        /* 1900 is the closest previous year divisible by 100 */
+        year += 68;
+        /* Subtract one day for each 100 years */
+        days -= year / 100;
+        /* 1600 is the closest previous year divisible by 400 */
+        year += 300;
+        /* Add one day for each 400 years */
+        days += year / 400;
+    }
+    else {
+        /*
+         * 1972 is the closest later year after 1970.
+         * Include the current year, so subtract 2.
+         */
+        year -= 2;
+        /* Subtract one day for each 4 years */
+        days += year / 4;
+        /* 2000 is the closest later year divisible by 100 */
+        year -= 28;
+        /* Add one day for each 100 years */
+        days -= year / 100;
+        /* 2000 is also the closest later year divisible by 400 */
+        /* Subtract one day for each 400 years */
+        days += year / 400;
     }
 
-    /* Now iterate to find the month */
-    monthoffset = month_offset[leap];
-    for (month = 1; month < 13; month++) {
-        if (monthoffset[month] >= dayoffset)
-            break;
+    month_lengths = days_in_month[is_leapyear(dts->year)];
+    month = dts->month - 1;
+
+    /* Add the months */
+    for (i = 0; i < month; ++i) {
+        days += month_lengths[i];
     }
-    day = dayoffset - month_offset[leap][month-1];
 
-    ymd.year  = year;
-    ymd.month = month;
-    ymd.day   = day;
+    /* Add the days */
+    days += dts->day - 1;
 
-    return ymd;
+    return days;
 }
 
 /*
- * Converts an integer number of seconds in a day to hours minutes seconds.
- * It assumes seconds is between 0 and 86399.
+ * Modifies '*days_' to be the day offset within the year,
+ * and returns the year.
  */
-
-static hmsstruct
-seconds_to_hmsstruct(npy_longlong dlong)
+static npy_int64
+days_to_yearsdays(npy_int64 *days_)
 {
-    int hour, minute, second;
-    hmsstruct hms;
-
-    hour   = dlong / 3600;
-    minute = (dlong % 3600) / 60;
-    second = dlong - (hour*3600 + minute*60);
+    const npy_int64 days_per_400years = (400*365 + 100 - 4 + 1);
+    /* Adjust so it's relative to the year 2000 (divisible by 400) */
+    npy_int64 days = (*days_) - (365*30 + 7), year;
 
-    hms.hour   = hour;
-    hms.min = minute;
-    hms.sec = second;
+    /* Break down the 400 year cycle to get the year and day within the year */
+    if (days >= 0) {
+        year = 400 * (days / days_per_400years);
+        days = days % days_per_400years;
+    }
+    else {
+        year = 400 * ((days - (days_per_400years - 1)) / days_per_400years);
+        days = days % days_per_400years;
+        if (days < 0) {
+            days += days_per_400years;
+        }
+    }
+
+    /* Work out the year/day within the 400 year cycle */
+    if (days >= 366) {
+        year += 100 * ((days-1) / (100*365 + 25 - 1));
+        days = (days-1) % (100*365 + 25 - 1);
+        if (days >= 365) {
+            year += 4 * ((days+1) / (4*365 + 1));
+            days = (days+1) % (4*365 + 1);
+            if (days >= 366) {
+                year += (days-1) / 365;
+                days = (days-1) % 365;
+            }
+        }
+    }
 
-    return hms;
+    *days_ = days;
+    return year + 2000;
 }
 
 /*
-  ====================================================
-  == End of section adapted from mx.DateTime       ==
-  ====================================================
-*/
-
+ * Fills in the year, month, day in 'dts' based on the days
+ * offset from 1970.
+ */
+static void
+set_datetimestruct_days(npy_int64 days, npy_datetimestruct *dts)
+{
+    int *month_lengths, i;
 
-/*==================================================
-// Parsing DateTime struct and returns a date-time number
-// =================================================
+    dts->year = days_to_yearsdays(&days);
 
- Structure is assumed to be already normalized
-*/
+    month_lengths = days_in_month[is_leapyear(dts->year)];
+    for (i = 0; i < 12; ++i) {
+        if (days < month_lengths[i]) {
+            dts->month = i + 1;
+            dts->day = days + 1;
+            return;
+        }
+        else {
+            days -= month_lengths[i];
+        }
+    }
+}
 
 /*
  * Converts a datetime from a datetimestruct to a datetime based
- * on some metadata.
+ * on some metadata. The date is assumed to be valid.
  *
  * TODO: If meta->num is really big, there could be overflow
  *
@@ -310,8 +289,8 @@ convert_datetimestruct_to_datetime(PyArray_DatetimeMetaData *meta,
     NPY_DATETIMEUNIT base = meta->base;
 
     /* If the datetimestruct is NaT, return NaT */
-    if (dts->year == NPY_MIN_INT64) {
-        *out = NPY_MIN_INT64;
+    if (dts->year == NPY_DATETIME_NAT) {
+        *out = NPY_DATETIME_NAT;
         return 0;
     }
 
@@ -332,7 +311,7 @@ convert_datetimestruct_to_datetime(PyArray_DatetimeMetaData *meta,
     }
     else {
         /* Otherwise calculate the number of days to start */
-        npy_int64 days = days_from_ymd(dts->year, dts->month, dts->day);
+        npy_int64 days = get_datetimestruct_days(dts);
 
         if (base == NPY_FR_W) {
             /* Truncate to weeks */
@@ -595,7 +574,6 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
                                     npy_datetime dt,
                                     npy_datetimestruct *out)
 {
-    ymdstruct ymd;
     hmsstruct hms;
     npy_int64 absdays;
     npy_int64 tmp, num1, num2, num3;
@@ -641,10 +619,7 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
 
         case NPY_FR_W:
             /* A week is 7 days */
-            ymd = days_to_ymdstruct(dt * 7);
-            out->year  = ymd.year;
-            out->month = ymd.month;
-            out->day   = ymd.day;
+            set_datetimestruct_days(dt * 7, out);
             break;
 
         case NPY_FR_B:
@@ -666,63 +641,48 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
                 /* Recall how C computes / and % with negative numbers */
                 absdays = 7 * ((dt - 1) / 5) + ((dt - 1) % 5) + 1;
             }
-            ymd = days_to_ymdstruct(absdays);
-            out->year  = ymd.year;
-            out->month = ymd.month;
-            out->day   = ymd.day;
+            set_datetimestruct_days(absdays, out);
             break;
 
         case NPY_FR_D:
-            ymd = days_to_ymdstruct(dt);
-            out->year  = ymd.year;
-            out->month = ymd.month;
-            out->day   = ymd.day;
+            set_datetimestruct_days(dt, out);
             break;
 
         case NPY_FR_h:
             if (dt >= 0) {
-                ymd  = days_to_ymdstruct(dt / 24);
+                set_datetimestruct_days(dt / 24, out);
                 out->hour  = dt % 24;
             }
             else {
-                ymd  = days_to_ymdstruct((dt - 23) / 24);
+                set_datetimestruct_days((dt - 23) / 24, out);
                 out->hour = 23 + (dt + 1) % 24;
             }
-            out->year  = ymd.year;
-            out->month = ymd.month;
-            out->day   = ymd.day;
             break;
 
         case NPY_FR_m:
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / 1440);
+                set_datetimestruct_days(dt / 1440, out);
                 out->min = dt % 1440;
             }
             else {
-                ymd = days_to_ymdstruct((dt - 1439) / 1440);
+                set_datetimestruct_days((dt - 1439) / 1440, out);
                 out->min = 1439 + (dt + 1) % 1440;
             }
             hms = seconds_to_hmsstruct(out->min * 60);
-            out->year   = ymd.year;
-            out->month  = ymd.month;
-            out->day    = ymd.day;
             out->hour   = hms.hour;
             out->min = hms.min;
             break;
 
         case NPY_FR_s:
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / 86400);
+                set_datetimestruct_days(dt / 86400, out);
                 out->sec = dt % 86400;
             }
             else {
-                ymd = days_to_ymdstruct((dt - 86399) / 86400);
+                set_datetimestruct_days((dt - 86399) / 86400, out);
                 out->sec = 86399 + (dt + 1) % 86400;
             }
             hms = seconds_to_hmsstruct(out->sec);
-            out->year   = ymd.year;
-            out->month  = ymd.month;
-            out->day    = ymd.day;
             out->hour   = hms.hour;
             out->min = hms.min;
             out->sec = hms.sec;
@@ -730,18 +690,15 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
 
         case NPY_FR_ms:
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / 86400000);
+                set_datetimestruct_days(dt / 86400000, out);
                 tmp  = dt % 86400000;
             }
             else {
-                ymd = days_to_ymdstruct((dt - 86399999) / 86400000);
+                set_datetimestruct_days((dt - 86399999) / 86400000, out);
                 tmp  = 86399999 + (dt + 1) % 86399999;
             }
             hms = seconds_to_hmsstruct(tmp / 1000);
             out->us  = (tmp % 1000)*1000;
-            out->year    = ymd.year;
-            out->month   = ymd.month;
-            out->day     = ymd.day;
             out->hour    = hms.hour;
             out->min     = hms.min;
             out->sec     = hms.sec;
@@ -752,18 +709,15 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
             num1 *= 1000;
             num2 = num1 - 1;
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / num1);
+                set_datetimestruct_days(dt / num1, out);
                 tmp = dt % num1;
             }
             else {
-                ymd = days_to_ymdstruct((dt - num2)/ num1);
+                set_datetimestruct_days((dt - num2)/ num1, out);
                 tmp = num2 + (dt + 1) % num1;
             }
             hms = seconds_to_hmsstruct(tmp / 1000000);
             out->us = tmp % 1000000;
-            out->year    = ymd.year;
-            out->month   = ymd.month;
-            out->day     = ymd.day;
             out->hour    = hms.hour;
             out->min     = hms.min;
             out->sec     = hms.sec;
@@ -776,20 +730,17 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
             num3 = 1000000;
             num3 *= 1000000;
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / num1);
+                set_datetimestruct_days(dt / num1, out);
                 tmp = dt % num1;
             }
             else {
-                ymd = days_to_ymdstruct((dt - num2)/ num1);
+                set_datetimestruct_days((dt - num2)/ num1, out);
                 tmp = num2 + (dt + 1) % num1;
             }
             hms = seconds_to_hmsstruct(tmp / 1000000000);
             tmp = tmp % 1000000000;
             out->us = tmp / 1000;
             out->ps = (tmp % 1000) * (npy_int64)(1000);
-            out->year    = ymd.year;
-            out->month   = ymd.month;
-            out->day     = ymd.day;
             out->hour    = hms.hour;
             out->min     = hms.min;
             out->sec     = hms.sec;
@@ -802,20 +753,17 @@ convert_datetime_to_datetimestruct(PyArray_DatetimeMetaData *meta,
             num2 = num1 - 1;
 
             if (dt >= 0) {
-                ymd = days_to_ymdstruct(dt / num1);
+                set_datetimestruct_days(dt / num1, out);
                 tmp = dt % num1;
             }
             else {
-                ymd = days_to_ymdstruct((dt - num2) / num1);
+                set_datetimestruct_days((dt - num2) / num1, out);
                 tmp = num2 + (dt + 1) % num1;
             }
             hms = seconds_to_hmsstruct(tmp / num3);
             tmp = tmp % num3;
             out->us = tmp / 1000000;
             out->ps = tmp % 1000000;
-            out->year    = ymd.year;
-            out->month   = ymd.month;
-            out->day     = ymd.day;
             out->hour    = hms.hour;
             out->min     = hms.min;
             out->sec     = hms.sec;
@@ -1460,6 +1408,9 @@ _uint64_euclidean_gcd(npy_uint64 x, npy_uint64 y)
  * To convert a npy_datetime or npy_timedelta, first the event number needs
  * to be divided away, then it needs to be scaled by num/denom, and
  * finally the event number can be added back in.
+ *
+ * If overflow occurs, both out_num and out_denom are set to 0, but
+ * no error is set.
  */
 NPY_NO_EXPORT void
 get_datetime_conversion_factor(PyArray_DatetimeMetaData *src_meta,
@@ -1498,7 +1449,7 @@ get_datetime_conversion_factor(PyArray_DatetimeMetaData *src_meta,
                 num *= (97 + 400*365);
                 denom *= 400;
                 /* Day -> dst_base */
-                denom *= get_datetime_units_factor(NPY_FR_D, dst_base);
+                num *= get_datetime_units_factor(NPY_FR_D, dst_base);
             }
         }
         else if (src_base == NPY_FR_M) {
@@ -1511,11 +1462,11 @@ get_datetime_conversion_factor(PyArray_DatetimeMetaData *src_meta,
                 num *= (97 + 400*365);
                 denom *= 400*12;
                 /* Day -> dst_base */
-                denom *= get_datetime_units_factor(NPY_FR_D, dst_base);
+                num *= get_datetime_units_factor(NPY_FR_D, dst_base);
             }
         }
         else {
-            denom *= get_datetime_units_factor(src_base, dst_base);
+            num *= get_datetime_units_factor(src_base, dst_base);
         }
     }
 
@@ -1546,7 +1497,8 @@ get_datetime_conversion_factor(PyArray_DatetimeMetaData *src_meta,
 NPY_NO_EXPORT PyObject *
 compute_datetime_metadata_greatest_common_divisor(
                         PyArray_Descr *type1,
-                        PyArray_Descr *type2)
+                        PyArray_Descr *type2,
+                        int strict_with_nonlinear_units)
 {
     PyArray_DatetimeMetaData *meta1, *meta2, *dt_data;
     NPY_DATETIMEUNIT base;
@@ -1572,15 +1524,12 @@ compute_datetime_metadata_greatest_common_divisor(
         return NULL;
     }
 
-    if (meta1->events != 1 || meta2->events != 1) {
-        /*
-         * When there are events specified, the metadata must
-         * match exactly.
-         */
-        if (meta1->base != meta2->base || meta1->events != meta2->events
-                || meta1->num != meta2->num) {
-            goto incompatible_units;
-        }
+    /* Take the maximum of the events */
+    if (meta1->events > meta2->events) {
+        events = meta1->events;
+    }
+    else {
+        events = meta2->events;
     }
 
     num1 = (npy_uint64)meta1->num;
@@ -1601,24 +1550,50 @@ compute_datetime_metadata_greatest_common_divisor(
                 base = NPY_FR_M;
                 num1 *= 12;
             }
-            else {
+            else if (strict_with_nonlinear_units) {
                 goto incompatible_units;
             }
+            else {
+                base = meta2->base;
+                /* Don't multiply num1 since there is no even factor */
+            }
         }
         else if (meta2->base == NPY_FR_Y) {
             if (meta1->base == NPY_FR_M) {
                 base = NPY_FR_M;
                 num2 *= 12;
             }
-            else {
+            else if (strict_with_nonlinear_units) {
                 goto incompatible_units;
             }
+            else {
+                base = meta1->base;
+                /* Don't multiply num2 since there is no even factor */
+            }
         }
         else if (meta1->base == NPY_FR_M ||
                             meta1->base == NPY_FR_B ||
                             meta2->base == NPY_FR_M ||
                             meta2->base == NPY_FR_B) {
-            goto incompatible_units;
+            if (strict_with_nonlinear_units) {
+                goto incompatible_units;
+            }
+            else {
+                if (meta1->base > meta2->base) {
+                    base = meta1->base;
+                }
+                else {
+                    base = meta2->base;
+                }
+
+                /*
+                 * When combining business days with other units, end
+                 * up with days instead of business days.
+                 */
+                if (base == NPY_FR_B) {
+                    base = NPY_FR_D;
+                }
+            }
         }
 
         /* Take the greater base (unit sizes are decreasing in enum) */
@@ -1668,7 +1643,7 @@ incompatible_units: {
                 PyObject_Repr((PyObject *)type2));
         PyUString_ConcatAndDel(&errmsg,
                 PyUString_FromString(" because they have "
-                    "incompatible base units or events"));
+                    "incompatible nonlinear base time units"));
         PyErr_SetObject(PyExc_TypeError, errmsg);
         return NULL;
     }
@@ -1697,9 +1672,13 @@ datetime_gcd_type_promotion(PyArray_Descr *type1, PyArray_Descr *type2)
     PyObject *gcdmeta;
     PyArray_Descr *dtype;
 
-    /* Get the metadata GCD */
+    /*
+     * Get the metadata GCD, being strict about nonlinear units for
+     * timedelta and relaxed for datetime.
+     */
     gcdmeta = compute_datetime_metadata_greatest_common_divisor(
-                                                type1, type2);
+                                            type1, type2,
+                                            type1->type_num == NPY_TIMEDELTA);
     if (gcdmeta == NULL) {
         return NULL;
     }
@@ -2065,7 +2044,7 @@ parse_iso_8601_date(char *str, int len, npy_datetimestruct *out)
                         tolower(str[0]) == 'n' &&
                         tolower(str[1]) == 'a' &&
                         tolower(str[2]) == 't')) {
-        out->year = NPY_MIN_INT64;
+        out->year = NPY_DATETIME_NAT;
         return 0;
     }
 
@@ -2831,7 +2810,7 @@ convert_datetime_to_pyobject(npy_datetime dt, PyArray_DatetimeMetaData *meta)
     npy_datetimestruct dts;
 
     /* Handle not-a-time */
-    if (dt == NPY_MIN_INT64) {
+    if (dt == NPY_DATETIME_NAT) {
         return PyUString_FromString("NaT");
     }
 
diff --git a/numpy/core/src/multiarray/dtype_transfer.c b/numpy/core/src/multiarray/dtype_transfer.c
index 0162bcc9e..73bb76ac1 100644
--- a/numpy/core/src/multiarray/dtype_transfer.c
+++ b/numpy/core/src/multiarray/dtype_transfer.c
@@ -19,6 +19,10 @@
 #include <numpy/ufuncobject.h>
 #include <numpy/npy_cpu.h>
 
+#include "numpy/npy_3kcompat.h"
+
+#include "_datetime.h"
+
 #include "lowlevel_strided_loops.h"
 
 #define NPY_LOWLEVEL_BUFFER_BLOCKSIZE  128
@@ -696,6 +700,241 @@ get_nbo_cast_numeric_transfer_function(int aligned,
     return NPY_SUCCEED;
 }
 
+/* Does a datetime->datetime or timedelta->timedelta cast */
+typedef struct {
+    free_strided_transfer_data freefunc;
+    copy_strided_transfer_data copyfunc;
+    /* The conversion fraction */
+    npy_int64 num, denom;
+    /* The number of events in the source and destination */
+    int src_events, dst_events;
+    /*
+     * The metadata for when dealing with Months, Years, or
+     * Business Days (all of which behave non-linearly).
+     */
+    PyArray_DatetimeMetaData src_meta, dst_meta;
+} _strided_datetime_cast_data;
+
+/* strided cast data copy function */
+void *_strided_datetime_cast_data_copy(void *data)
+{
+    _strided_datetime_cast_data *newdata =
+            (_strided_datetime_cast_data *)PyArray_malloc(
+                                        sizeof(_strided_datetime_cast_data));
+    if (newdata == NULL) {
+        return NULL;
+    }
+
+    memcpy(newdata, data, sizeof(_strided_datetime_cast_data));
+
+    return (void *)newdata;
+}
+
+static void
+_strided_to_strided_datetime_general_cast(char *dst, npy_intp dst_stride,
+                        char *src, npy_intp src_stride,
+                        npy_intp N, npy_intp src_itemsize,
+                        void *data)
+{
+    _strided_datetime_cast_data *d = (_strided_datetime_cast_data *)data;
+    npy_int64 dt;
+    npy_datetimestruct dts;
+
+    while (N > 0) {
+        memcpy(&dt, src, sizeof(dt));
+
+        if (convert_datetime_to_datetimestruct(&d->src_meta,
+                                               dt, &dts) < 0) {
+            dt = NPY_DATETIME_NAT;
+        }
+        else {
+            dts.event = dts.event % d->dst_meta.events;
+            if (convert_datetimestruct_to_datetime(&d->dst_meta,
+                                                   &dts, &dt) < 0) {
+                dt = NPY_DATETIME_NAT;
+            }
+        }
+
+        memcpy(dst, &dt, sizeof(dt));
+
+        dst += dst_stride;
+        src += src_stride;
+        --N;
+    }
+}
+
+static void
+_strided_to_strided_datetime_cast(char *dst, npy_intp dst_stride,
+                        char *src, npy_intp src_stride,
+                        npy_intp N, npy_intp src_itemsize,
+                        void *data)
+{
+    _strided_datetime_cast_data *d = (_strided_datetime_cast_data *)data;
+    npy_int64 num = d->num, denom = d->denom;
+    npy_int64 dt;
+    int event = 0, src_events = d->src_events, dst_events = d->dst_events;
+
+    while (N > 0) {
+        memcpy(&dt, src, sizeof(dt));
+
+        if (dt != NPY_DATETIME_NAT) {
+            /* Remove the event number from the value */
+            if (src_events > 1) {
+                event = (int)(dt % src_events);
+                dt = dt / src_events;
+                if (event < 0) {
+                    --dt;
+                    event += src_events;
+                }
+            }
+
+            /* Apply the scaling */
+            if (dt < 0) {
+                dt = (dt * num - (denom - 1)) / denom;
+            }
+            else {
+                dt = dt * num / denom;
+            }
+
+            /* Add the event number back in */
+            if (dst_events > 1) {
+                event = event % dst_events;
+                dt = dt * dst_events + event;
+            }
+        }
+
+        memcpy(dst, &dt, sizeof(dt));
+
+        dst += dst_stride;
+        src += src_stride;
+        --N;
+    }
+}
+
+static void
+_aligned_strided_to_strided_datetime_cast_no_events(char *dst,
+                        npy_intp dst_stride,
+                        char *src, npy_intp src_stride,
+                        npy_intp N, npy_intp src_itemsize,
+                        void *data)
+{
+    _strided_datetime_cast_data *d = (_strided_datetime_cast_data *)data;
+    npy_int64 num = d->num, denom = d->denom;
+    npy_int64 dt;
+
+    while (N > 0) {
+        dt = *(npy_int64 *)src;
+
+        if (dt != NPY_DATETIME_NAT) {
+            /* Apply the scaling */
+            if (dt < 0) {
+                dt = (dt * num - (denom - 1)) / denom;
+            }
+            else {
+                dt = dt * num / denom;
+            }
+        }
+
+        *(npy_int64 *)dst = dt;
+
+        dst += dst_stride;
+        src += src_stride;
+        --N;
+    }
+}
+
+/*
+ * Assumes src_dtype and dst_dtype are both datetimes or both timedeltas
+ */
+static int
+get_nbo_cast_datetime_transfer_function(int aligned,
+                            npy_intp src_stride, npy_intp dst_stride,
+                            PyArray_Descr *src_dtype, PyArray_Descr *dst_dtype,
+                            PyArray_StridedTransferFn **out_stransfer,
+                            void **out_transferdata)
+{
+    PyArray_DatetimeMetaData *src_meta, *dst_meta;
+    npy_int64 num = 0, denom = 0;
+    _strided_datetime_cast_data *data;
+
+    src_meta = get_datetime_metadata_from_dtype(src_dtype);
+    if (src_meta == NULL) {
+        return NPY_FAIL;
+    }
+    dst_meta = get_datetime_metadata_from_dtype(dst_dtype);
+    if (dst_meta == NULL) {
+        return NPY_FAIL;
+    }
+
+    get_datetime_conversion_factor(src_meta, dst_meta, &num, &denom);
+
+    if (num == 0) {
+        PyObject *errmsg;
+        errmsg = PyUString_FromString("Integer overflow "
+                    "getting a conversion factor between types ");
+        PyUString_ConcatAndDel(&errmsg,
+                PyObject_Repr((PyObject *)src_dtype));
+        PyUString_ConcatAndDel(&errmsg,
+                PyUString_FromString(" and "));
+        PyUString_ConcatAndDel(&errmsg,
+                PyObject_Repr((PyObject *)dst_dtype));
+        PyErr_SetObject(PyExc_OverflowError, errmsg);
+        return NPY_FAIL;
+    }
+
+    /* Allocate the data for the casting */
+    data = (_strided_datetime_cast_data *)PyArray_malloc(
+                                    sizeof(_strided_datetime_cast_data));
+    if (data == NULL) {
+        PyErr_NoMemory();
+        *out_stransfer = NULL;
+        *out_transferdata = NULL;
+        return NPY_FAIL;
+    }
+    data->freefunc = &PyArray_free;
+    data->copyfunc = &_strided_datetime_cast_data_copy;
+    data->num = num;
+    data->denom = denom;
+    data->src_events = src_meta->events;
+    data->dst_events = dst_meta->events;
+
+    /*
+     * Special case the datetime (but not timedelta) with the nonlinear
+     * units (years, months, business days). For timedelta, an average
+     * years and months value is used.
+     */
+    if (src_dtype->type_num == NPY_DATETIME &&
+            (src_meta->base == NPY_FR_Y ||
+             src_meta->base == NPY_FR_M ||
+             src_meta->base == NPY_FR_B ||
+             dst_meta->base == NPY_FR_Y ||
+             dst_meta->base == NPY_FR_M ||
+             dst_meta->base == NPY_FR_B)) {
+        memcpy(&data->src_meta, src_meta, sizeof(data->src_meta));
+        memcpy(&data->dst_meta, dst_meta, sizeof(data->dst_meta));
+        *out_stransfer = &_strided_to_strided_datetime_general_cast;
+    }
+    else if (aligned && data->src_events == 1 && data->dst_events == 1) {
+        *out_stransfer = &_aligned_strided_to_strided_datetime_cast_no_events;
+    }
+    else {
+        *out_stransfer = &_strided_to_strided_datetime_cast;
+    }
+    *out_transferdata = data;
+
+#if NPY_DT_DBG_TRACING
+    printf("Dtype transfer from ");
+    PyObject_Print((PyObject *)src_dtype, stdout, 0);
+    printf(" to ");
+    PyObject_Print((PyObject *)dst_dtype, stdout, 0);
+    printf("\n");
+    printf("has conversion fraction %lld/%lld\n", num, denom);
+#endif
+
+
+    return NPY_SUCCEED;
+}
+
 static int
 get_nbo_cast_transfer_function(int aligned,
                             npy_intp src_stride, npy_intp dst_stride,
@@ -722,6 +961,19 @@ get_nbo_cast_transfer_function(int aligned,
                                     out_stransfer, out_transferdata);
     }
 
+    /* As a parameterized type, datetime->datetime sometimes needs casting */
+    if ((src_dtype->type_num == NPY_DATETIME &&
+                dst_dtype->type_num == NPY_DATETIME) ||
+            (src_dtype->type_num == NPY_TIMEDELTA &&
+                dst_dtype->type_num == NPY_TIMEDELTA)) {
+        *out_needs_wrap = !PyArray_ISNBO(src_dtype->byteorder) ||
+                          !PyArray_ISNBO(dst_dtype->byteorder);
+        return get_nbo_cast_datetime_transfer_function(aligned,
+                                    src_stride, dst_stride,
+                                    src_dtype, dst_dtype,
+                                    out_stransfer, out_transferdata);
+    }
+
     *out_needs_wrap = !aligned ||
                       !PyArray_ISNBO(src_dtype->byteorder) ||
                       !PyArray_ISNBO(dst_dtype->byteorder);
@@ -854,7 +1106,9 @@ get_cast_transfer_function(int aligned,
     npy_intp src_itemsize = src_dtype->elsize,
             dst_itemsize = dst_dtype->elsize;
 
-    if (src_dtype->type_num == dst_dtype->type_num) {
+    if (src_dtype->type_num == dst_dtype->type_num &&
+            src_dtype->type_num != NPY_DATETIME &&
+            src_dtype->type_num != NPY_TIMEDELTA) {
         PyErr_SetString(PyExc_ValueError,
                 "low level cast function is for unequal type numbers");
         return NPY_FAIL;
diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index d885909d5..fe2718620 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -1786,7 +1786,8 @@ PyUFunc_BinaryComparisonTypeResolution(PyUFuncObject *ufunc,
     if (ufunc->nin != 2 || ufunc->nout != 1) {
         PyErr_Format(PyExc_RuntimeError, "ufunc %s is configured "
                 "to use binary comparison type resolution but has "
-                "the wrong number of inputs or outputs");
+                "the wrong number of inputs or outputs",
+                ufunc_name);
         return -1;
     }
 
@@ -4389,11 +4390,11 @@ ufunc_generic_call(PyUFuncObject *self, PyObject *args, PyObject *kwds)
         for (i = 0; i < self->nargs; i++) {
             PyArray_XDECREF_ERR(mps[i]);
         }
-        if (errval == -1)
+        if (errval == -1) {
             return NULL;
+        }
         else if (self->nin == 2 && self->nout == 1) {
-          /* To allow the other argument to be given a chance
-           */
+            /* To allow the other argument to be given a chance */
             Py_INCREF(Py_NotImplemented);
             return Py_NotImplemented;
         }
diff --git a/numpy/core/tests/test_datetime.py b/numpy/core/tests/test_datetime.py
index 67e0fa1b4..fde277d4f 100644
--- a/numpy/core/tests/test_datetime.py
+++ b/numpy/core/tests/test_datetime.py
@@ -26,6 +26,68 @@ class TestDateTime(TestCase):
         assert_raises(TypeError, np.dtype, 'M16')
         assert_raises(TypeError, np.dtype, 'm16')
 
+    def test_days_creation(self):
+        assert_equal(np.array('1599', dtype='M8[D]').astype('i8'),
+                (1600-1970)*365 - (1972-1600)/4 + 3 - 365)
+        assert_equal(np.array('1600', dtype='M8[D]').astype('i8'),
+                (1600-1970)*365 - (1972-1600)/4 + 3)
+        assert_equal(np.array('1601', dtype='M8[D]').astype('i8'),
+                (1600-1970)*365 - (1972-1600)/4 + 3 + 366)
+        assert_equal(np.array('1900', dtype='M8[D]').astype('i8'),
+                (1900-1970)*365 - (1970-1900)/4)
+        assert_equal(np.array('1901', dtype='M8[D]').astype('i8'),
+                (1900-1970)*365 - (1970-1900)/4 + 365)
+        assert_equal(np.array('1967', dtype='M8[D]').astype('i8'), -3*365 - 1)
+        assert_equal(np.array('1968', dtype='M8[D]').astype('i8'), -2*365 - 1)
+        assert_equal(np.array('1969', dtype='M8[D]').astype('i8'), -1*365)
+        assert_equal(np.array('1970', dtype='M8[D]').astype('i8'), 0*365)
+        assert_equal(np.array('1971', dtype='M8[D]').astype('i8'), 1*365)
+        assert_equal(np.array('1972', dtype='M8[D]').astype('i8'), 2*365)
+        assert_equal(np.array('1973', dtype='M8[D]').astype('i8'), 3*365 + 1)
+        assert_equal(np.array('1974', dtype='M8[D]').astype('i8'), 4*365 + 1)
+        assert_equal(np.array('2000', dtype='M8[D]').astype('i8'),
+                 (2000 - 1970)*365 + (2000 - 1972)/4)
+        assert_equal(np.array('2001', dtype='M8[D]').astype('i8'),
+                 (2000 - 1970)*365 + (2000 - 1972)/4 + 366)
+        assert_equal(np.array('2400', dtype='M8[D]').astype('i8'),
+                 (2400 - 1970)*365 + (2400 - 1972)/4 - 3)
+        assert_equal(np.array('2401', dtype='M8[D]').astype('i8'),
+                 (2400 - 1970)*365 + (2400 - 1972)/4 - 3 + 366)
+
+        assert_equal(np.array('1600-02-29', dtype='M8[D]').astype('i8'),
+                (1600-1970)*365 - (1972-1600)/4 + 3 + 31 + 28)
+        assert_equal(np.array('1600-03-01', dtype='M8[D]').astype('i8'),
+                (1600-1970)*365 - (1972-1600)/4 + 3 + 31 + 29)
+        assert_equal(np.array('2000-02-29', dtype='M8[D]').astype('i8'),
+                 (2000 - 1970)*365 + (2000 - 1972)/4 + 31 + 28)
+        assert_equal(np.array('2000-03-01', dtype='M8[D]').astype('i8'),
+                 (2000 - 1970)*365 + (2000 - 1972)/4 + 31 + 29)
+        assert_equal(np.array('2001-03-22', dtype='M8[D]').astype('i8'),
+                 (2000 - 1970)*365 + (2000 - 1972)/4 + 366 + 31 + 28 + 21)
+
+    def test_days_to_pydatetime(self):
+        assert_equal(np.array('1599', dtype='M8[D]').astype('O'),
+                    datetime.date(1599, 1, 1))
+        assert_equal(np.array('1600', dtype='M8[D]').astype('O'),
+                    datetime.date(1600, 1, 1))
+        assert_equal(np.array('1601', dtype='M8[D]').astype('O'),
+                    datetime.date(1601, 1, 1))
+        assert_equal(np.array('1900', dtype='M8[D]').astype('O'),
+                    datetime.date(1900, 1, 1))
+        assert_equal(np.array('1901', dtype='M8[D]').astype('O'),
+                    datetime.date(1901, 1, 1))
+        assert_equal(np.array('2000', dtype='M8[D]').astype('O'),
+                    datetime.date(2000, 1, 1))
+        assert_equal(np.array('2001', dtype='M8[D]').astype('O'),
+                    datetime.date(2001, 1, 1))
+        assert_equal(np.array('1600-02-29', dtype='M8[D]').astype('O'),
+                    datetime.date(1600, 2, 29))
+        assert_equal(np.array('1600-03-01', dtype='M8[D]').astype('O'),
+                    datetime.date(1600, 3, 1))
+        assert_equal(np.array('2001-03-22', dtype='M8[D]').astype('O'),
+                    datetime.date(2001, 3, 22))
+                
+
     def test_dtype_comparison(self):
         assert_(not (np.dtype('M8[us]') == np.dtype('M8[ms]')))
         assert_(np.dtype('M8[us]') != np.dtype('M8[ms]'))
@@ -134,11 +196,11 @@ class TestDateTime(TestCase):
                      np.array('1980-02-29T23:59:59.999999Z', dtype='M8[M]'))
 
     def test_different_unit_comparison(self):
-        # Check some years
+        # Check some years with units that won't overflow
         for unit1 in ['Y', 'M', 'D', '6h', 'h', 'm', 's', '10ms',
-                                                        'ms', 'us', 'ps']:
+                                                        'ms', 'us']:
             dt1 = np.dtype('M8[%s]' % unit1)
-            for unit2 in ['Y', 'M', 'D', 'h', 'm', 's', 'ms', 'us', 'ps']:
+            for unit2 in ['Y', 'M', 'D', 'h', 'm', 's', 'ms', 'us']:
                 dt2 = np.dtype('M8[%s]' % unit2)
                 assert_equal(np.array('1945', dtype=dt1),
                              np.array('1945', dtype=dt2))
@@ -148,15 +210,15 @@ class TestDateTime(TestCase):
                              np.array('9999', dtype=dt2))
                 assert_equal(np.array('10000', dtype=dt1),
                              np.array('10000-01-01', dtype=dt2))
-        # Check some days
-        for unit1 in ['D', '12h', 'h', 'm', 's', '4s', 'ms', 'us', 'ps']:
+        # Check some days with units that won't overflow
+        for unit1 in ['D', '12h', 'h', 'm', 's', '4s', 'ms', 'us']:
             dt1 = np.dtype('M8[%s]' % unit1)
-            for unit2 in ['D', 'h', 'm', 's', 'ms', 'us', 'ps']:
+            for unit2 in ['D', 'h', 'm', 's', 'ms', 'us']:
                 dt2 = np.dtype('M8[%s]' % unit2)
                 assert_equal(np.array('1932-02-17', dtype=dt1),
-                             np.array('1932-02-17T00:00:00', dtype=dt2))
+                             np.array('1932-02-17T00:00:00Z', dtype=dt2))
                 assert_equal(np.array('10000-04-27', dtype=dt1),
-                             np.array('10000-04-27T00:00:00', dtype=dt2))
+                             np.array('10000-04-27T00:00:00Z', dtype=dt2))
                 assert_equal(np.array('today', dtype=dt1),
                              np.array('today', dtype=dt2))