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