LibRan  0.1
Pseudo-random number distribution generator
LRuinvcdf.c
Go to the documentation of this file.
1 
109 /*
110  * Copyright 2019 R.K. Owen, Ph.D.
111  * License see lgpl.md (Gnu Lesser General Public License)
112  */
113 #ifdef __cplusplus
114 extern "C" {
115 #endif
116 
117 #include "libran.h"
118 #include <math.h>
119 
120 #define Abs(a) (((a)<0)?-(a):(a))
121 #define Signabs(a,b) ((b)<0?-(a):(a))
122 #define Sign(a,b) ((b)<0?-Abs(a):Abs(a))
123 
131 int LRd_uinvcdf(LR_obj *o, double (*cdf)(double)) {
132  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
133  if (o->d != LR_double) {
134  o->errno = LRerr_BadDataType;
135  return o->errno;
136  }
137  aux->dcdf = cdf;
138 
139  return LRerr_OK;
140 }
141 
170 double LRd_zeroin(
171  double ax, double bx,
172  double U, double (*f)(double), double tol) {
173 
174  const double zero = 0.0, one = 1.0, two = 2.0, three = 3.0,
175  half = 0.5;
176  static double eps;
177  static char FIRST = (0 == 0);
178  double a,b,c,e,d,fa,fb,fc,tol1,xm,p,q,r,s,tmp1,tmp2;
179 
180  if (FIRST) {
181  eps = LR_dgetval("LR_DEPS");
182  FIRST = !FIRST;
183  }
184 /* initialization */
185  a = ax;
186  b = bx;
187  fa = (*f)(a) - U;
188  fb = (*f)(b) - U;
189 /* begin step */
190  c = a;
191  fc = fa;
192  d = b - a;
193  e = d;
194 
195  while (1) {
196 /* tighten interval */
197  if (Abs(fc) < Abs(fb)) {
198  a = b;
199  b = c;
200  c = a;
201  fa = fb;
202  fb = fc;
203  fc = fa;
204  }
205 /* convergence test */
206  tol1 = eps*Abs(b) + half*tol;
207  xm = half*(c - b);
208  if (Abs(xm) < tol1 || fb == zero) return (b); /* end */
209 
210  if (Abs(e) < tol1 || Abs(fa) <= Abs(fb)) {
211 /* bisection */
212  d = xm;
213  e = d;
214  } else {
215  if (a == c) {
216 /* linear interpolation */
217  s = fb/fa;
218  p = two*xm*s;
219  q = one - s;
220  } else {
221 /* inverse quadratic interpolation */
222  q = fa/fc;
223  r = fb/fc;
224  s = fb/fa;
225  p = s*(two*xm*q*(q - r) - (b - a)*(r - one));
226  q = (q - one)*(r - one)*(s - one);
227  }
228 /* adjusts sign */
229  if (p > zero) q = -q;
230  p = Abs(p);
231 
232  tmp1 = tol1*q;
233  tmp2 = half*e*q;
234  if ((two*p) >= (three*xm*q - Abs(tmp1))
235  || (p >= Abs(tmp2))) {
236 /* bisection */
237  d = xm;
238  e = d;
239  } else {
240 /* interpolation is acceptable */
241  e = d;
242  d = p/q;
243  }
244  }
245 /* complete step */
246  a = b;
247  fa = fb;
248  if (Abs(d) > tol1) b += d;
249  else b += Signabs(tol1,xm);
250  fb = (*f)(b) - U;
251 
252  if ((fb*Signabs(one,fc)) > zero) {
253 /* begin step again */
254  c = a;
255  fc = fa;
256  d = b - a;
257  e = d;
258  }
259  }
260 }
261 
269 double LRd_uinvcdf_RAN(LR_obj *o) {
270  double u, ax,bx, fax, fbx;
271  double zero = 0.0, one = 1.0, two = 2.0;
272  if (!o->aux) {
273  o->errno = LRerr_NoAuxiliaryObject;
274  return NAN;
275  }
276  if (!(((LR_uinvcdf *) o->aux)->dcdf)) {
277  o->errno = LRerr_BadAuxSetup;
278  return NAN;
279  }
280 
281  u = o->ud(o);
282  if (!isnan(o->a.d)) {
283  ax = o->a.d;
284  } else {
285  double scale = one;
286  ax = o->m.d;
287  fax = ((LR_uinvcdf *) o->aux)->dcdf(ax);
288  if (fax < zero) {
289  o->errno = LRerr_InvalidCDF;
290  return NAN;
291  }
292  while (fax > u) {
293  ax = o->m.d - scale*o->s.d;
294  fax = ((LR_uinvcdf *) o->aux)->dcdf(ax);
295  scale *= two;
296  }
297  }
298  if (!isnan(o->b.d)) {
299  bx = o->b.d;
300  } else {
301  double scale = one;
302  bx = o->m.d;
303  fbx = ((LR_uinvcdf *) o->aux)->dcdf(bx);
304  if (fbx > one) {
305  o->errno = LRerr_InvalidCDF;
306  return NAN;
307  }
308  while (fbx < u) {
309  bx = o->m.d + scale*o->s.d;
310  fbx = ((LR_uinvcdf *) o->aux)->dcdf(bx);
311  scale *= two;
312  }
313  }
314  return LRd_zeroin(ax, bx, u, ((LR_uinvcdf *) o->aux)->dcdf, zero);
315 }
316 
329 double LRd_uinvcdf_PDF(LR_obj *o, double x) {
330  static double sqeps = NAN, isqeps, nearzero;
331  double half = 0.5, zero = 0.0, one = 1.0, xp, xm, fp, fm, f0, ret;
332  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
333 
334  if (!aux) {
335  o->errno = LRerr_NoAuxiliaryObject;
336  return NAN;
337  }
338  if (!(aux->dcdf)) {
339  o->errno = LRerr_BadAuxSetup;
340  return NAN;
341  }
342 
343  if (isnan(sqeps)) {
344  sqeps = LR_dgetval("LR_DSQEPS");
345  isqeps = one/sqeps;
346  nearzero = sqeps*sqrt(sqeps);
347  }
348 
349  /* use three values around x */
350  if (-nearzero < x && x < nearzero) {
351  xp = sqeps;
352  xm = -sqeps;
353  } else {
354  xp = x*(one + sqeps);
355  xm = x*(one - sqeps);
356  }
357  fm = aux->dcdf(xm);
358  fp = aux->dcdf(xp);
359 
360  if (fp == one) {
361  f0 = aux->dcdf(x);
362  ret = (f0-fm)/(x-xm);
363  } else if (fm == zero) {
364  f0 = aux->dcdf(x);
365  ret = (fp-f0)/(xp-x);
366  } else {
367  ret = (fp-fm)/(xp-xm);
368  }
369  if (ret < zero) {
370  o->errno = LRerr_InvalidCDF;
371  return NAN;
372  }
373  return ret;
374 }
375 
384 double LRd_uinvcdf_CDF(LR_obj *o, double x) {
385  double zero = 0.0, one = 1.0, ret;
386  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
387 
388  if (!aux) {
389  o->errno = LRerr_NoAuxiliaryObject;
390  return NAN;
391  }
392  if (!(aux->dcdf)) {
393  o->errno = LRerr_BadAuxSetup;
394  return NAN;
395  }
396 
397  ret = aux->dcdf(x);
398 
399  if (ret < zero || ret > one) {
400  o->errno = LRerr_InvalidCDF;
401  return NAN;
402  }
403  return ret;
404 }
405 
413 int LRf_uinvcdf(LR_obj *o, float (*cdf)(float)) {
414  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
415  if (o->d != LR_float) {
416  o->errno = LRerr_BadDataType;
417  return o->errno;
418  }
419  aux->fcdf = cdf;
420 
421  return LRerr_OK;
422 }
423 
452 float LRf_zeroin(
453  float ax, float bx,
454  float U, float (*f)(float), float tol) {
455 
456  const float zero = 0.0, one = 1.0, two = 2.0, three = 3.0,
457  half = 0.5;
458  static float eps;
459  static char FIRST = (0 == 0);
460  float a,b,c,e,d,fa,fb,fc,tol1,xm,p,q,r,s,tmp1,tmp2;
461 
462  if (FIRST) {
463  eps = LR_fgetval("LR_FEPS");
464  FIRST = !FIRST;
465  }
466 /* initialization */
467  a = ax;
468  b = bx;
469  fa = (*f)(a) - U;
470  fb = (*f)(b) - U;
471 /* begin step */
472  c = a;
473  fc = fa;
474  d = b - a;
475  e = d;
476 
477  while (1) {
478 /* tighten interval */
479  if (Abs(fc) < Abs(fb)) {
480  a = b;
481  b = c;
482  c = a;
483  fa = fb;
484  fb = fc;
485  fc = fa;
486  }
487 /* convergence test */
488  tol1 = eps*Abs(b) + half*tol;
489  xm = half*(c - b);
490  if (Abs(xm) < tol1 || fb == zero) return (b); /* end */
491 
492  if (Abs(e) < tol1 || Abs(fa) <= Abs(fb)) {
493 /* bisection */
494  d = xm;
495  e = d;
496  } else {
497  if (a == c) {
498 /* linear interpolation */
499  s = fb/fa;
500  p = two*xm*s;
501  q = one - s;
502  } else {
503 /* inverse quadratic interpolation */
504  q = fa/fc;
505  r = fb/fc;
506  s = fb/fa;
507  p = s*(two*xm*q*(q - r) - (b - a)*(r - one));
508  q = (q - one)*(r - one)*(s - one);
509  }
510 /* adjusts sign */
511  if (p > zero) q = -q;
512  p = Abs(p);
513 
514  tmp1 = tol1*q;
515  tmp2 = half*e*q;
516  if ((two*p) >= (three*xm*q - Abs(tmp1))
517  || (p >= Abs(tmp2))) {
518 /* bisection */
519  d = xm;
520  e = d;
521  } else {
522 /* interpolation is acceptable */
523  e = d;
524  d = p/q;
525  }
526  }
527 /* complete step */
528  a = b;
529  fa = fb;
530  if (Abs(d) > tol1) b += d;
531  else b += Signabs(tol1,xm);
532  fb = (*f)(b) - U;
533 
534  if ((fb*Signabs(one,fc)) > zero) {
535 /* begin step again */
536  c = a;
537  fc = fa;
538  d = b - a;
539  e = d;
540  }
541  }
542 }
543 
551 float LRf_uinvcdf_RAN(LR_obj *o) {
552  float u, ax,bx, fax, fbx;
553  float zero = 0.0, one = 1.0, two = 2.0;
554  if (!o->aux) {
555  o->errno = LRerr_NoAuxiliaryObject;
556  return NAN;
557  }
558  if (!(((LR_uinvcdf *) o->aux)->fcdf)) {
559  o->errno = LRerr_BadAuxSetup;
560  return NAN;
561  }
562 
563  u = o->ud(o);
564  if (!isnan(o->a.f)) {
565  ax = o->a.f;
566  } else {
567  float scale = one;
568  ax = o->m.f;
569  fax = ((LR_uinvcdf *) o->aux)->fcdf(ax);
570  if (fax < zero) {
571  o->errno = LRerr_InvalidCDF;
572  return NAN;
573  }
574  while (fax > u) {
575  ax = o->m.f - scale*o->s.f;
576  fax = ((LR_uinvcdf *) o->aux)->fcdf(ax);
577  scale *= two;
578  }
579  }
580  if (!isnan(o->b.f)) {
581  bx = o->b.f;
582  } else {
583  float scale = one;
584  bx = o->m.f;
585  fbx = ((LR_uinvcdf *) o->aux)->fcdf(bx);
586  if (fbx > one) {
587  o->errno = LRerr_InvalidCDF;
588  return NAN;
589  }
590  while (fbx < u) {
591  bx = o->m.f + scale*o->s.f;
592  fbx = ((LR_uinvcdf *) o->aux)->fcdf(bx);
593  scale *= two;
594  }
595  }
596  return LRf_zeroin(ax, bx, u, ((LR_uinvcdf *) o->aux)->fcdf, zero);
597 }
598 
611 float LRf_uinvcdf_PDF(LR_obj *o, float x) {
612  static float sqeps = NAN, isqeps, nearzero;
613  float half = 0.5, zero = 0.0, one = 1.0, xp, xm, fp, fm, f0, ret;
614  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
615 
616  if (!aux) {
617  o->errno = LRerr_NoAuxiliaryObject;
618  return NAN;
619  }
620  if (!(aux->fcdf)) {
621  o->errno = LRerr_BadAuxSetup;
622  return NAN;
623  }
624 
625  if (isnan(sqeps)) {
626  sqeps = LR_fgetval("LR_FSQEPS");
627  isqeps = one/sqeps;
628  nearzero = sqeps*sqrtf(sqeps);
629  }
630 
631  /* use three values around x */
632  if (-nearzero < x && x < nearzero) {
633  xp = sqeps;
634  xm = -sqeps;
635  } else {
636  xp = x*(one + sqeps);
637  xm = x*(one - sqeps);
638  }
639  fm = aux->fcdf(xm);
640  fp = aux->fcdf(xp);
641 
642  if (fp == one) {
643  f0 = aux->fcdf(x);
644  ret = (f0-fm)/(x-xm);
645  } else if (fm == zero) {
646  f0 = aux->fcdf(x);
647  ret = (fp-f0)/(xp-x);
648  } else {
649  ret = (fp-fm)/(xp-xm);
650  }
651  if (ret < zero) {
652  o->errno = LRerr_InvalidCDF;
653  return NAN;
654  }
655  return ret;
656 }
657 
666 float LRf_uinvcdf_CDF(LR_obj *o, float x) {
667  float zero = 0.0, one = 1.0, ret;
668  LR_uinvcdf *aux = (LR_uinvcdf *) o->aux;
669 
670  if (!aux) {
671  o->errno = LRerr_NoAuxiliaryObject;
672  return NAN;
673  }
674  if (!(aux->fcdf)) {
675  o->errno = LRerr_BadAuxSetup;
676  return NAN;
677  }
678 
679  ret = aux->fcdf(x);
680 
681  if (ret < zero || ret > one) {
682  o->errno = LRerr_InvalidCDF;
683  return NAN;
684  }
685  return ret;
686 }
687 
688 #ifdef __cplusplus
689 }
690 #endif
LR_obj::b
LR_val b
Definition: libran.h:139
LRerr_OK
#define LRerr_OK
Definition: libran.h:32
LR_obj::d
LR_data_type d
Definition: libran.h:137
LR_fgetval
float LR_fgetval(char *str)
LR_fgetval(char *val) - return the configure value.
Definition: urand.c:464
LRf_zeroin
float LRf_zeroin(float ax, float bx, float U, float(*f)(float), float tol)
LRf_zeroin() Routine numerically finds solution to f(x)-U = 0 specialized to CDF()s.
Definition: LRuinvcdf.c:452
LRerr_BadAuxSetup
#define LRerr_BadAuxSetup
Definition: libran.h:38
LR_obj::aux
void * aux
Definition: libran.h:168
LR_dgetval
double LR_dgetval(char *str)
LR_dgetval(char *val) - return the configure value.
Definition: urand.c:483
LRf_uinvcdf_CDF
float LRf_uinvcdf_CDF(LR_obj *o, float x)
LRf_uinvcdf_CDF(LR_obj *o, float x) - float User supplied cumulative distribution function...
Definition: LRuinvcdf.c:666
LRd_zeroin
double LRd_zeroin(double ax, double bx, double U, double(*f)(double), double tol)
LRd_zeroin() Routine numerically finds solution to f(x)-U = 0 specialized to CDF()s.
Definition: LRuinvcdf.c:170
LRd_uinvcdf_RAN
double LRd_uinvcdf_RAN(LR_obj *o)
LRd_uinvcdf_RAN(LR_obj *o) - double random variate via inverse method of the UserCDF() fn...
Definition: LRuinvcdf.c:269
LRf_uinvcdf
int LRf_uinvcdf(LR_obj *o, float(*cdf)(float))
LRf_uinvcdf() - set the user defined CDF for this variate distribution.
Definition: LRuinvcdf.c:413
LRf_uinvcdf_RAN
float LRf_uinvcdf_RAN(LR_obj *o)
LRf_uinvcdf_RAN(LR_obj *o) - float random variate via inverse method of the UserCDF() fn...
Definition: LRuinvcdf.c:551
LRerr_BadDataType
#define LRerr_BadDataType
Definition: libran.h:34
LR_obj::ud
double(* ud)(LR_obj *)
Definition: libran.h:154
LRd_uinvcdf_PDF
double LRd_uinvcdf_PDF(LR_obj *o, double x)
LRd_uinvcdf_PDF(LR_obj *o, double x) - double probability distribution function approximated from the...
Definition: LRuinvcdf.c:329
LR_uinvcdf
A special object for using a user defined CDF.
Definition: libran.h:238
LR_obj::s
LR_val s
Definition: libran.h:141
LRd_uinvcdf
int LRd_uinvcdf(LR_obj *o, double(*cdf)(double))
LRd_uinvcdf() - set the user defined CDF for this variate distribution.
Definition: LRuinvcdf.c:131
LR_obj::errno
int errno
Definition: libran.h:169
LR_obj::a
LR_val a
Definition: libran.h:138
LR_val::d
double d
Definition: libran.h:87
LR_val::f
float f
Definition: libran.h:86
LRerr_NoAuxiliaryObject
#define LRerr_NoAuxiliaryObject
Definition: libran.h:36
LRf_uinvcdf_PDF
float LRf_uinvcdf_PDF(LR_obj *o, float x)
LRf_uinvcdf_PDF(LR_obj *o, float x) - float probability distribution function approximated from the u...
Definition: LRuinvcdf.c:611
LRd_uinvcdf_CDF
double LRd_uinvcdf_CDF(LR_obj *o, double x)
LRd_uinvcdf_CDF(LR_obj *o, double x) - double User supplied cumulative distribution function...
Definition: LRuinvcdf.c:384
libran.h
The LibRan common header file.
LRerr_InvalidCDF
#define LRerr_InvalidCDF
Definition: libran.h:46
LR_obj
the fundamental LibRan random variate distribution object
Definition: libran.h:134
LR_obj::m
LR_val m
Definition: libran.h:140
Generated by   doxygen 1.8.13