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