src/language/expressions/helpers.c

   1 #include <config.h>
   2 #include "helpers.h"
   3 #include <gsl/gsl_roots.h>
   4 #include <gsl/gsl_sf.h>
   5 #include <libpspp/pool.h>
   6 #include "private.h"
   7
   8 const struct substring empty_string = {NULL, 0};
   9
  10 static void
  11 expr_error (void *aux UNUSED, const char *format, ...)
  12 {
  13   struct msg m;
  14   va_list args;
  15
  16   /* FIXME: we can do better about saying where the error
  17      occurred. */
  18   m.category = MSG_SYNTAX;
  19   m.severity = MSG_ERROR;
  20   msg_location (&m.where);
  21   va_start (args, format);
  22   m.text = xvasprintf (format, args);
  23   va_end (args);
  24
  25   msg_emit (&m);
  26 }
  27
  28 double
  29 expr_ymd_to_ofs (double year, double month, double day)
  30 {
  31   int y = year;
  32   int m = month;
  33   int d = day;
  34
  35   if (y != year || m != month || d != day)
  36     {
  37       msg (SE, _("One of the arguments to a DATE function is not an integer.  "
  38                  "The result will be system-missing."));
  39       return SYSMIS;
  40     }
  41
  42   return calendar_gregorian_to_offset (y, m, d, expr_error, NULL);
  43 }
  44
  45 double
  46 expr_ymd_to_date (double year, double month, double day)
  47 {
  48   double ofs = expr_ymd_to_ofs (year, month, day);
  49   return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
  50 }
  51
  52 double
  53 expr_wkyr_to_date (double week, double year)
  54 {
  55   int w = week;
  56
  57   if (w != week)
  58     {
  59       msg (SE, _("The week argument to DATE.WKYR is not an integer.  "
  60                  "The result will be system-missing."));
  61       return SYSMIS;
  62     }
  63   else if (w < 1 || w > 53)
  64     {
  65       msg (SE, _("The week argument to DATE.WKYR is outside the acceptable "
  66                  "range of 1 to 53.  "
  67                  "The result will be system-missing."));
  68       return SYSMIS;
  69     }
  70   else
  71     {
  72       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
  73       if (yr_1_1 != SYSMIS)
  74         return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
  75       else
  76         return SYSMIS;
  77     }
  78 }
  79
  80 double
  81 expr_yrday_to_date (double year, double yday)
  82 {
  83   int yd = yday;
  84
  85   if (yd != yday)
  86     {
  87       msg (SE, _("The day argument to DATE.YRDAY is not an integer.  "
  88                  "The result will be system-missing."));
  89       return SYSMIS;
  90     }
  91   else if (yd < 1 || yd > 366)
  92     {
  93       msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable "
  94                  "range of 1 to 366.  "
  95                  "The result will be system-missing."));
  96       return SYSMIS;
  97     }
  98   else
  99     {
 100       double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
 101       if (yr_1_1 != SYSMIS)
 102         return DAY_S * (yr_1_1 + yd - 1.);
 103       else
 104         return SYSMIS;
 105     }
 106 }
 107
 108 double
 109 expr_yrmoda (double year, double month, double day)
 110 {
 111   if (year >= 0 && year <= 99)
 112     year += 1900;
 113   else if (year != (int) year && year > 47516)
 114     {
 115       msg (SE, _("The year argument to YRMODA is greater than 47516.  "
 116                  "The result will be system-missing."));
 117       return SYSMIS;
 118     }
 119
 120   return expr_ymd_to_ofs (year, month, day);
 121 }
 122
 123 int
 124 compare_string (const struct substring *a, const struct substring *b)
 125 {
 126   size_t i;
 127
 128   for (i = 0; i < a->length && i < b->length; i++)
 129     if (a->string[i] != b->string[i])
 130       return a->string[i] < b->string[i] ? -1 : 1;
 131   for (; i < a->length; i++)
 132     if (a->string[i] != ' ')
 133       return 1;
 134   for (; i < b->length; i++)
 135     if (b->string[i] != ' ')
 136       return -1;
 137   return 0;
 138 }
 139
 140 size_t
 141 count_valid (double *d, size_t d_cnt)
 142 {
 143   size_t valid_cnt;
 144   size_t i;
 145
 146   valid_cnt = 0;
 147   for (i = 0; i < d_cnt; i++)
 148     valid_cnt += is_valid (d[i]);
 149   return valid_cnt;
 150 }
 151
 152 struct substring
 153 alloc_string (struct expression *e, size_t length)
 154 {
 155   struct substring s;
 156   s.length = length;
 157   s.string = pool_alloc (e->eval_pool, length);
 158   return s;
 159 }
 160
 161 struct substring
 162 copy_string (struct expression *e, const char *old, size_t length)
 163 {
 164   struct substring s = alloc_string (e, length);
 165   memcpy (s.string, old, length);
 166   return s;
 167 }
 168
 169 /* Returns the noncentral beta cumulative distribution function
 170    value for the given arguments.
 171
 172    FIXME: The accuracy of this function is not entirely
 173    satisfactory.  We only match the example values given in AS
 174    310 to the first 5 significant digits. */
 175 double
 176 ncdf_beta (double x, double a, double b, double lambda)
 177 {
 178   double c;
 179
 180   if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
 181     return SYSMIS;
 182
 183   c = lambda / 2.;
 184   if (lambda < 54.)
 185     {
 186       /* Algorithm AS 226. */
 187       double x0, a0, beta, temp, gx, q, ax, sumq, sum;
 188       double err_max = 2 * DBL_EPSILON;
 189       double err_bound;
 190       int iter_max = 100;
 191       int iter;
 192
 193       x0 = floor (c - 5.0 * sqrt (c));
 194       if (x0 < 0.)
 195         x0 = 0.;
 196       a0 = a + x0;
 197       beta = (gsl_sf_lngamma (a0)
 198               + gsl_sf_lngamma (b)
 199               - gsl_sf_lngamma (a0 + b));
 200       temp = gsl_sf_beta_inc (a0, b, x);
 201       gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
 202       if (a0 >= a)
 203         q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
 204       else
 205         q = exp (-c);
 206       ax = q * temp;
 207       sumq = 1. - q;
 208       sum = ax;
 209
 210       iter = 0;
 211       do
 212         {
 213           iter++;
 214           temp -= gx;
 215           gx = x * (a + b + iter - 1.) * gx / (a + iter);
 216           q *= c / iter;
 217           sumq -= q;
 218           ax = temp * q;
 219           sum += ax;
 220
 221           err_bound = (temp - gx) * sumq;
 222         }
 223       while (iter < iter_max && err_bound > err_max);
 224
 225       return sum;
 226     }
 227   else
 228     {
 229       /* Algorithm AS 310. */
 230       double m, m_sqrt;
 231       int iter, iter_lower, iter_upper, iter1, iter2, j;
 232       double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
 233       double err_bound;
 234       double err_max = 2 * DBL_EPSILON;
 235
 236       iter = 0;
 237
 238       m = floor (c + .5);
 239       m_sqrt = sqrt (m);
 240       iter_lower = m - 5. * m_sqrt;
 241       iter_upper = m + 5. * m_sqrt;
 242
 243       t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
 244       q = exp (t);
 245       r = q;
 246       psum = q;
 247       beta = (gsl_sf_lngamma (a + m)
 248               + gsl_sf_lngamma (b)
 249               - gsl_sf_lngamma (a + m + b));
 250       s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
 251       fx = gx = exp (s1);
 252       ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
 253       iter++;
 254       sum = q * temp;
 255       iter1 = m;
 256
 257       while (iter1 >= iter_lower && q >= err_max)
 258         {
 259           q = q * iter1 / c;
 260           iter++;
 261           gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
 262           iter1--;
 263           temp += gx;
 264           psum += q;
 265           sum += q * temp;
 266         }
 267
 268       t0 = (gsl_sf_lngamma (a + b)
 269             - gsl_sf_lngamma (a + 1.)
 270             - gsl_sf_lngamma (b));
 271       s0 = a * log (x) + b * log (1. - x);
 272
 273       s = 0.;
 274       for (j = 0; j < iter1; j++)
 275         {
 276           double t1;
 277           s += exp (t0 + s0 + j * log (x));
 278           t1 = log (a + b + j) - log (a + 1. + j) + t0;
 279           t0 = t1;
 280         }
 281
 282       err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
 283       q = r;
 284       temp = ftemp;
 285       gx = fx;
 286       iter2 = m;
 287       for (;;)
 288         {
 289           double ebd = err_bound + (1. - psum) * temp;
 290           if (ebd < err_max || iter >= iter_upper)
 291             break;
 292
 293           iter2++;
 294           iter++;
 295           q = q * c / iter2;
 296           psum += q;
 297           temp -= gx;
 298           gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
 299           sum += q * temp;
 300         }
 301
 302       return sum;
 303     }
 304 }
 305
 306 double
 307 cdf_bvnor (double x0, double x1, double r)
 308 {
 309   double z = x0 * x0 - 2. * r * x0 * x1 + x1 * x1;
 310   return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
 311 }
 312
 313 double
 314 idf_fdist (double P, double df1, double df2)
 315 {
 316   double temp = gslextras_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
 317   return temp * df2 / ((1. - temp) * df1);
 318 }
 319
 320 /*
 321  *  Mathlib : A C Library of Special Functions
 322  *  Copyright (C) 1998 Ross Ihaka
 323  *  Copyright (C) 2000 The R Development Core Team
 324  *
 325  *  This program is free software; you can redistribute it and/or
 326  *  modify
 327  *  it under the terms of the GNU General Public License as
 328  *  published by
 329  *  the Free Software Foundation; either version 2 of the
 330  *  License, or
 331  *  (at your option) any later version.
 332  *
 333  *  This program is distributed in the hope that it will be
 334  *  useful,
 335  *  but WITHOUT ANY WARRANTY; without even the implied warranty
 336  *  of
 337  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 338  *  GNU General Public License for more details.
 339  *
 340  *  You should have received a copy of the GNU General Public
 341  *  License
 342  *  along with this program; if not, write to the Free Software
 343  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 344  *  02110-1301 USA.
 345  */
 346
 347 /* Returns the density of the noncentral beta distribution with
 348    noncentrality parameter LAMBDA. */
 349 double
 350 npdf_beta (double x, double a, double b, double lambda)
 351 {
 352   if (lambda < 0. || a <= 0. || b <= 0.)
 353     return SYSMIS;
 354   else if (lambda == 0.)
 355     return gsl_ran_beta_pdf (x, a, b);
 356   else
 357     {
 358       double max_error = 2 * DBL_EPSILON;
 359       int max_iter = 200;
 360       double term = gsl_ran_beta_pdf (x, a, b);
 361       double lambda2 = 0.5 * lambda;
 362       double weight = exp (-lambda2);
 363       double sum = weight * term;
 364       double psum = weight;
 365       int k;
 366       for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
 367         weight *= lambda2 / k;
 368         term *= x * (a + b) / a;
 369         sum += weight * term;
 370         psum += weight;
 371         a += 1;
 372       }
 373       return sum;
 374     }
 375 }