1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: fmarsa.c
5 * Environment: ANSI C
6 *
7 * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
8 * e-mail: lecuyer@iro.umontreal.ca
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted without a fee for private, research,
13 * academic, or other non-commercial purposes.
14 * Any use of this software in a commercial environment requires a
15 * written licence from the copyright owner.
16 *
17 * Any changes made to this package must be clearly identified as such.
18 *
19 * In scientific publications which used this software, a reference to it
20 * would be appreciated.
21 *
22 * Redistributions of source code must retain this copyright notice
23 * and the following disclaimer.
24 *
25 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
27 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
28 *
29 \*************************************************************************/
30
31
32 #include "util.h"
33 #include "gofs.h"
34 #include "num.h"
35
36 #include "fmarsa.h"
37 #include "fcho.h"
38 #include "ffam.h"
39 #include "fres.h"
40 #include "ftab.h"
41 #include "smarsa.h"
42 #include "unif01.h"
43
44 #include <string.h>
45 #include <limits.h>
46 #include <math.h>
47
48 long fmarsa_Maxn = 1024 * 1024 * 32;
49 long fmarsa_MaxL = 1024 * 4;
50
51
52
53
54 /*------------------------------ Functions --------------------------------*/
55
56
InitRes2(ffam_Fam * fam,fmarsa_Res2 * res,int N,int Nr,int j1,int j2,int jstep,char * name1,char * name2)57 static void InitRes2 (
58 ffam_Fam *fam,
59 fmarsa_Res2 *res, /* Results holder */
60 int N, /* Number of replications */
61 int Nr,
62 int j1, int j2, int jstep,
63 char *name1,
64 char *name2
65 )
66 /*
67 * Initializes the fmarsa_Res2 structure
68 */
69 {
70 fres_InitCont (fam, res->GCD, N, Nr, j1, j2, jstep, name1);
71 fres_InitCont (fam, res->NumIter, N, Nr, j1, j2, jstep, name2);
72 }
73
74
75 /*-------------------------------------------------------------------------*/
76
fmarsa_CreateRes2(void)77 fmarsa_Res2 * fmarsa_CreateRes2 (void)
78 {
79 fmarsa_Res2 *res;
80 res = util_Malloc (sizeof (fmarsa_Res2));
81 res->NumIter = fres_CreateCont ();
82 res->GCD = fres_CreateCont ();
83 return res;
84 }
85
86
87 /*-------------------------------------------------------------------------*/
88
fmarsa_DeleteRes2(fmarsa_Res2 * res)89 void fmarsa_DeleteRes2 (fmarsa_Res2 *res)
90 {
91 if (res == NULL)
92 return;
93 fres_DeleteCont (res->GCD);
94 fres_DeleteCont (res->NumIter);
95 util_Free (res);
96 }
97
98
99 /*=========================================================================*/
100
PrintHead(char * test,ffam_Fam * fam,long N,long n,int r,int s,int L,int t,int p,int Nr,int j1,int j2,int jstep)101 static void PrintHead (char *test, ffam_Fam * fam,
102 long N, long n, int r, int s, int L, int t, int p,
103 int Nr, int j1, int j2, int jstep)
104 {
105 printf
106 ("\n\n================================================================\n");
107 printf ("Family: %s\n\n", fam->name);
108 printf ("Test: %s\n", test);
109 printf (" N = %ld,", N);
110 if (n)
111 printf (" n = %ld,", n);
112 printf (" r = %d,", r);
113 if (s)
114 printf (" s = %d,", s);
115 if (L)
116 printf (" L = %d", L);
117 if (t)
118 printf (" t = %d,", t);
119 if (p)
120 printf (" p = %d", p);
121 printf ("\n Nr = %d, j1 = %d, j2 = %d, jstep = %d\n\n",
122 Nr, j1, j2, jstep);
123 }
124
125
126 /*=========================================================================*/
127
CheckParamMat(int prec,void * cho,long * pn,int * pr,int * ps,long * pL,long LMin,int i,int j)128 static int CheckParamMat (int prec, void *cho,
129 long *pn, int *pr, int *ps, long *pL, long LMin, int i, int j)
130 /*
131 * Set the values of the parameters for the test. If a parameter is < 0,
132 * will call a choose function to set it. Otherwise, will accept it as is.
133 * Returns 0 if parameters are ok for the test, returns -1 if the test
134 * should not be done for these parameters.
135 */
136 {
137 fcho_Cho2 *cho2 = cho;
138 fcho_Cho *chon;
139 fcho_Cho *choL;
140
141 util_Assert (cho, "fmarsa: cho is NULL");
142 chon = cho2->Chon;
143 choL = cho2->Chop2;
144 if (*pn < 0) {
145 util_Assert (chon, "fmarsa: n < 0 and chon is NULL");
146 *pn = chon->Choose (chon->param, i, j);
147
148 if (*pn <= 3.0 * gofs_MinExpected) {
149 printf ("n is too small\n\n");
150 return -1;
151 }
152 if (*pn > fmarsa_Maxn) {
153 printf ("n > %2ld\n\n", fmarsa_Maxn);
154 return -1;
155 }
156 }
157
158 *ps = fcho_Chooses (*pr, *ps, prec);
159 if (*ps <= 0)
160 return -1;
161
162 if (*pL < 0) {
163 util_Assert (choL, "fmarsa: L < 0 and chop2 is NULL");
164 *pL = choL->Choose (choL->param, i, j);
165
166 if (*pL <= LMin) {
167 printf ("L is too small\n\n");
168 return -1;
169 }
170 if (*pL > fmarsa_MaxL) {
171 printf ("L > %2ld\n\n", fmarsa_MaxL);
172 return -1;
173 }
174 }
175 return 0;
176 }
177
178
179 /*=========================================================================*/
180
181
TabMatrixR(ffam_Fam * fam,void * res1,void * cho,void * par1,int i,int j,int irow,int icol)182 static void TabMatrixR (ffam_Fam * fam, void *res1, void *cho,
183 void *par1, int i, int j, int irow, int icol)
184 {
185 int r, s;
186 long N, n, L;
187 const long *Par = par1;
188 fres_Cont *fres = res1;
189 sres_Chi2 *sres;
190
191 N = Par[0];
192 n = Par[1];
193 r = Par[2];
194 s = Par[3];
195 L = Par[4];
196
197 if (CheckParamMat (fam->Resol[irow], cho, &n, &r, &s, &L, 1, i, j))
198 return;
199
200 sres = sres_CreateChi2 ();
201 smarsa_MatrixRank (fam->Gen[irow], sres, N, n, r, s, L, L);
202 fres_FillTableEntryC (fres, sres->pVal2, N, irow, icol);
203 sres_DeleteChi2 (sres);
204 }
205
206
207 /*-------------------------------------------------------------------------*/
208
fmarsa_MatrixR1(ffam_Fam * fam,fres_Cont * res,fcho_Cho2 * cho,long N,long n,int r,int s,int L,int Nr,int j1,int j2,int jstep)209 void fmarsa_MatrixR1 (ffam_Fam * fam, fres_Cont * res, fcho_Cho2 * cho,
210 long N, long n, int r, int s, int L, int Nr, int j1, int j2, int jstep)
211 {
212 long Par[5] = { 0 };
213 lebool localRes;
214
215 Par[0] = N;
216 Par[1] = n;
217 Par[2] = r;
218 Par[3] = s;
219 Par[4] = L;
220 if (res == NULL) {
221 localRes = TRUE;
222 res = fres_CreateCont ();
223 } else
224 localRes = FALSE;
225
226 util_Assert (n < 0 || L < 0,
227 "fmarsa_MatrixR1: Either n or L must be < 0" );
228 PrintHead ("fmarsa_MatrixR1", fam, N, n, r, s, L, 0, 0, Nr, j1, j2,
229 jstep);
230 fres_InitCont (fam, res, N, Nr, j1, j2, jstep, "fmarsa_MatrixR1");
231 ftab_MakeTables (fam, res, cho, Par, TabMatrixR, Nr, j1, j2, jstep);
232 fres_PrintCont (res);
233 if (localRes)
234 fres_DeleteCont (res);
235 }
236
237
238 /*=========================================================================*/
239
WriteBirthEC(void * vpar,long junk1,long junk2)240 static void WriteBirthEC (void *vpar, long junk1, long junk2)
241 {
242 double *Par = vpar;
243 double EC = Par[2];
244 printf ("Choose d such that EC = %f\n\n", EC);
245 }
246
247 /*-------------------------------------------------------------------------*/
248
ChooseBirthEC(void * vpar,long n,long junk)249 static double ChooseBirthEC (void *vpar, long n, long junk)
250 {
251 double *Par = vpar;
252 long N = Par[0];
253 int t = Par[1];
254 double EC = Par[2];
255 long d;
256 double k, dr;
257 double Mu;
258
259 WriteBirthEC (vpar, 0, 0);
260 k = (N * (double) n * n * n) / (4.0 * EC);
261 if (k >= smarsa_Maxk) {
262 printf ("k >= %2.0f\n\n", smarsa_Maxk);
263 return -1.0;
264 }
265 d = dr = pow (k, 1.0 / t);
266 if (dr > LONG_MAX) {
267 printf ("d > LONG_MAX\n\n");
268 return -1.0;
269 }
270
271 k = pow ((double) d, (double) t);
272 Mu = N * (double) n * n * n / (4.0 * k);
273 if (8.0 * Mu > sqrt (sqrt (k))) {
274 printf ("8 EC > k^(1/4)\n\n");
275 return -1.0;
276 }
277 return (double) d;
278 }
279
280 /*-------------------------------------------------------------------------*/
281
fmarsa_CreateBirthEC(long N,int t,double EC)282 fcho_Cho *fmarsa_CreateBirthEC (long N, int t, double EC)
283 {
284 fcho_Cho *cho;
285 double *Par;
286
287 cho = util_Malloc (sizeof (fcho_Cho));
288 Par = util_Calloc (3, sizeof (double));
289 Par[0] = N;
290 Par[1] = t;
291 Par[2] = EC;
292 cho->param = Par;
293 cho->Write = WriteBirthEC;
294 cho->Choose = ChooseBirthEC;
295 cho->name = util_Calloc (2, sizeof (char));
296 strcpy (cho->name, "d");
297 return cho;
298 }
299
300 /*-------------------------------------------------------------------------*/
301
fmarsa_DeleteBirthEC(fcho_Cho * cho)302 void fmarsa_DeleteBirthEC (fcho_Cho * cho)
303 {
304 if (NULL == cho)
305 return;
306 cho->name = util_Free (cho->name);
307 cho->param = util_Free (cho->param);
308 util_Free (cho);
309 }
310
311
312 /*=========================================================================*/
313
CheckParamBirth(int prec,void * cho,long * pn,int * pr,long * pd,int i,int j)314 static int CheckParamBirth (int prec, void *cho,
315 long *pn, int *pr, long *pd, int i, int j)
316 /*
317 * Set the values of the parameters for the test.
318 * Returns 0 if parameters are ok for the test, returns -1 if the test
319 * should not be done for these parameters.
320 */
321 {
322 fcho_Cho2 *cho2 = cho;
323 fcho_Cho *chon;
324 fcho_Cho *chod;
325 int s;
326
327 util_Assert (cho, "fmarsa: cho is NULL");
328 chon = cho2->Chon;
329 chod = cho2->Chop2;
330 util_Assert (chon, "fmarsa: chon is NULL");
331 *pn = chon->Choose (chon->param, i, j);
332 if (*pn > fmarsa_Maxn) {
333 printf ("n > %2ld\n\n", fmarsa_Maxn);
334 return -1;
335 }
336
337 util_Assert (chod, "fmarsa: chop2 is NULL");
338 *pd = chod->Choose (chod->param, *pn, 0);
339 if (*pd <= 1.0)
340 return -1;
341
342 s = num_Log2 ((double) *pd);
343 if (*pr + s > prec) {
344 printf ("r + Lg(d) > Resolution of generator\n\n");
345 return -1;
346 }
347
348 return 0;
349 }
350
351
352 /*=========================================================================*/
353
TabBirthdayS(ffam_Fam * fam,void * vres,void * cho,void * vpar,int i,int j,int irow,int icol)354 static void TabBirthdayS (ffam_Fam * fam, void *vres, void *cho,
355 void *vpar, int i, int j, int irow, int icol)
356 {
357 int r, t, p;
358 long N, n, d;
359 const long *Par = vpar;
360 fres_Poisson *fres = vres;
361 sres_Poisson *sres;
362
363 N = Par[0];
364 r = Par[1];
365 t = Par[2];
366 p = Par[3];
367
368 if (CheckParamBirth (fam->Resol[irow], cho, &n, &r, &d, i, j))
369 return;
370
371 sres = sres_CreatePoisson ();
372 smarsa_BirthdaySpacings (fam->Gen[irow], sres, N, n, r, d, t, p);
373 fres_FillTableEntryPoisson (fres, sres->Mu, sres->sVal2, sres->pLeft,
374 sres->pRight, sres->pVal2, irow, icol);
375 sres_DeletePoisson (sres);
376 }
377
378
379 /*-------------------------------------------------------------------------*/
380
fmarsa_BirthdayS1(ffam_Fam * fam,fres_Poisson * res,fcho_Cho2 * cho,long N,int r,int t,int p,int Nr,int j1,int j2,int jstep)381 void fmarsa_BirthdayS1 (ffam_Fam * fam, fres_Poisson * res, fcho_Cho2 * cho,
382 long N, int r, int t, int p, int Nr, int j1, int j2, int jstep)
383 {
384 long Par[4] = { 0 };
385 lebool localRes;
386
387 Par[0] = N;
388 Par[1] = r;
389 Par[2] = t;
390 Par[3] = p;
391
392 if (res == NULL) {
393 localRes = TRUE;
394 res = fres_CreatePoisson ();
395 } else
396 localRes = FALSE;
397
398 PrintHead ("fmarsa_BirthdayS1",
399 fam, N, 0, r, 0, 0, t, p, Nr, j1, j2, jstep);
400 fres_InitPoisson (fam, res, Nr, j1, j2, jstep, "fmarsa_BirthdayS1");
401 ftab_MakeTables (fam, res, cho, Par, TabBirthdayS, Nr, j1, j2, jstep);
402 ftab_PrintTable2 (res->Exp, res->Obs, FALSE);
403 ftab_PrintTable (res->PVal2);
404 if (localRes)
405 fres_DeletePoisson (res);
406 }
407
408
409 /*========================================================================*/
410
fmarsa_SerialOver1(void)411 void fmarsa_SerialOver1 (void)
412 {
413 util_Error ("fmarsa_SerialOver1: use fmultin_SerialOver1 instead");
414 }
415
416 /*========================================================================*/
417
fmarsa_CollisionOver1(void)418 void fmarsa_CollisionOver1 (void)
419 {
420 util_Error ("fmarsa_CollisionOver1: use fmultin_SerialOver1 instead");
421 }
422
423
424 /*=========================================================================*/
425
TabGCD(ffam_Fam * fam,void * res1,void * cho,void * par1,int i,int j,int irow,int icol)426 static void TabGCD (ffam_Fam * fam, void *res1, void *cho,
427 void *par1, int i, int j, int irow, int icol)
428 {
429 int r, s;
430 long N, n;
431 const long *Par = par1;
432 fmarsa_Res2 *fres = res1;
433 smarsa_Res2 *sres;
434
435 N = Par[0];
436 r = Par[1];
437 s = Par[2];
438
439 n = fcho_ChooseParamL (cho, (long) (3.0 * gofs_MinExpected),
440 fmarsa_Maxn, i, j);
441 if (n <= 0)
442 return;
443 s = fcho_Chooses (r, s, fam->Resol[irow]);
444 if (s <= 0)
445 return;
446
447 sres = smarsa_CreateRes2 ();
448 smarsa_GCD (fam->Gen[irow], sres, N, n, r, s);
449 fres_FillTableEntryC (fres->GCD, sres->GCD->pVal2, N, irow, icol);
450 fres_FillTableEntryC (fres->NumIter, sres->NumIter->pVal2, N, irow, icol);
451 smarsa_DeleteRes2 (sres);
452 }
453
454
455 /*-------------------------------------------------------------------------*/
456
fmarsa_GCD1(ffam_Fam * fam,fmarsa_Res2 * res,fcho_Cho * cho,long N,int r,int s,int Nr,int j1,int j2,int jstep)457 void fmarsa_GCD1 (ffam_Fam *fam, fmarsa_Res2 *res, fcho_Cho *cho,
458 long N, int r, int s, int Nr, int j1, int j2, int jstep)
459 {
460 long Par[3] = { 0 };
461 lebool localRes;
462
463 Par[0] = N;
464 Par[1] = r;
465 Par[2] = s;
466 if (res == NULL) {
467 localRes = TRUE;
468 res = fmarsa_CreateRes2 ();
469 } else
470 localRes = FALSE;
471
472 PrintHead ("fmarsa_GCD1", fam, N, 0, r, s, 0, 0, 0, Nr, j1, j2, jstep);
473 InitRes2 (fam, res, N, Nr, j1, j2, jstep,
474 "fmarsa_GCD1, GCD", "fmarsa_GCD1, NumIter");
475 ftab_MakeTables (fam, res, cho, Par, TabGCD, Nr, j1, j2, jstep);
476 fres_PrintCont (res->GCD);
477 /* fres_PrintCont (res->NumIter); */
478 if (localRes)
479 fmarsa_DeleteRes2 (res);
480 }
481
482
483 /*=========================================================================*/
484