"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020b/lib/src/calendar.c" (18 Mar 2020, 28783 Bytes) of package /linux/misc/gretl-2020b.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "calendar.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2020a_vs_2020b.

    1 /*
    2  *  gretl -- Gnu Regression, Econometrics and Time-series Library
    3  *  Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
    4  *
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  *
   10  *  This program is distributed in the hope that it will be useful,
   11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13  *  GNU General Public License for more details.
   14  *
   15  *  You should have received a copy of the GNU General Public License
   16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   17  *
   18  */
   19 
   20 #include "libgretl.h"
   21 
   22 /**
   23  * SECTION:calendar
   24  * @short_description: functions for working with dates
   25  * @title: Calendar
   26  * @include: libgretl.h
   27  *
   28  * Here we have various functions dealing with calendar dates;
   29  * for the most part these are designed for handling daily
   30  * time-series data.
   31  */
   32 
   33 static int days_in_month[2][13] = {
   34     {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
   35     {0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
   36 };
   37 
   38 /* Jan 1, 0001, was a Monday on the proleptic Gregorian calendar */
   39 #define DAY1 1
   40 
   41 /* Note on the GLib API used below: where "julian" occurs in the names
   42    of GLib calendrical functions it refers to the "Julian day"; that is,
   43    the number of days since some fixed starting point, as used by
   44    astronomers. This is quite distinct from the Julian calendar.
   45    However, GLib takes Julian day 1 to be the first of January in
   46    AD 1, as opposed to the astronomical starting point in 4713 BC,
   47    so these are not strictly Julian days, and in our own functions
   48    which use the same concept we call them "epoch days".
   49 */
   50 
   51 static int leap_year (int yr)
   52 {
   53     return (!(yr % 4) && (yr % 100)) || !(yr % 400);
   54 }
   55 
   56 static int valid_ymd (int y, int m, int d, int julian)
   57 {
   58     int ok = 1;
   59 
   60     if (!g_date_valid_dmy(d, m, y)) {
   61     ok = 0;
   62     if (julian && y > 0 && m == 2 && d == 29 && y%4 == 0) {
   63         ok = 1;
   64     }
   65     }
   66 
   67     return ok;
   68 }
   69 
   70 static int day_of_week_from_ymd (int y, int m, int d, int julian)
   71 {
   72     GDate date;
   73     int wd;
   74 
   75     if (!valid_ymd(y, m, d, julian)) {
   76     return -1;
   77     }
   78 
   79     g_date_clear(&date, 1);
   80 
   81     if (julian) {
   82     guint32 ed = epoch_day_from_julian_ymd(y, m, d);
   83 
   84     g_date_set_julian(&date, ed);
   85     } else {
   86     g_date_set_dmy(&date, d, m, y);
   87     }
   88 
   89     wd = g_date_get_weekday(&date);
   90 
   91     /* switch to Sunday == 0 */
   92     return wd == G_DATE_SUNDAY ? 0 : wd;
   93 }
   94 
   95 /**
   96  * epoch_day_from_ymd:
   97  * @y: year (1 <= y <= 9999).
   98  * @m: month (1 <= m <= 12).
   99  * @d: day of month (1 <= d <= 31).
  100  *
  101  * Returns: the epoch day number, which equals 1 for the first of
  102  * January in the year AD 1 on the proleptic Gregorian calendar,
  103  * or 0 on error.
  104  */
  105 
  106 guint32 epoch_day_from_ymd (int y, int m, int d)
  107 {
  108     GDate date;
  109 
  110     if (!g_date_valid_dmy(d, m, y)) {
  111     return 0;
  112     }
  113 
  114     g_date_clear(&date, 1);
  115     g_date_set_dmy(&date, d, m, y);
  116 
  117     return g_date_get_julian(&date);
  118 }
  119 
  120 /**
  121  * ymd_bits_from_epoch_day:
  122  * @ed: epoch day (ed >= 1).
  123  * @y: location to receive year.
  124  * @m: location to receive month.
  125  * @m: location to receive day.
  126  *
  127  * Returns: 0 on success, non-zero on error.
  128  */
  129 
  130 int ymd_bits_from_epoch_day (guint32 ed, int *y, int *m, int *d)
  131 {
  132     GDate date;
  133 
  134     if (!g_date_valid_julian(ed)) {
  135     return E_INVARG;
  136     }
  137 
  138     g_date_clear(&date, 1);
  139     g_date_set_julian(&date, ed);
  140 
  141     *y = g_date_get_year(&date);
  142     *m = g_date_get_month(&date);
  143     *d = g_date_get_day(&date);
  144 
  145     return 0;
  146 }
  147 
  148 /**
  149  * julian_ymd_bits_from_epoch_day:
  150  * @ed: epoch day (ed >= 1).
  151  * @y: location to receive year.
  152  * @m: location to receive month.
  153  * @m: location to receive day.
  154  *
  155  * Follows the algorithm of E.G. Richards (2013), "Calendars," In S.E.
  156  * Urban & P.K. Seidelmann, eds. Explanatory Supplement to the Astronomical
  157  * Almanac, 3rd ed. (pp. 585-624), Mill Valley, CA: University Science Books
  158  * (as set out on https://en.wikipedia.org/wiki/Julian_day).
  159  *
  160  * There are other algorithms for this purpose on the internet but they are
  161  * mostly wrong (at least, not right for all dates); many of them fail
  162  * the round-trip test (date -> epoch day -> date) for some dates.
  163  *
  164  * Returns: 0 on success, non-zero on error.
  165  */
  166 
  167 int julian_ymd_bits_from_epoch_day (guint32 ed, int *py,
  168                     int *pm, int *pd)
  169 {
  170     int y = 4716;
  171     int p = 1461;
  172     int f = ed + 1721425 + 1401;
  173     int e = 4 * f + 3;
  174     int g = (e % p)/4;
  175     int h = 5 * g + 2;
  176 
  177     /* The addition of 1721425 above translates from our
  178        "epoch day" to Julian Day Number; the addition of 1401
  179        to the JDN is specified by Richards' algorithm.
  180     */
  181 
  182     *pd = (h % 153)/5 + 1;
  183     *pm = (h/153 + 2) % 12 + 1;
  184     *py = e/p - y + (14 - *pm)/12;
  185 
  186     return 0;
  187 }
  188 
  189 /**
  190  * epoch_day_from_julian_ymd:
  191  * @y: year (y >= 1).
  192  * @m: month (1 <= m <= 12).
  193  * @d: day of month (1 <= d <= 31).
  194  *
  195  * The @y, @m and @d arguments are assumed to refer to a date on
  196  * the Julian calendar. The conversion algorithm is taken from
  197  * https://en.wikipedia.org/wiki/Julian_day, where it appears to
  198  * be credited to the Department of Computer Science at UT, Austin.
  199  *
  200  * Returns: the epoch day number, which equals 1 for the first of
  201  * January in the year AD 1 on the proleptic Gregorian calendar,
  202  * or 0 on error.
  203  */
  204 
  205 guint32 epoch_day_from_julian_ymd (int y, int m, int d)
  206 {
  207     int a = (14 - m)/12;
  208     int jd;
  209 
  210     y = y + 4800 - a;
  211     m = m + 12*a - 3;
  212 
  213     jd = d + (153*m + 2)/5 + 365*y + y/4 - 32083;
  214 
  215     if (jd <= 1721425) {
  216     /* prior to AD 1 */
  217     return 0;
  218     } else {
  219     return (guint32) jd - 1721425;
  220     }
  221 }
  222 
  223 /**
  224  * ymd_extended_from_epoch_day:
  225  * @ed: epoch day (ed >= 1).
  226  * @julian: non-zero to use Julian calendar, otherwise Gregorian.
  227  * @err: location to receive error code.
  228  *
  229  * Returns: a string on the pattern YYYY-MM-DD (ISO 8601 extended
  230  * date format) given the epoch day number, which equals 1 for the
  231  * first of January in the year 1 AD, or NULL on error.
  232  */
  233 
  234 char *ymd_extended_from_epoch_day (guint32 ed, int julian, int *err)
  235 {
  236     char *ret = NULL;
  237     int y, m, d;
  238     int myerr;
  239 
  240     if (julian) {
  241     myerr = julian_ymd_bits_from_epoch_day(ed, &y, &m, &d);
  242     } else {
  243     myerr = ymd_bits_from_epoch_day(ed, &y, &m, &d);
  244     }
  245 
  246     if (!myerr) {
  247     ret = calloc(12, 1);
  248     if (ret == NULL) {
  249         myerr = E_ALLOC;
  250     } else {
  251         sprintf(ret, "%04d-%02d-%02d", y, m, d);
  252     }
  253     }
  254 
  255     if (err != NULL) {
  256     *err = myerr;
  257     }
  258 
  259     return ret;
  260 }
  261 
  262 /**
  263  * ymd_basic_from_epoch_day:
  264  * @ed: epoch day (ed >= 1).
  265  * @julian: non-zero to use Julian calendar, otherwise Gregorian.
  266  * @err: location to receive error code.
  267  *
  268  * Returns: an 8-digit number on the pattern YYYYMMDD (ISO 8601 basic
  269  * date format) given the epoch day number, which equals 1 for the
  270  * first of January in the year 1 AD, or #NADBL on error.
  271  */
  272 
  273 double ymd_basic_from_epoch_day (guint32 ed, int julian, int *err)
  274 {
  275     int y = 0, m = 0, d = 0;
  276 
  277     if (julian) {
  278     *err = julian_ymd_bits_from_epoch_day(ed, &y, &m, &d);
  279     } else {
  280     *err = ymd_bits_from_epoch_day(ed, &y, &m, &d);
  281     }
  282 
  283     if (*err) {
  284     return NADBL;
  285     } else {
  286     return 10000*y + 100*m + d;
  287     }
  288 }
  289 
  290 /**
  291  * epoch_day_from_ymd_basic:
  292  * @ymd: number that is supposed to represent YYYYMMDD.
  293  *
  294  * Returns: the epoch day number corresponding to @ymd, interpreted
  295  * as YYYYMMDD (ISO 8601 "basic") if possible, or 0 on error.
  296  */
  297 
  298 guint32 epoch_day_from_ymd_basic (double ymd)
  299 {
  300     gchar tmp[9];
  301     int y, m, d, n;
  302 
  303     n = g_snprintf(tmp, 9, "%.0f", ymd);
  304     if (n != 8) {
  305     return 0;
  306     }
  307 
  308     y = floor(ymd / 10000);
  309     ymd -= y * 10000;
  310     m = floor(ymd / 100);
  311     d = ymd - m * 100;
  312 
  313     return epoch_day_from_ymd(y, m, d);
  314 }
  315 
  316 /**
  317  * weekday_from_epoch_day:
  318  * @ed: epoch day (ed >= 1).
  319  *
  320  * Returns: the weekday (Sunday = 0) corrsponding to @ed,
  321  * or -1 on error.
  322  */
  323 
  324 int weekday_from_epoch_day (guint32 ed)
  325 {
  326     GDate date;
  327     int wd;
  328 
  329     if (!g_date_valid_julian(ed)) {
  330     return -1;
  331     }
  332 
  333     g_date_clear(&date, 1);
  334     g_date_set_julian(&date, ed);
  335     wd = g_date_get_weekday(&date);
  336 
  337     return wd == G_DATE_SUNDAY ? 0 : wd;
  338 }
  339 
  340 /**
  341  * get_epoch_day:
  342  * @datestr: string representation of calendar date, in form
  343  * YY[YY]-MM-DD.
  344  *
  345  * Returns: the epoch day number, or -1 on failure.
  346  */
  347 
  348 guint32 get_epoch_day (const char *datestr)
  349 {
  350     GDate date;
  351     int y, m, d, nf = 0;
  352     int ydigits = 0;
  353 
  354     if (strchr(datestr, '-')) {
  355     ydigits = strcspn(datestr, "-");
  356     nf = sscanf(datestr, YMD_READ_FMT, &y, &m, &d);
  357     } else if (strchr(datestr, '/')) {
  358     ydigits = strcspn(datestr, "/");
  359     nf = sscanf(datestr, "%d/%d/%d", &y, &m, &d);
  360     } else if (strlen(datestr) == 8) {
  361     ydigits = 4;
  362     nf = sscanf(datestr, "%4d%2d%2d", &y, &m, &d);
  363     }
  364 
  365     if (nf != 3 || (ydigits != 4 && ydigits != 2)) {
  366     return 0;
  367     }
  368 
  369     if (ydigits == 2) {
  370     y = FOUR_DIGIT_YEAR(y);
  371     }
  372 
  373     if (!g_date_valid_dmy(d, m, y)) {
  374     return 0;
  375     }
  376 
  377     g_date_clear(&date, 1);
  378     g_date_set_dmy(&date, d, m, y);
  379 
  380     return g_date_get_julian(&date);
  381 }
  382 
  383 /**
  384  * calendar_obs_number:
  385  * @datestr: string representation of calendar date, in form
  386  * YY[YY]/MM/DD.
  387  * @dset: pointer to dataset information.
  388  *
  389  * Returns: The zero-based observation number for the given
  390  * date within the current data set.
  391  */
  392 
  393 int calendar_obs_number (const char *datestr, const DATASET *dset)
  394 {
  395     guint32 ed0 = (guint32) dset->sd0;
  396     guint32 t = get_epoch_day(datestr);
  397 
  398 #ifdef CAL_DEBUG
  399     fprintf(stderr, "calendar_obs_number: '%s' gave epoch day = %u\n",
  400         datestr, t);
  401 #endif
  402 
  403     if (t <= 0 || t < ed0) {
  404     return -1;
  405     } else if (t == ed0) {
  406     return 0;
  407     }
  408 
  409     /* subtract starting day for dataset */
  410     t -= ed0;
  411 
  412     if (t <= 0) {
  413     return -1;
  414     }
  415 
  416     if (dset->pd == 52) {
  417     /* weekly data */
  418     t /= 7;
  419     } else if (dset->pd == 5 || dset->pd == 6) {
  420     /* daily, 5- or 6-day week: subtract number of irrelevant days */
  421     int startday = (ed0 - 1 + DAY1) % 7;
  422     int wkends = (t + startday - 1) / 7;
  423 
  424 #ifdef CAL_DEBUG
  425     printf("calendar_obs_number: ed0=%d, date=%s, t=%d, startday=%d, wkends=%d\n",
  426            (int) ed0, date, (int) t, startday, wkends);
  427 #endif
  428 
  429     if (dset->pd == 5) {
  430         t -= (2 * wkends);
  431     } else {
  432         t -= wkends;
  433     }
  434     }
  435 
  436     return (int) t;
  437 }
  438 
  439 /* convert from $t in dataset to epoch day */
  440 
  441 static int t_to_epoch_day (int t, guint32 start, int wkdays)
  442 {
  443     int startday = (start - 1 + DAY1) % 7;
  444     int wkends = (t + startday - 1) / wkdays;
  445 
  446     if (wkdays == 5) {
  447     wkends *= 2;
  448     }
  449 
  450     return start + t + wkends;
  451 }
  452 
  453 /**
  454  * epoch_day_from_t:
  455  * @t: 0-based observation index.
  456  * @dset: pointer to dataset.
  457  *
  458  * Returns: the epoch day based on calendrical information
  459  * in @dset. In case of error 0 is returned.
  460  */
  461 
  462 guint32 epoch_day_from_t (int t, const DATASET *dset)
  463 {
  464     guint32 d0 = (guint32) dset->sd0;
  465     guint32 dt = 0;
  466 
  467     if (dset->pd == 52) {
  468     dt = d0 + 7 * t;
  469     } else if (dset->pd == 7) {
  470     dt = d0 + t;
  471     } else {
  472     dt = t_to_epoch_day(t, d0, dset->pd);
  473     }
  474 
  475     return dt;
  476 }
  477 
  478 /**
  479  * calendar_date_string:
  480  * @targ: string to be filled out.
  481  * @t: zero-based index of observation.
  482  * @dset: pointer to dataset.
  483  *
  484  * Writes to @targ the calendar representation of the date of
  485  * observation @t, in the form YY[YY]/MM/DD.
  486  *
  487  * Returns: 0 on success, non-zero on error.
  488  */
  489 
  490 int calendar_date_string (char *targ, int t, const DATASET *dset)
  491 {
  492     guint32 d0, dt = 0;
  493     int y, m, d;
  494     int err = 0;
  495 
  496     *targ = '\0';
  497     d0 = (guint32) dset->sd0;
  498 
  499     if (dset->pd == 52) {
  500     dt = d0 + 7 * t;
  501     } else if (dset->pd == 7) {
  502     dt = d0 + t;
  503     } else {
  504     /* 5- or 6-day data */
  505     if (t == 0 && dset->pd == 5) {
  506         int wd = weekday_from_epoch_day(d0);
  507 
  508         if (wd == 0 || wd == 6) {
  509         gretl_errmsg_sprintf(_("Invalid starting date for %d-day data"), dset->pd);
  510         return E_DATA;
  511         }
  512     }
  513     dt = t_to_epoch_day(t, d0, dset->pd);
  514     }
  515 
  516     err = ymd_bits_from_epoch_day(dt, &y, &m, &d);
  517 
  518     if (!err) {
  519     if (strlen(dset->stobs) == 8) {
  520         sprintf(targ, YMD_WRITE_Y2_FMT, y % 100, m, d);
  521     } else {
  522         sprintf(targ, YMD_WRITE_Y4_FMT, y, m, d);
  523     }
  524     }
  525 
  526     return err;
  527 }
  528 
  529 /**
  530  * MS_excel_date_string:
  531  * @targ: date string to be filled out.
  532  * @mst: MS Excel-type date code: days since base.
  533  * @pd: periodicity of data (or 0 if unknown).
  534  * @d1904: set to 1 if the base is 1904/01/01; otherwise
  535  * the base is assumed to be 1899/12/31.
  536  *
  537  * Writes to @targ the calendar representation of the date of
  538  * observation @mst, in the form YYYY-MM-DD if @pd is 0, 5,
  539  * 6, 7 or 52 (unknown, daily, or weekly frequency), otherwise
  540  * in the appropriate format for annual, quarterly or monthly
  541  * according to @pd.
  542  *
  543  * Returns: 0.
  544  */
  545 
  546 int MS_excel_date_string (char *targ, int mst, int pd, int d1904)
  547 {
  548     int y = (d1904)? 1904 : 1900;
  549     int d = (d1904)? 2 : 1;
  550     int m = 1;
  551     int leap, drem;
  552 
  553     *targ = '\0';
  554 
  555     if (mst == 0) {
  556     /* date coincident with base */
  557     if (d1904) {
  558         d = 1;
  559     } else {
  560         y = 1899;
  561         m = 12;
  562         d = 31;
  563     }
  564     } else if (mst > 0) {
  565     /* date subsequent to base */
  566     drem = mst + d1904;
  567 
  568     while (1) {
  569         int yd = 365 + leap_year(y);
  570 
  571         /* MS nincompoopery */
  572         if (y == 1900) yd++;
  573 
  574         if (drem > yd) {
  575         drem -= yd;
  576         y++;
  577         } else {
  578         break;
  579         }
  580     }
  581 
  582     leap = leap_year(y) + (y == 1900);
  583 
  584     for (m=1; m<=12; m++) {
  585         int md = days_in_month[leap][m];
  586 
  587         if (drem > md) {
  588         drem -= md;
  589         } else {
  590         d = drem;
  591         break;
  592         }
  593     }
  594     } else {
  595     /* mst < 0, date prior to base */
  596     drem = -(mst + d1904);
  597 
  598     y = (d1904)? 1903 : 1899;
  599 
  600     while (1) {
  601         int yd = 365 + leap_year(y);
  602 
  603         if (drem > yd) {
  604         drem -= yd;
  605         y--;
  606         } else {
  607         break;
  608         }
  609     }
  610 
  611     leap = leap_year(y);
  612 
  613     for (m=12; m>0; m--) {
  614         int md = days_in_month[leap][m];
  615 
  616         if (drem >= md) {
  617         drem -= md;
  618         } else {
  619         d = md - drem;
  620         break;
  621         }
  622     }
  623     }
  624 
  625     if (pd == 1) {
  626     sprintf(targ, "%d", y);
  627     } else if (pd == 12) {
  628     sprintf(targ, "%d:%02d", y, m);
  629     } else if (pd == 4) {
  630     int q = 1 + m / 3.25;
  631 
  632     sprintf(targ, "%d:%d", y, q);
  633     } else {
  634     sprintf(targ, YMD_WRITE_Y4_FMT, y, m, d);
  635     }
  636 
  637     return 0;
  638 }
  639 
  640 /**
  641  * get_dec_date:
  642  * @datestr: calendar representation of date: YYYY-MM-DD.
  643  *
  644  * Returns: representation of date as year plus fraction of year.
  645  */
  646 
  647 double get_dec_date (const char *datestr)
  648 {
  649     GDate date;
  650     int y, m, d, nf;
  651     int ydigits = 0;
  652     double yrday, frac;
  653 
  654     nf = sscanf(datestr, YMD_READ_FMT, &y, &m, &d);
  655 
  656     if (nf == 3) {
  657     ydigits = strcspn(datestr, "-");
  658     } else if (strchr(datestr, '/') != NULL) {
  659     /* backward compatibility */
  660     ydigits = strcspn(datestr, "/");
  661     nf = sscanf(datestr, "%d/%d/%d", &y, &m, &d);
  662     }
  663 
  664     if (nf != 3 || (ydigits != 4 && ydigits != 2)) {
  665     return NADBL;
  666     }
  667 
  668     if (ydigits == 2) {
  669     y = FOUR_DIGIT_YEAR(y);
  670     }
  671 
  672     if (!g_date_valid_dmy(d, m, y)) {
  673     return NADBL;
  674     }
  675 
  676     g_date_clear(&date, 1);
  677     g_date_set_dmy(&date, d, m, y);
  678     yrday = g_date_get_day_of_year(&date);
  679     frac = yrday / (365 + g_date_is_leap_year(y));
  680 
  681     return y + frac;
  682 }
  683 
  684 /**
  685  * day_of_week:
  686  * @y: year.
  687  * @m: month, 1 to 12.
  688  * @d: day in month, 1 to 31.
  689  * @julian: non-zero to use Julian calendar, otherwise Gregorian.
  690  * @err: location to receive error code.
  691  *
  692  * Returns: the day of the week for the supplied date
  693  * (Sunday = 0, Monday = 1, ...) or %NADBL on failure
  694  * (the date is invalid).
  695  */
  696 
  697 double day_of_week (int y, int m, int d, int julian, int *err)
  698 {
  699     int wd = day_of_week_from_ymd(y, m, d, julian);
  700 
  701     return wd < 0 ? NADBL : wd;
  702 }
  703 
  704 #define day_in_calendar(w, d) (((w) == 6 && d != 0) || \
  705                    ((w) == 5 && d != 0 && d != 6))
  706 
  707 /**
  708  * weekday_from_date:
  709  * @datestr: calendar representation of date, [YY]YY/MM/DD
  710  *
  711  * Returns: day of week as integer, Sunday = 0.
  712  */
  713 
  714 int weekday_from_date (const char *datestr)
  715 {
  716     int y, m, d;
  717     int ydigits;
  718 
  719     if (sscanf(datestr, YMD_READ_FMT, &y, &m, &d) != 3) {
  720     return -1;
  721     }
  722 
  723     ydigits = strcspn(datestr, "-");
  724     if (ydigits != 4 && ydigits != 2) {
  725     return -1;
  726     }
  727 
  728     if (ydigits == 2) {
  729     y = FOUR_DIGIT_YEAR(y);
  730     }
  731 
  732     return day_of_week_from_ymd(y, m, d, 0);
  733 }
  734 
  735 /**
  736  * day_starts_month:
  737  * @d: day of month, 1-based
  738  * @m: month number, 1-based
  739  * @y: 4-digit year
  740  * @wkdays: number of days in week (7, 6 or 5)
  741  * @pad: location to receive 1 if the first day of the month
  742  * can reasonably be padded by a missing value (Jan 1), or NULL.
  743  *
  744  * Returns: 1 if the day is the "first day of the month",
  745  * allowance made for the possibility of a 5- or 6-day week,
  746  * else 0.
  747  */
  748 
  749 int day_starts_month (int d, int m, int y, int wkdays, int *pad)
  750 {
  751     int ret = 0;
  752 
  753     if (wkdays == 7) {
  754     if (d == 1) {
  755         ret = 1;
  756     } else if (pad != NULL && m == 1 && d == 2) {
  757         /* second of January */
  758         *pad = 1;
  759         ret = 1;
  760     }
  761     } else {
  762     /* 5- or 6-day week: check for first weekday or non-Sunday */
  763     int i, idx = day_of_week_from_ymd(y, m, 1, 0);
  764 
  765     for (i=1; i<6; i++) {
  766        if (day_in_calendar(wkdays, idx % 7)) {
  767            break;
  768        }
  769        idx++;
  770     }
  771     if (d == i) {
  772         ret = 1;
  773     } else if (pad != NULL && m == 1 && d == i + 1) {
  774         /* January again */
  775         *pad = 1;
  776         ret = 1;
  777     }
  778     }
  779 
  780     return ret;
  781 }
  782 
  783 /**
  784  * day_ends_month:
  785  * @d: day of month, 1-based
  786  * @m: month number, 1-based
  787  * @y: 4-digit year
  788  * @wkdays: number of days in week (7, 6 or 5)
  789  *
  790  * Returns: 1 if the day is the "last day of the month",
  791  * allowance made for the possibility of a 5- or 6-day week, else 0.
  792  */
  793 
  794 int day_ends_month (int d, int m, int y, int wkdays)
  795 {
  796     int ret = 0;
  797     int leap = (m == 2)? leap_year(y) : 0;
  798     int dm = days_in_month[leap][m];
  799 
  800     if (wkdays == 7) {
  801     ret = (d == dm);
  802     } else {
  803     /* 5- or 6-day week: check for last weekday or non-Sunday */
  804     int i, idx = day_of_week_from_ymd(y, m, dm, 0);
  805 
  806     for (i=dm; i>0; i--) {
  807         if (day_in_calendar(wkdays, idx % 7)) {
  808         break;
  809         }
  810         idx--;
  811     }
  812     ret = (d == i);
  813     }
  814 
  815     return ret;
  816 }
  817 
  818 /**
  819  * get_days_in_month:
  820  * @m: month number, 1-based
  821  * @y: 4-digit year
  822  * @wkdays: number of days in week (7, 6 or 5)
  823  * @julian: non-zero for Julian calendar, otherwise Gregorian.
  824  *
  825  * Returns: the number of (relevant) days in the month, allowance
  826  * made for the possibility of a 5- or 6-day week.
  827  */
  828 
  829 int get_days_in_month (int m, int y, int wkdays, int julian)
  830 {
  831     int dm, leap = 0;
  832     int ret = 0;
  833 
  834     if (m == 2) {
  835     leap = julian ? (y%4 == 0) : leap_year(y);
  836     }
  837     dm = days_in_month[leap][m];
  838 
  839     if (wkdays == 7) {
  840     ret = dm;
  841     } else {
  842     int i, idx = day_of_week_from_ymd(y, m, 1, julian);
  843 
  844     for (i=0; i<dm; i++) {
  845         if (day_in_calendar(wkdays, idx % 7)) {
  846         ret++;
  847         }
  848         idx++;
  849     }
  850     }
  851 
  852     return ret;
  853 }
  854 
  855 /**
  856  * month_day_index:
  857  * @y: 4-digit year
  858  * @m: month number, 1-based
  859  * @d: day in month.
  860  * @wkdays: number of days in week (7, 6 or 5)
  861  *
  862  * Returns: the 1-based index of calendar day @d in month
  863  * @m of year @y, allowing for the possibility of a 5- or
  864  * 6-day week.
  865  */
  866 
  867 int month_day_index (int y, int m, int d, int wkdays)
  868 {
  869     int ret = 0;
  870 
  871     if (wkdays == 7) {
  872     ret = d;
  873     } else {
  874     int i, idx = day_of_week_from_ymd(y, m, 1, 0);
  875 
  876     for (i=1; i<=d; i++) {
  877         if (day_in_calendar(wkdays, idx)) {
  878         ret++;
  879         }
  880         idx = (idx == 6)? 0 : idx + 1;
  881     }
  882     }
  883 
  884     return ret;
  885 }
  886 
  887 /**
  888  * days_in_month_before:
  889  * @y: 4-digit year
  890  * @m: month number, 1-based
  891  * @d: day in month.
  892  * @wkdays: number of days in week (7, 6 or 5)
  893  *
  894  * Returns: the number of relevant days in the month prior to
  895  * the supplied date, allowing for the possibility of a 5- or
  896  * 6-day week.
  897  */
  898 
  899 int days_in_month_before (int y, int m, int d, int wkdays)
  900 {
  901     int ret = 0;
  902 
  903     if (wkdays == 7) {
  904     ret = d - 1;
  905     } else {
  906     int i, idx = day_of_week_from_ymd(y, m, 1, 0);
  907 
  908     for (i=1; i<d; i++) {
  909         if (day_in_calendar(wkdays, idx % 7)) {
  910         ret++;
  911         }
  912         idx++;
  913     }
  914     }
  915 
  916     return ret;
  917 }
  918 
  919 /**
  920  * days_in_month_after:
  921  * @y: 4-digit year
  922  * @m: month number, 1-based
  923  * @d: day in month.
  924  * @wkdays: number of days in week (7, 6 or 5)
  925  *
  926  * Returns: the number of relevant days in the month after
  927  * the supplied date, allowing for the possibility of a 5- or
  928  * 6-day week.
  929  */
  930 
  931 int days_in_month_after (int y, int m, int d, int wkdays)
  932 {
  933     int leap = (m == 2)? leap_year(y) : 0;
  934     int dm = days_in_month[leap][m];
  935     int ret = 0;
  936 
  937     if (wkdays == 7) {
  938     ret = dm - d;
  939     } else {
  940     int i, wd = day_of_week_from_ymd(y, m, dm, 0);
  941 
  942     for (i=dm; i>d; i--) {
  943         if (day_in_calendar(wkdays, wd)) {
  944         ret++;
  945         }
  946         if (wd > 0) {
  947         wd--;
  948         } else {
  949         wd = 6;
  950         }
  951     }
  952     }
  953 
  954     return ret;
  955 }
  956 
  957 /**
  958  * date_to_daily_index:
  959  * @datestr: date in format YYYY-MM-DD.
  960  * @wkdays: number of days in week (7, 6 or 5)
  961  *
  962  * Returns: the zero-based index of the specified day
  963  * within the specified month and year. In the case
  964  * of 5- or 6-day data index zero does not necessarily
  965  * correspond to the first day of the month but rather
  966  * to the first relevant day.
  967  */
  968 
  969 int date_to_daily_index (const char *datestr, int wkdays)
  970 {
  971     int y, m, d, seq = 0;
  972 
  973     if (sscanf(datestr, YMD_READ_FMT, &y, &m, &d) != 3) {
  974     return -1;
  975     }
  976 
  977     if (wkdays == 7) {
  978     seq = d - 1;
  979     } else {
  980     int leap = leap_year(y);
  981     int n = days_in_month[leap][m];
  982     int i, idx = day_of_week_from_ymd(y, m, 1, 0);
  983 
  984     for (i=1; i<=n; i++) {
  985         if (d == i) {
  986         break;
  987         }
  988         if (day_in_calendar(wkdays, idx % 7)) {
  989         seq++;
  990         }
  991         idx++;
  992     }
  993     }
  994 
  995     return seq;
  996 }
  997 
  998 /**
  999  * daily_index_to_date:
 1000  * @targ: location to receive the date (YYYY-MM-DD).
 1001  * @y: year.
 1002  * @m: month.
 1003  * @idx: zero-based index of day within month.
 1004  * @wkdays: number of days in week (7, 6 or 5)
 1005  *
 1006  * Fills out @targ with the calendar data implied by
 1007  * the specification of @y, @m, @seq and @wkdays,
 1008  * provided this specification corresponds to an actual
 1009  * calendar date.
 1010 
 1011  * Returns: 0 on successful completion, non-zero if
 1012  * there is no such calendar date.
 1013  */
 1014 
 1015 int daily_index_to_date (char *targ, int y, int m, int idx,
 1016              int wkdays)
 1017 {
 1018     int day = 0;
 1019 
 1020     *targ = '\0';
 1021 
 1022     if (m < 1 || m > 12 || idx < 0 || idx > 30) {
 1023     fprintf(stderr, "daily_index_to_date: y=%d, m=%d, seq=%d\n",
 1024         y, m, idx);
 1025     return E_INVARG;
 1026     }
 1027 
 1028     if (wkdays == 7) {
 1029     day = idx + 1;
 1030     } else {
 1031     int leap = leap_year(y);
 1032     int n = days_in_month[leap][m];
 1033     int wd = day_of_week_from_ymd(y, m, 1, 0);
 1034     int i, seq = 0;
 1035 
 1036     for (i=1; i<=n; i++) {
 1037         if (day_in_calendar(wkdays, wd % 7)) {
 1038         if (seq == idx) {
 1039             day = i;
 1040             break;
 1041         }
 1042         seq++;
 1043         }
 1044         wd++;
 1045     }
 1046     }
 1047 
 1048     if (day <= 0) {
 1049     return E_DATA;
 1050     } else {
 1051     sprintf(targ, YMD_WRITE_FMT, y, m, day);
 1052     return 0;
 1053     }
 1054 }
 1055 
 1056 /**
 1057  * n_hidden_missing_obs:
 1058  * @dset: dataset information.
 1059  * @t1: first observation.
 1060  * @t2: last observation.
 1061  *
 1062  * For daily data with user-supplied data strings,
 1063  * determine the number of "hidden" missing observations
 1064  * in the range @t1 to @t2 inclusive. This is
 1065  * the difference between the actual number of
 1066  * observations and the number that should be there,
 1067  * according to the calendar. Allowance is made for
 1068  * 5- or 6-day data, via the data frequency given
 1069  * in @dset.
 1070  *
 1071  * Returns: number of hidden observations.
 1072  */
 1073 
 1074 int n_hidden_missing_obs (const DATASET *dset, int t1, int t2)
 1075 {
 1076     int n_present = t2 - t1 + 1;
 1077     int cal_n;
 1078 
 1079     if (!dated_daily_data(dset) || dset->S == NULL) {
 1080     return 0;
 1081     }
 1082 
 1083     t1 = calendar_obs_number(dset->S[t1], dset);
 1084     t2 = calendar_obs_number(dset->S[t2], dset);
 1085 
 1086     cal_n = t2 - t1 + 1;
 1087 
 1088     return cal_n - n_present;
 1089 }
 1090 
 1091 /**
 1092  * guess_daily_pd:
 1093  * @dset: dataset information.
 1094  *
 1095  * Based on user-supplied daily date strings recorded in
 1096  * @dset, try to guess whether the number of observations
 1097  * per week is 5, 6 or 7 (given that some observations
 1098  * may be missing).
 1099  *
 1100  * Returns: best quess at data frequency.
 1101  */
 1102 
 1103 int guess_daily_pd (const DATASET *dset)
 1104 {
 1105     int t, wd, pd = 5;
 1106     int wdbak = -1;
 1107     int havesat = 0;
 1108     int gotsat = 0, gotsun = 0;
 1109     int contig = 0;
 1110 
 1111     wd = weekday_from_date(dset->S[0]);
 1112     if (6 - wd < dset->n) {
 1113     havesat = 1;
 1114     }
 1115 
 1116     for (t=0; t<dset->n && t<28; t++) {
 1117     wd = weekday_from_date(dset->S[t]);
 1118     if (wd == 0) {
 1119         gotsun = 1;
 1120     } else if (wd == 6) {
 1121         gotsat = 1;
 1122     }
 1123     if ((wdbak + 1) % 7 == wd) {
 1124         contig++;
 1125     }
 1126     wdbak = wd;
 1127     }
 1128 
 1129     if (gotsat && gotsun) {
 1130     pd = 7;
 1131     } else if (contig > 10) {
 1132     if (gotsun) pd = 7;
 1133     else if (gotsat) pd = 6;
 1134     } else if (dset->n > 7) {
 1135     if (!gotsun && !gotsat) {
 1136         pd = 5;
 1137     } else if (!gotsun) {
 1138         pd = 6;
 1139     }
 1140     } else if (havesat && !gotsat) {
 1141     pd = 5;
 1142     } else {
 1143     pd = 7;
 1144     }
 1145 
 1146     return pd;
 1147 }
 1148 
 1149 /**
 1150  * iso_basic_to_extended:
 1151  * @b: source array of YYYYMMDD values.
 1152  * @y: array to hold year values.
 1153  * @m: array to hold month values.
 1154  * @d: array to hold day-of-week values.
 1155  * @n: length of all the above arrays.
 1156  *
 1157  * Given the array @b of ISO 8601 "basic" daily dates (YYYYMMDD as
 1158  * doubles), fill out the arrays @y, @m and @d with year, month
 1159  * and day.
 1160  *
 1161  * Returns: 0.
 1162  */
 1163 
 1164 int iso_basic_to_extended (const double *b, double *y, double *m,
 1165                double *d, int n)
 1166 {
 1167     int bi, yi, mi, di;
 1168     int i, julian;
 1169 
 1170     for (i=0; i<n; i++) {
 1171     julian = 0;
 1172     if (na(b[i])) {
 1173         y[i] = m[i] = NADBL;
 1174         if (d != NULL) {
 1175         d[i] = NADBL;
 1176         }
 1177     } else {
 1178         bi = (int) b[i];
 1179         if (bi < 0) {
 1180         julian = 1;
 1181         bi = -bi;
 1182         }
 1183         yi = bi / 10000;
 1184         mi = (bi - 10000*yi) / 100;
 1185         di = bi - 10000*yi - 100*mi;
 1186         /* now check for legit date */
 1187         if (!valid_ymd(yi, mi, di, julian)) {
 1188         y[i] = m[i] = NADBL;
 1189         if (d != NULL) {
 1190             d[i] = NADBL;
 1191         }
 1192         } else {
 1193         y[i] = yi;
 1194         m[i] = mi;
 1195         if (d != NULL) {
 1196             d[i] = di;
 1197         }
 1198         }
 1199     }
 1200     }
 1201 
 1202     return 0;
 1203 }
 1204 
 1205 /**
 1206  * easterdate:
 1207  * @year: year for which we want Easter date (Gregorian).
 1208  *
 1209  * Algorithm taken from Wikipedia page
 1210  * https://en.wikipedia.org/wiki/Computus
 1211  * under the heading "Anonymous Gregorian algorithm".
 1212  *
 1213  * Returns: the date of Easter in the Gregorian calendar as
 1214  * (month + day/100). Note that April the 10th is, under
 1215  * this convention, 4.1; hence, 4.2 is April the 20th, not
 1216  * April the 2nd (which would be 4.02).
 1217  */
 1218 
 1219 double easterdate (int year)
 1220 {
 1221     int a = year % 19;
 1222     int b = year / 100;
 1223     int c = year % 100;
 1224     int d = b / 4;
 1225     int e = b % 4;
 1226     int f = (b + 8) / 25;
 1227     int g = (b - f + 1) / 3;
 1228     int h = (19 * a + b - d - g + 15) % 30;
 1229     int i = c / 4;
 1230     int k = c % 4;
 1231     int L = (32 + 2 * e + 2 * i - h - k) % 7;
 1232     int m = (a + 11 * h + 22 * L) / 451 ;
 1233 
 1234     int month = (h + L - 7 * m + 114) / 31;
 1235     int day = ((h + L - 7 * m + 114) % 31) + 1;
 1236 
 1237     return month + day * 0.01;
 1238 }
 1239 
 1240 /**
 1241  * dayspan:
 1242  * @ed1: first epoch day.
 1243  * @ed2: last epoch day.
 1244  * @wkdays: relevant days per week (5, 6 or 7).
 1245  * @err: location to receive error code.
 1246  *
 1247  * Returns: The number of days in the interval @ed1 to
 1248  * @ed2, inclusive, taking account of the number of daily
 1249  * observations per week, @wkdays. If @wkdays = 6 Sundays
 1250  * are disregarded; if @wkdays = 5 both Saturdays and
 1251  * Sundays are disregarded.
 1252  */
 1253 
 1254 int day_span (guint32 ed1, guint32 ed2, int wkdays, int *err)
 1255 {
 1256     int n = 0;
 1257 
 1258     if (!g_date_valid_julian(ed1) ||
 1259     !g_date_valid_julian(ed2) ||
 1260     ed2 < ed1) {
 1261     *err = E_INVARG;
 1262     } else if (wkdays == 7) {
 1263     /* simple! */
 1264     n = ed2 - ed1 + 1;
 1265     } else {
 1266     GDate date;
 1267     guint32 i;
 1268     int idx;
 1269 
 1270     g_date_clear(&date, 1);
 1271     g_date_set_julian(&date, ed1);
 1272     idx = g_date_get_weekday(&date) - 1;
 1273 
 1274     for (i=ed1; i<=ed2; i++) {
 1275         if (day_in_calendar(wkdays, idx % 7)) {
 1276         n++;
 1277         }
 1278         idx++;
 1279     }
 1280     }
 1281 
 1282     return n;
 1283 }
 1284 
 1285 /**
 1286  * iso_week_number:
 1287  * @y: year.
 1288  * @m: month (1 to 12).
 1289  * @d: day on month.
 1290  * @err: location to receive error code.
 1291  *
 1292  * Returns: The ISO 8601 week number (1 to 53) corresponding
 1293  * to the Gregorian date specified by the first 3 arguments,
 1294  * or -1 on invalid input.
 1295  */
 1296 
 1297 int iso_week_number (int y, int m, int d, int *err)
 1298 {
 1299     GDate date;
 1300     int wnum = -1;
 1301 
 1302     if (!g_date_valid_dmy(d, m, y)) {
 1303     *err = E_INVARG;
 1304     } else {
 1305     g_date_clear(&date, 1);
 1306     g_date_set_dmy(&date, d, m, y);
 1307     wnum = g_date_get_iso8601_week_of_year(&date);
 1308     }
 1309 
 1310     return wnum;
 1311 }
 1312 
 1313 /**
 1314  * iso_week_from_date:
 1315  * @datestr: date in ISO 8601 format, YYYY-MM-DD.
 1316  *
 1317  * Returns: The ISO 8601 week number (1 to 53) corresponding
 1318  * to the Gregorian date specified by @datestr, or -1 on
 1319  * invalid input.
 1320  */
 1321 
 1322 int iso_week_from_date (const char *datestr)
 1323 {
 1324     int y, m, d;
 1325     int ydigits;
 1326     int err = 0;
 1327 
 1328     if (sscanf(datestr, YMD_READ_FMT, &y, &m, &d) != 3) {
 1329     return -1;
 1330     }
 1331 
 1332     ydigits = strcspn(datestr, "-");
 1333 
 1334     if (ydigits != 4 && ydigits != 2) {
 1335     return -1;
 1336     } else if (ydigits == 2) {
 1337     y = FOUR_DIGIT_YEAR(y);
 1338     }
 1339 
 1340     return iso_week_number(y, m, d, &err);
 1341 }