LibRan  0.1
Pseudo-random number distribution generator
gsn.c
Go to the documentation of this file.
1 
94 #include <stdio.h>
95 #include <limits.h>
96 
97 #define NINT 15
98 
103 typedef long LINT;
104 
105 /* poly order -> , interval downward v */
106 
111 LINT pn[NINT][NINT] = {
112  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
113  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
114  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
115  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
116  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
117  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
118  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
119  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
120  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
121  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
122  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
123  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
124  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
125  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
126  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
127 };
128 
133 LINT pd[NINT][NINT] = {
134  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
135  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
136  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
137  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
138  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
139  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
140  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
141  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
142  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
143  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
144  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
145  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
146  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
147  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
148  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
149 };
150 
155 LINT cn[NINT][NINT] = {
156  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
157  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
158  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
159  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
160  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
161  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
162  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
163  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
164  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
165  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
166  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
167  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
168  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
169  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
170  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
171 };
172 
177 LINT cd[NINT][NINT] = {
178  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
179  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
180  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
181  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
182  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
183  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
184  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
185  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
186  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
187  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
188  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
189  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
190  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
191  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
192  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
193 };
194 
198 LINT igcd(LINT a, LINT b) {
199  LINT t;
200  if (a == 0 || b == 0) return 0;
201  if (a < 0) a = -a;
202  if (b < 0) b = -b;
203  if (a < b) { /* swap */ t = a; a = b; b = t; }
204  while ((t = a%b) > 1) { a = b; b = t; }
205  if (t) return t; else return b;
206 }
207 
211 LINT ifac(int n) {
212  LINT i = 1, t = 1;
213  if (n == 0 || n == 1) return t;
214 
215  do {
216  i++;
217  t *= i;
218  } while (i < n);
219 
220  return t;
221 }
222 
226 LINT icomb(int n, int k) {
227  return ifac(n) / ifac(k) / ifac(n-k);
228 }
229 
242 void iadd(LINT an, LINT ad, LINT bn, LINT bd, LINT x, LINT *on, LINT *od) {
243  LINT xx;
244 
245  if (an == 0 || ad == 0) ad = 1;
246  if (bn == 0 || bd == 0) bd = 1;
247 
248  *on = an*bd + bn*ad*x;
249  *od = ad*bd;
250  xx = igcd(*on, *od);
251  if (xx != 0 && xx != 1) {
252  *on /= xx;
253  *od /= xx;
254  }
255  /* set denominator to 1 if numerator is 0 */
256  if (*on == 0)
257  *od = 1;
258 }
259 
272 void imult(LINT an, LINT ad, LINT bn, LINT bd, LINT x, LINT *on, LINT *od) {
273  LINT xx;
274 
275  if (an == 0 || ad == 0) ad = 1;
276  if (bn == 0 || bd == 0) bd = 1;
277 
278  *on = an*bn*x;
279  *od = ad*bd;
280  xx = igcd(*on, *od);
281  if (xx != 0 && xx != 1) {
282  *on /= xx;
283  *od /= xx;
284  }
285  /* set denominator to 1 if numerator is 0 */
286  if (*on == 0)
287  *od = 1;
288 }
289 
298 void showdf(char *DFstr, int order, LINT xn[][NINT], LINT xd[][NINT]) {
299  printf("%s : %d >>> \n", DFstr, order);
300  for (int j = 0; j <= order; j++) {
301  printf("% 2d = ",j);
302  for (int i = 0; i <= order; i++) {
303  if (xn[j][i])
304  printf("%+ld/%ld*x^%d",xn[j][i],xd[j][i],i);
305  }
306  printf("\n");
307  }
308 }
309 
317 void writedf(char *DFstr, int order, LINT xn[][NINT]) {
318  LINT maxcoef = 0;
319 
320  printf("long %s%d[%d][%d] = {\n", DFstr, order, order+1, order+1);
321  for (int j = 0; j <= order; j++) {
322  printf(" {");
323  for (int i = 0; i <= order; i++) {
324  printf("%ldl",xn[j][i]);
325  if (i < order) printf(",");
326  LINT cc = (xn[j][i] < 0 ? -xn[j][i] : xn[j][i]);
327  if (cc > maxcoef)
328  maxcoef = cc;
329  }
330  printf("}");
331  if (j < order) printf(",");
332  printf("\n");
333  }
334  printf("};\n");
335  printf("/* integer type must support max value = %ld */\n\n", maxcoef);
336 }
337 
346 void newcdf(int nn) {
347  extern LINT pn[NINT][NINT], pd[NINT][NINT],
348  cn[NINT][NINT], cd[NINT][NINT];
349 
350  LINT ss = 1;
351  for (int i = 0; i <= nn; i++) {
352  /* do (x-i)^n for all intervals > i */
353  LINT coef = ss*icomb((LINT) nn, (LINT) i);
354  ss *= -1;
355  for (int j = i; j <= nn; j++) {
356  /* binomial expansion*/
357  cn[j][nn] += coef;
358  LINT s = (LINT) (-i);
359  for (int k = 1; k <= nn; k++) {
360  cn[j][nn-k]+=(s*coef*icomb((LINT) nn,(LINT) k));
361  s *= (LINT) (-i);
362  }
363  }
364  }
365  /* normalize with 1/n! */
366  LINT coef = ifac(nn);
367  LINT cc, maxint = 0, maxcoef = 0;
368  for (int j = 0; j <= nn; j++) {
369  for (int i = 0; i <= nn; i++) {
370  cc = (cn[j][i] < 0 ? -cn[j][i] : cn[j][i]);
371  if (cc > maxint)
372  maxint = cc;
373  imult(cn[j][i],cd[j][i],(LINT) 1,coef,(LINT) 1,
374  &(cn[j][i]),&(cd[j][i]));
375  cc = (cn[j][i] < 0 ? -cn[j][i] : cn[j][i]);
376  if (cc > maxcoef)
377  maxcoef = cc;
378  cc = (cd[j][i] < 0 ? -cd[j][i] : cd[j][i]);
379  if (cc > maxcoef)
380  maxcoef = cc;
381  }
382  }
383  printf("max int = %ld\nmax coef = %ld\n"
384  "INT_MAX = %d\nLONG_MAX = %ld\nLLONG_MAX = %lld\n"
385  "%d! = %ld\n",
386  maxint, maxcoef,
387  INT_MAX, LONG_MAX, LLONG_MAX, nn, ifac(nn));
388 }
389 
399 void newpdf(int nn) {
400  extern LINT pn[NINT][NINT], pd[NINT][NINT],
401  cn[NINT][NINT], cd[NINT][NINT];
402 
403  for (int j = 0; j <= nn; j++) {
404  for (int i = 1; i <= nn; i++) {
405  imult((LINT) i,(LINT) 1,cn[j][i],cd[j][i],(LINT) 1,
406  &(pn[j][i-1]), &(pd[j][i-1]));
407  }
408  }
409 }
410 
440 void transint(int nn, LINT xn[][NINT], LINT xd[][NINT]) {
441  LINT tmpn[NINT], tmpd[NINT], zero = 0, one = 1;
442  for (int j = 0; j <= nn; j++) { /* interval */
443  /* zero tmp */
444  for (int i = 0; i <= nn; i++) {
445  tmpn[i] = zero;
446  tmpd[i] = one;
447  }
448 
449  for (int m = 0; m <= nn; m++) {
450  /* use Horners rule and binomial expansion*/
451  LINT ct = j;
452  LINT cn = icomb(nn,nn-m)*xn[j][nn], cd = xd[j][nn];
453 
454  for (int k = nn-1; k >= m; k--) {
455 /*
456  printf("icomb(%d,%d) = %ld\txn = %ld\tcn = %ld\n",
457  k,k-m, icomb(k,k-m),xn[j][k], cn);
458 */
459  iadd(icomb(k,k-m)*xn[j][k], xd[j][k], cn, cd,
460  ct, &cn, &cd);
461  }
462  tmpn[m] = cn;
463  tmpd[m] = cd;
464  }
465  /* transfer results back */
466  for (int i = 0; i <= nn; i++) {
467  xn[j][i] = tmpn[i];
468  xd[j][i] = tmpd[i];
469  }
470  }
471  if (0) {
472 /* This code bit can be enabled for checking that the coefficient
473  * values make sense. On each local interval sum of the coefs is the
474  * same as evaluating the function where x=1 this
475  * should give the a_0 coef on the next interval or the x=0 evaluation.
476  */
477  for (int j = 0; j <= nn; j++) { /* interval */
478  /* sum coefs (same as x=1 on local interval */
479  /* should equal the constant on next interval */
480  LINT sn = zero, sd = one;
481  for (int i = 0; i <= nn; i++) {
482  iadd(sn,sd, xn[j][i], xd[j][i], one, &sn, &sd);
483  }
484  printf("*** %d %ld/%ld\n", j, sn, sd);
485  if (j < nn)
486  printf("+++ %d %ld/%ld\n", j+1,
487  xn[j+1][0], xd[j+1][0]);
488  }
489  }
490 }
491 
492 int main() {
493  int nn = 12;
494  newcdf(nn);
495  showdf("CDF",nn, cn, cd);
496  newpdf(nn);
497  showdf("PDF",nn, pn, pd);
498 if (1) {
499  transint(nn, cn, cd);
500  showdf("CDF",nn, cn, cd);
501  transint(nn, pn, pd);
502  showdf("PDF",nn, pn, pd);
503 }
504 
505 if (0) {
506  writedf("gscdfn", nn, cn);
507  writedf("gscdfd", nn, cd);
508  writedf("gspdfn", nn, pn);
509  writedf("gspdfd", nn, pd);
510 }
511 }
pd
LINT pd[NINT][NINT]
PDF denominator.
Definition: gsn.c:133
newpdf
void newpdf(int nn)
void newpdf(nn) calculate the new PDF of order nn by differentiating the CDF.
Definition: gsn.c:399
cn
LINT cn[NINT][NINT]
CDF numerator.
Definition: gsn.c:155
pn
LINT pn[NINT][NINT]
PDF numerator.
Definition: gsn.c:111
imult
void imult(LINT an, LINT ad, LINT bn, LINT bd, LINT x, LINT *on, LINT *od)
void imult(an,ad,bn,bd, x, *on, *od) return a*b*x where a,b are represented as rationals ...
Definition: gsn.c:272
igcd
LINT igcd(LINT a, LINT b)
LINT igcd(LINT a, LINT b) return Greatest Common Denominator of a & b.
Definition: gsn.c:198
newcdf
void newcdf(int nn)
void newcdf(nn) calculate the new CDF of order nn
Definition: gsn.c:346
LINT
long LINT
long integer type used for calculations.
Definition: gsn.c:103
icomb
LINT icomb(int n, int k)
LINT icomb(LINT n, LINT k) return binomial factor n!/[k!*(n-k)!].
Definition: gsn.c:226
writedf
void writedf(char *DFstr, int order, LINT xn[][NINT])
void writedf(DFstr, order, xn[][NINT]) output array as a C declaration
Definition: gsn.c:317
cd
LINT cd[NINT][NINT]
CDF denominator.
Definition: gsn.c:177
showdf
void showdf(char *DFstr, int order, LINT xn[][NINT], LINT xd[][NINT])
void showdf(DFstr, order, xn[][NINT], xd[][NINT]) display the rational values of &#39;x&#39;.
Definition: gsn.c:298
ifac
LINT ifac(int n)
LINT ifac(LINT n) return the factorial of n (n!)
Definition: gsn.c:211
transint
void transint(int nn, LINT xn[][NINT], LINT xd[][NINT])
void transint(nn, xn[][NINT], xd[]{NINT]) translate each interval to [0,1)
Definition: gsn.c:440
iadd
void iadd(LINT an, LINT ad, LINT bn, LINT bd, LINT x, LINT *on, LINT *od)
void iadd(an,ad,bn,bd, x, *on, *od) return a + b*x where a,b are represented as rationals ...
Definition: gsn.c:242
Generated by   doxygen 1.8.13