1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: smarsa.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 #include "gdef.h"
32 #include "util.h"
33 #include "tables.h"
34 #include "chrono.h"
35 #include "num.h"
36
37 #include "smarsa.h"
38 #include "smultin.h"
39 #include "wdist.h"
40 #include "swrite.h"
41 #include "unif01.h"
42
43 #include "vectorsF2.h"
44
45 #include "gofs.h"
46 #include "gofw.h"
47 #include "fdist.h"
48 #include "fbar.h"
49 #include "fmass.h"
50 #include "statcoll.h"
51
52 #include <math.h>
53 #include <stdio.h>
54 #include <stdlib.h>
55
56
57 #define LENGTH 200
58
59 /* MAXK = 2^64 */
60 #define STR_MAXK "18446744073709551616"
61
62
63
64 /*---------------------------- Extern variables ---------------------------*/
65
66 #ifdef USE_LONGLONG
67 double smarsa_Maxk = 18446744073709551616.0; /* 2^64 */
68 #else
69 double smarsa_Maxk = num_MaxIntDouble; /* 2^53 */
70 #endif
71
72
73
74
75 /*------------------------------- Functions -------------------------------*/
76
WriteResultsPoisson(sres_Poisson * res,long N)77 static void WriteResultsPoisson (sres_Poisson *res, long N)
78 {
79 printf ("\n----------------------------------------------------"
80 "\nTotal expected number = N*Lambda : ");
81 num_WriteD (N * res->Lambda, 10, 2, 2);
82 printf ("\nTotal observed number : %7ld\n",
83 (long) res->sVal2);
84 gofw_Writep1 (res->pVal2);
85 printf ("\n");
86 }
87
88
89 /*=========================================================================*/
90
InitRes(smarsa_Res * res,long N,double Lambda,char * nam)91 static void InitRes (
92 smarsa_Res *res, /* Results holder */
93 long N, /* Number of replications */
94 double Lambda, /* Poisson mean */
95 char *nam /* Test name */
96 )
97 /*
98 * Initializes res
99 */
100 {
101 sres_InitBasic (res->Bas, N, nam);
102 sres_InitPoisson (res->Pois, N, Lambda, nam);
103 }
104
105
106 /*-------------------------------------------------------------------------*/
107
smarsa_CreateRes(void)108 smarsa_Res * smarsa_CreateRes (void)
109 {
110 smarsa_Res *res;
111 res = util_Malloc (sizeof (smarsa_Res));
112 res->Bas = sres_CreateBasic ();
113 res->Pois = sres_CreatePoisson ();
114 res->Pois->pLeft = -1.0;
115 res->Pois->pRight = -1.0;
116 return res;
117 }
118
119
120 /*-------------------------------------------------------------------------*/
121
smarsa_DeleteRes(smarsa_Res * res)122 void smarsa_DeleteRes (smarsa_Res *res)
123 {
124 if (res == NULL)
125 return;
126 sres_DeleteBasic (res->Bas);
127 sres_DeletePoisson (res->Pois);
128 util_Free (res);
129 }
130
131
132 /*=========================================================================*/
133
InitRes2(smarsa_Res2 * res,long N,int jmax,int tmax)134 static void InitRes2 (
135 smarsa_Res2 *res, /* Results holder */
136 long N, /* Number of replications */
137 int jmax, /* Max class index for GCD */
138 int tmax /* Max class index for NumIter */
139 )
140 /*
141 * Initializes the smarsa_Res2 structure
142 */
143 {
144 sres_InitChi2 (res->GCD, N, jmax, "smarsa_GCD: GCD");
145 sres_InitChi2 (res->NumIter, N, tmax, "smarsa_GCD: NumIter");
146 }
147
148
149 /*-------------------------------------------------------------------------*/
150
smarsa_CreateRes2(void)151 smarsa_Res2 *smarsa_CreateRes2 (void)
152 {
153 smarsa_Res2 *res;
154 res = util_Malloc (sizeof (smarsa_Res2));
155 res->GCD = sres_CreateChi2 ();
156 res->NumIter = sres_CreateChi2 ();
157 return res;
158 }
159
160
161 /*-------------------------------------------------------------------------*/
162
smarsa_DeleteRes2(smarsa_Res2 * res)163 void smarsa_DeleteRes2 (smarsa_Res2 *res)
164 {
165 if (res == NULL)
166 return;
167 sres_DeleteChi2 (res->GCD);
168 sres_DeleteChi2 (res->NumIter);
169 util_Free (res);
170 }
171
172
173 /*=========================================================================*/
174
smarsa_SerialOver(unif01_Gen * gen,sres_Basic * res,long N,long n,int r,long d,int t)175 void smarsa_SerialOver (unif01_Gen *gen, sres_Basic *res,
176 long N, long n, int r, long d, int t)
177 {
178 double ValDelta[] = { 1.0 };
179 smultin_Param *par;
180
181 if (swrite_Basic)
182 printf ("***********************************************************\n"
183 "Test smarsa_SerialOver calling smultin_MultinomialOver\n\n");
184
185 par = smultin_CreateParam (1, ValDelta, smultin_GenerCellSerial, 0);
186 if (NULL == res) {
187 smultin_MultinomialOver (gen, par, NULL, N, n, r, d, t, FALSE);
188 } else {
189 smultin_Res *resm;
190 resm = smultin_CreateRes (par);
191 smultin_MultinomialOver (gen, par, resm, N, n, r, d, t, FALSE);
192 sres_InitBasic (res, N, "smarsa_SerialOver");
193 statcoll_SetDesc (res->sVal1, "SerialOver sVal1");
194 res->sVal1->NObs = resm->Collector[0]->NObs;
195 tables_CopyTabD (resm->Collector[0]->V, res->sVal1->V, 1, N);
196 tables_CopyTabD (resm->sVal2[0], res->sVal2, 0, gofw_NTestTypes - 1);
197 tables_CopyTabD (resm->pVal2[0], res->pVal2, 0, gofw_NTestTypes - 1);
198 smultin_DeleteRes (resm);
199 }
200 smultin_DeleteParam (par);
201 }
202
203
204 /*=========================================================================*/
205
smarsa_CollisionOver(unif01_Gen * gen,smarsa_Res * res,long N,long n,int r,long d,int t)206 void smarsa_CollisionOver (unif01_Gen *gen, smarsa_Res *res,
207 long N, long n, int r, long d, int t)
208 {
209 double ValDelta[] = { -1.0 };
210 smultin_Param *par;
211
212 if (swrite_Basic)
213 printf ("***********************************************************\n"
214 "Test smarsa_CollisionOver calling smultin_MultinomialOver\n\n");
215
216 par = smultin_CreateParam (1, ValDelta, smultin_GenerCellSerial, 3);
217 if (NULL == res) {
218 smultin_MultinomialOver (gen, par, NULL, N, n, r, d, t, TRUE);
219 } else {
220 smultin_Res *resm;
221 resm = smultin_CreateRes (par);
222 smultin_MultinomialOver (gen, par, resm, N, n, r, d, t, TRUE);
223 InitRes (res, N, resm->Mu[0], "smarsa_CollisionOver");
224 statcoll_SetDesc (res->Bas->sVal1, "CollisionOver sVal1");
225 statcoll_SetDesc (res->Pois->sVal1, "CollisionOver sVal1");
226 res->Pois->sVal1->NObs = resm->Collector[0]->NObs;
227 res->Bas->sVal1->NObs = resm->Collector[0]->NObs;
228 tables_CopyTabD (resm->Collector[0]->V, res->Bas->sVal1->V, 1, N);
229 tables_CopyTabD (resm->Collector[0]->V, res->Pois->sVal1->V, 1, N);
230 res->Pois->pVal2 = resm->pColl;
231 if (resm->CollApprox == smultin_CollPoissonSparse)
232 res->Pois->sVal2 = resm->NbCollisions;
233 else
234 res->Pois->sVal2 = resm->NbCells[0];
235 tables_CopyTabD (resm->sVal2[0], res->Bas->sVal2, 0,
236 gofw_NTestTypes - 1);
237 tables_CopyTabD (resm->pVal2[0], res->Bas->pVal2, 0,
238 gofw_NTestTypes - 1);
239 smultin_DeleteRes (resm);
240 }
241 smultin_DeleteParam (par);
242 }
243
244
245 /*=========================================================================*/
246
smarsa_Opso(unif01_Gen * gen,smarsa_Res * res,long N,int r,int p)247 void smarsa_Opso (unif01_Gen * gen, smarsa_Res * res, long N, int r, int p)
248 {
249 int d;
250 long NBalls;
251
252 switch (p) {
253 case 1:
254 NBalls = 2097152;
255 d = 1024;
256 break;
257 case 2:
258 NBalls = 4194304;
259 d = 2048;
260 break;
261 case 3:
262 NBalls = 8388608;
263 d = 2048;
264 break;
265 default:
266 util_Error ("smarsa_Opso: p must be in {1, 2, 3}");
267 }
268
269 if (swrite_Basic)
270 printf ("***********************************************************\n"
271 "Test smarsa_Opso calling smarsa_CollisionOver\n\n");
272 smarsa_CollisionOver (gen, res, N, NBalls, r, d, 2);
273 }
274
275
276 /*=========================================================================*/
277 /*
278 * The CPU time needed for BirthdaySpacings is 6 times longer when I used
279 * the standard function qsort of stdlib.h. Thus we use our own QuickSort.
280 */
281
282 #undef QSORT
283 #ifdef QSORT
compareD(const void * p0,const void * q0)284 static int compareD (const void *p0, const void *q0)
285 {
286 double x = *((const double *) p0);
287 double y = *((const double *) q0);
288 return (x < y) ? -1 : (x > y) ? 1 : 0;
289 }
290 /* qsort ((void *)(DatDiff + 1), (size_t) n, sizeof (double), compareD); */
291 #endif
292
293
294 /*=========================================================================*/
295
WriteDataBirth(unif01_Gen * gen,char * TestName,long N,long n,int r,long d,int t,int p,double k,smultin_CellType kc,double Lambda)296 static void WriteDataBirth (unif01_Gen * gen, char *TestName, long N, long n,
297 int r, long d, int t, int p, double k, smultin_CellType kc,
298 double Lambda)
299 {
300 swrite_Head (gen, TestName, N, n, r);
301 printf (", d = %1ld, t = %1d, p = %1d\n\n", d, t, p);
302 #ifdef USE_LONGLONG
303 if (kc == 0 && d > 1) /* kc = 2^64 */
304 printf ("\n Number of cells = d^t = " STR_MAXK "\n");
305 else
306 printf ("\n Number of cells = d^t = %18" PRIuLEAST64 "\n", kc);
307 #else
308 printf ("\n Number of cells = d^t = %16.0f\n", k);
309 #endif
310 printf (" Lambda = Poisson mean = ");
311 num_WriteD (Lambda, 12, 4, 2);
312 printf ("\n\n");
313 }
314
315
316 /*-------------------------------------------------------------------------*/
317
smarsa_BirthdaySpacings(unif01_Gen * gen,sres_Poisson * res,long N,long n,int r,long d,int t,int Order)318 void smarsa_BirthdaySpacings (unif01_Gen *gen, sres_Poisson *res,
319 long N, long n, int r, long d, int t, int Order)
320 {
321 long Seq; /* Replication number */
322 long j;
323 long Sum;
324 double Y; /* Number of collisions */
325 double k;
326 smultin_CellType kc;
327 double Lambda; /* Poisson mean */
328 smultin_CellType *Dates, *DatDiff;
329 fmass_INFO Mass;
330 char str[LENGTH + 1];
331 lebool localRes = FALSE;
332 chrono_Chrono *Timer;
333 char *TestName = "smarsa_BirthdaySpacings test";
334
335 Timer = chrono_Create ();
336 kc = k = d;
337 for (j = 2; j <= t; j++) {
338 k *= d;
339 kc *= d;
340 }
341 Lambda = (double) n * n / k * (n / 4.0);
342
343 if (swrite_Basic)
344 WriteDataBirth (gen, TestName, N, n, r, d, t, Order, k, kc, Lambda);
345
346 if (d <= 1) {
347 util_Warning (TRUE,
348 "smarsa_BirthdaySpacings: d <= 1. The test is not done.");
349 return;
350 }
351 if (k > smarsa_Maxk) {
352 util_Warning (TRUE,
353 "smarsa_BirthdaySpacings: d^t > smarsa_Maxk. The test is not done.");
354 return;
355 }
356 if (8.0 * N * Lambda > sqrt (sqrt (k))) {
357 util_Warning (TRUE,
358 "smarsa_BirthdaySpacings: 8N Lambda > k^(1/4). The test is not done.");
359 return;
360 }
361 if (res == NULL) {
362 localRes = TRUE;
363 res = sres_CreatePoisson ();
364 }
365 sres_InitPoisson (res, N, Lambda, "smarsa_BirthdaySpacings");
366
367 Dates = util_Calloc (1 + (size_t) n, sizeof (smultin_CellType));
368 DatDiff = util_Calloc (1 + (size_t) n, sizeof (smultin_CellType));
369
370 sprintf (str, "The N statistic values (a Poisson with mean %g):", Lambda);
371 statcoll_SetDesc (res->sVal1, str);
372
373 Sum = 0;
374 for (Seq = 1; Seq <= N; Seq++) {
375 /* Generate and sort the "birth dates" */
376 if (Order == 2) {
377 for (j = 1; j <= n; j++) {
378 Dates[j] = smultin_GenerCellSerial2 (gen, r, t, d);
379 }
380 } else {
381 for (j = 1; j <= n; j++) {
382 Dates[j] = smultin_GenerCellSerial (gen, r, t, d);
383 }
384 }
385 #ifdef USE_LONGLONG
386 tables_QuickSortULL (Dates, 1, n);
387 /* Compute the differences between adjacent dates */
388 gofs_DiffULL (Dates, DatDiff, 1, n, 0ULL, 1ULL);
389 /* The last cell is a special case */
390 DatDiff[n] = kc - Dates[n] + Dates[1];
391 tables_QuickSortULL (DatDiff, 1, n);
392 #else
393 tables_QuickSortD (Dates, 1, n);
394 /* Compute the differences between adjacent dates */
395 gofs_DiffD (Dates, DatDiff, 1, n, 0.0, 1.0);
396 /* The last cell is a special case */
397 DatDiff[n] = kc - Dates[n] + Dates[1];
398 tables_QuickSortD (DatDiff, 1, n);
399 #endif
400
401 /* Count the number of collisions in DatDiff */
402 Y = 0.0;
403 for (j = 2; j <= n; j++) {
404 if (DatDiff[j] == DatDiff[j - 1])
405 Y += 1.0;
406 }
407 Sum += Y;
408 statcoll_AddObs (res->sVal1, Y);
409 if (swrite_Counters) {
410 #ifdef USE_LONGLONG
411 tables_WriteTabULL (Dates, 1, n, 3, 21, "Birthdates:");
412 tables_WriteTabULL (DatDiff, 1, n, 3, 21, "Birthdate differences:");
413 #else
414 tables_WriteTabD (Dates, 1, n, 4, 17, 0, 0, "Birthdates:");
415 tables_WriteTabD (DatDiff, 1, n, 4, 17, 0, 0,
416 "Birthdate differences:");
417 #endif
418 }
419 }
420
421 res->sVal2 = Sum;
422 Mass = fmass_CreatePoisson (N * Lambda);
423 res->pLeft = fdist_Poisson2 (Mass, Sum);
424 res->pRight = fbar_Poisson2 (Mass, Sum);
425 fmass_DeletePoisson (Mass);
426 res->pVal2 = gofw_pDisc (res->pLeft, res->pRight);
427
428 if (swrite_Collectors)
429 statcoll_Write (res->sVal1, 5, 14, 1, 1);
430 if (swrite_Basic) {
431 WriteResultsPoisson (res, N);
432 swrite_Final (gen, Timer);
433 }
434 util_Free (Dates);
435 util_Free (DatDiff);
436 if (localRes)
437 sres_DeletePoisson (res);
438 chrono_Delete (Timer);
439 }
440
441
442 /*=========================================================================*/
443
WriteDataCAT(unif01_Gen * gen,char * TestName,long N,long n,int r,long d,int t,long S[],double Lambda)444 static void WriteDataCAT (unif01_Gen *gen, char *TestName,
445 long N, long n, int r, long d, int t, long S[], double Lambda)
446 {
447 int i;
448 swrite_Head (gen, TestName, N, n, r);
449 printf (", d = %1ld, t = %1d\n\n", d, t);
450 for (i = 0; i < t; i++) {
451 printf (" S[%1d] = %1ld\n", i, S[i]);
452 }
453 printf ("\n Lambda = Poisson mean = ");
454 num_WriteD (Lambda, 12, 4, 2);
455 printf ("\n\n");
456 }
457
458
459 /*-------------------------------------------------------------------------*/
460
TestCATData(long d,int t,long S1[])461 static void TestCATData (long d, int t, long S1[])
462 /*
463 * Test that the key to search for has no overlap, that is cannot be
464 * written as ABA, where A and B are parts of the key.
465 */
466 {
467 int i, j, s;
468 long k1, k2;
469 i = 0;
470 j = t - 1;
471 k1 = k2 = 0;
472 while (i < j) {
473 k1 = k1 * d + S1[i];
474 k2 = 0;
475 for (s = j; s < t; s++)
476 k2 = k2 * d + S1[s];
477 util_Assert (k1 != k2,
478 "CATData: target cell number of the form ABA");
479 i++;
480 j--;
481 }
482 }
483
484
485 /*-------------------------------------------------------------------------*/
486 #if 0
487 static void CATGenere1 (
488 unif01_Gen *gen,
489 long n, /* Number of points */
490 int r, /* Drop the first r bits of each U01 */
491 long d, /* Number of segments on the 1-dim. line */
492 int t, /* Dimension */
493 long Key, /* Key to search for */
494 long k1, /* = d^(t-1) */
495 long *Count /* Number of times Key appears */
496 )
497 /*
498 * Generate the n points in the dense case and count the number of times
499 * cell Key appears. This is the circular version with n points. It also
500 * correspond to the case of aperiodic Key.
501 */
502 {
503 int j, i;
504 long Indice = 0;
505 long Y = 0; /* Counter */
506 long Premier[32];
507
508 util_Assert (t <= 32, "smarsa_CAT.Genere: t > 32");
509
510 /* Generation of the first (t - 1) elements of the first tuple */
511 for (j = 1; j < t; j++) {
512 Premier[j] = unif01_StripL (gen, r, d);
513 Indice = Indice * d + Premier[j];
514 }
515
516 /* Generation of the n - (tt - 1) tuples */
517 for (j = 1; j <= n - (t - 1); j++) {
518 /* Remove the leftmost component ... */
519 Indice %= k1;
520 /* ... shift and get another for the rightmost one */
521 Indice = Indice * d + unif01_StripL (gen, r, d);
522 if (Indice == Key) {
523 ++Y;
524 /* Key found: jump over the whole Indice and restart */
525 Indice = 0;
526 for (i = 1; i < t; i++) {
527 Indice = Indice * d + unif01_StripL (gen, r, d);
528 j++;
529 }
530 }
531 }
532
533 /* Generation of the last (t - 1) tuples. We use numbers in array */
534 /* Premier[] so that the sequence is in fact circular */
535 for (j = 1; j < t; j++) {
536 Indice %= k1;
537 Indice = Indice * d + Premier[j];
538 if (Indice == Key)
539 ++Y;
540 }
541
542 *Count = Y;
543 }
544 #endif
545
546 /*-------------------------------------------------------------------------*/
547
CATGenere(unif01_Gen * gen,long n,int r,long d,int t,long Key,long k1,long * Count)548 static void CATGenere (
549 unif01_Gen *gen,
550 long n, /* Number of points */
551 int r, /* Drop the first r bits of each U01 */
552 long d, /* Number of segments on the 1-dim. line */
553 int t, /* Dimension */
554 long Key, /* Key to search for */
555 long k1, /* = d^(t-1) */
556 long *Count /* Number of times Key appears */
557 )
558 /*
559 * Generate the n points in the dense case and count the number of times
560 * cell Key appears. This is the non-circular version with n - t + 1 points.
561 * It also correspond to the case of aperiodic Key.
562 */
563 {
564 int j, i;
565 long Indice;
566 long Y = 0; /* Counter */
567
568 /* Generation of the first (t - 1) elements of the first tuple */
569 Indice = 0;
570 for (j = 1; j < t; j++)
571 Indice = Indice * d + unif01_StripL (gen, r, d);
572
573 /* Generation of the n - (tt - 1) tuples */
574 for (j = 1; j <= n - (t - 1); j++) {
575 /* Remove the leftmost component ... */
576 Indice %= k1;
577 /* ... shift and get another for the rightmost one */
578 Indice = Indice * d + unif01_StripL (gen, r, d);
579 if (Indice == Key) {
580 ++Y;
581 /* Key found: jump over the whole Indice and restart */
582 Indice = 0;
583 for (i = 1; i < t; i++) {
584 Indice = Indice * d + unif01_StripL (gen, r, d);
585 j++;
586 }
587 }
588 }
589
590 *Count = Y;
591 }
592
593
594 /*-------------------------------------------------------------------------*/
595
smarsa_CAT(unif01_Gen * gen,sres_Poisson * res,long N,long n,int r,long d,int t,long S[])596 void smarsa_CAT (unif01_Gen *gen, sres_Poisson *res,
597 long N, long n, int r, long d, int t, long S[])
598 {
599 long Seq;
600 long i;
601 double k;
602 long k1; /* d^(t-1) */
603 long Key; /* Cell number to search for */
604 double Lambda; /* Poisson mean */
605 long Sum;
606 long Co;
607 fmass_INFO Mass;
608 char str[LENGTH + 1];
609 lebool localRes = FALSE;
610 chrono_Chrono *Timer;
611 char *TestName = "smarsa_CAT test";
612
613 Timer = chrono_Create ();
614 k1 = d;
615 for (i = 2; i < t; i++)
616 k1 *= d;
617 k = k1 * d;
618 Lambda = (n - t + 1) / k;
619 if (swrite_Basic)
620 WriteDataCAT (gen, TestName, N, n, r, d, t, S, Lambda);
621 util_Assert (d > 1, "smarsa_CAT: d <= 1");
622
623 Key = 0;
624 for (i = 0; i < t; i++) {
625 if (S[i] < 0 || S[i] >= d) {
626 util_Error ("smarsa_CAT: S[i] must be in [0, d - 1]");
627 }
628 Key = Key * d + S[i];
629 }
630 TestCATData (d, t, S);
631 if (res == NULL) {
632 localRes = TRUE;
633 res = sres_CreatePoisson ();
634 }
635 sres_InitPoisson (res, N, Lambda, "smarsa_CAT");
636 sprintf (str, "The N statistic values (a Poisson with mean %g):", Lambda);
637 statcoll_SetDesc (res->sVal1, str);
638
639 Sum = 0;
640 for (Seq = 1; Seq <= N; Seq++) {
641 CATGenere (gen, n, r, d, t, Key, k1, &Co);
642 statcoll_AddObs (res->sVal1, (double) Co);
643 Sum += Co;
644 }
645
646 res->sVal2 = Sum;
647 Mass = fmass_CreatePoisson (res->Mu);
648 res->pLeft = fdist_Poisson2 (Mass, Sum);
649 res->pRight = fbar_Poisson2 (Mass, Sum);
650 fmass_DeletePoisson (Mass);
651 res->pVal2 = gofw_pDisc (res->pLeft, res->pRight);
652
653 if (swrite_Collectors)
654 statcoll_Write (res->sVal1, 5, 14, 1, 1);
655 if (swrite_Basic) {
656 WriteResultsPoisson (res, N);
657 swrite_Final (gen, Timer);
658 }
659 if (localRes)
660 sres_DeletePoisson (res);
661 chrono_Delete (Timer);
662 }
663
664
665 /*=========================================================================*/
666
WriteDataCATBits(unif01_Gen * gen,char * TestName,long N,long n,int r,int s,int L,unsigned long Key,double Lambda)667 static void WriteDataCATBits (unif01_Gen *gen, char *TestName,
668 long N, long n, int r, int s, int L, unsigned long Key, double Lambda)
669 {
670 swrite_Head (gen, TestName, N, n, r);
671 printf (", s = %1d, L = %1d, Key = %lu\n\n", s, L, Key);
672 printf (" Lambda = Poisson mean = ");
673 num_WriteD (Lambda, 12, 4, 2);
674 printf ("\n\n");
675 }
676
677
678 /*-------------------------------------------------------------------------*/
679
TestCATBitsData(int L,unsigned long Key)680 static void TestCATBitsData (int L, unsigned long Key)
681 /*
682 * Test that the key to search for has no overlap, that is cannot be
683 * written as ABA, where A and B are parts of the key.
684 */
685 {
686 int i;
687 unsigned long mask = 1, shift = L - 1;
688 i = 0;
689 while (i < L / 2) {
690 if ((Key & mask) == (Key >> shift)) {
691 bitset_WriteSet ("Key = ", Key, L);
692 util_Error ("CATBitsData: Key of the form ABA");
693 }
694 i++;
695 shift--;
696 mask = num_TwoExp[i + 1] - 1.0;
697 }
698 }
699
700
701 /*-------------------------------------------------------------------------*/
702
CATGenerBits(unif01_Gen * gen,long n,int r,int s,int L,unsigned long KEY0,long * Count)703 static void CATGenerBits (unif01_Gen *gen, long n, int r, int s,
704 int L, unsigned long KEY0, long *Count)
705 {
706 /*
707 * Generate the bits in the CATBits test. Points are generated with
708 * overlapping. We have a window of size L bits, and we slide it 1 bit
709 * forward at each step to generate a point. We then check whether it
710 * equals the L bits Key.
711 */
712 const unsigned long MASK0 = num_TwoExp[L] - 1.0;
713 unsigned long Mask, Key, Z0, Z;
714 int j0, j, k;
715 long i;
716 long co;
717
718 util_Assert (L <= 32, "CATBits: GenerBits: L > 32");
719 co = 0;
720
721 if ((s >= L) && (L <= 16)) {
722 const int q = s - L;
723
724 /* Make sure to skip the first half of the loop for the first number
725 since there is no previous Z */
726 j0 = L;
727 Z = 0;
728
729 for (i = 0; i < n / s; i++) {
730 Z0 = unif01_StripB (gen, r, s);
731
732 /* The last L - j0 bits of the previous number */
733 Mask = MASK0 << (L - j0);
734 Key = KEY0 << (L - j0);
735 Z |= (Z0 >> q);
736 j = j0;
737 while (j < L) {
738 if (Key == (Z & Mask)) {
739 co++;
740 j += L;
741 Mask >>= L;
742 Key >>= L;
743 } else {
744 j++;
745 Mask >>= 1;
746 Key >>= 1;
747 }
748 }
749 j0 = j % L;
750
751 /* The first s - L bits of the current number */
752 Z = Z0;
753 Mask = MASK0 << (q - j0);
754 Key = KEY0 << (q - j0);
755 j = j0;
756 while (j < q) {
757 if (Key == (Z & Mask)) {
758 co++;
759 j += L;
760 Mask >>= L;
761 Key >>= L;
762 } else {
763 j++;
764 Mask >>= 1;
765 Key >>= 1;
766 }
767 }
768 j0 = j - q;
769 Z = Z0 << L;
770 }
771
772 } else if (s >= L) {
773 #ifdef USE_LONGLONG
774 const ulonglong MASK0 = num_TwoExp[L] - 1.0;
775 ulonglong Z, Z0;
776 ulonglong Mask, Key;
777 const int q = s - L;
778
779 /* Make sure to skip the first half of the loop for the first number
780 since there is no previous Z */
781 j0 = L;
782 Z = 0;
783
784 for (i = 0; i < n / s; i++) {
785 Z0 = unif01_StripB (gen, r, s);
786
787 /* The last L - j0 bits of the previous number */
788 Mask = MASK0 << (L - j0);
789 Key = KEY0 << (L - j0);
790 Z |= (Z0 >> q);
791 j = j0;
792 while (j < L) {
793 if (Key == (Z & Mask)) {
794 co++;
795 j += L;
796 Mask >>= L;
797 Key >>= L;
798 } else {
799 j++;
800 Mask >>= 1;
801 Key >>= 1;
802 }
803 }
804 j0 = j % L;
805
806 /* The first s - L bits of the current number */
807 Z = Z0;
808 Mask = MASK0 << (q - j0);
809 Key = KEY0 << (q - j0);
810 j = j0;
811 while (j < q) {
812 if (Key == (Z & Mask)) {
813 co++;
814 j += L;
815 Mask >>= L;
816 Key >>= L;
817 } else {
818 j++;
819 Mask >>= 1;
820 Key >>= 1;
821 }
822 }
823 j0 = j - q;
824 Z = Z0 << L;
825 }
826 #else
827 if (L <= s)
828 util_Error ("CATGenerBits: L <= s and L > 16");
829 #endif
830
831 } else if ((s < L) && (L + s <= 32)) {
832 const int t = L / s;
833 util_Assert (L % s == 0, "CATBits: L > s but L % s not 0");
834
835 /* Generation of the first L random bits */
836 Z = 0;
837 for (j = 0; j < t; j++) {
838 Z <<= s;
839 Z |= unif01_StripB (gen, r, s);
840 }
841 j0 = 0;
842
843 /* Generation of the rest of the random bits */
844 for (i = 0; i < (n - L) / s; i++) {
845 Z = (Z << s) | unif01_StripB (gen, r, s);
846 Mask = MASK0 << (s - j0);
847 Key = KEY0 << (s - j0);
848 j = j0;
849 while (j < s) {
850 if (Key == (Z & Mask)) {
851 co++;
852 j += L;
853 i += t - 1;
854 for (k = 1; k < t; k++) {
855 Z <<= s;
856 Z |= unif01_StripB (gen, r, s);
857 }
858 } else {
859 j++;
860 Mask >>= 1;
861 Key >>= 1;
862 }
863 }
864 j0 = j % s;
865 }
866
867 } else {
868 #ifdef USE_LONGLONG
869 const ulonglong MASK0 = num_TwoExp[L] - 1.0;
870 const int t = L / s;
871 ulonglong Z;
872 ulonglong Mask, Key, Key0 = KEY0;
873
874 if (L > s) {
875 util_Assert (L % s == 0, "CATBits: L > s but L % s not 0");
876 }
877
878 /* Generation of the first L random bits */
879 Z = 0;
880 for (j = 0; j < t; j++) {
881 Z <<= s;
882 Z |= unif01_StripB (gen, r, s);
883 }
884 j0 = 0;
885
886 /* Generation of the rest of the random bits */
887 for (i = 0; i < (n - L) / s; i++) {
888 Z = (Z << s) | unif01_StripB (gen, r, s);
889 Mask = MASK0 << (s - j0);
890 Key = Key0 << (s - j0);
891 j = j0;
892 while (j < s) {
893 if (Key == (Z & Mask)) {
894 co++;
895 j += L;
896 i += t - 1;
897 for (k = 1; k < t; k++) {
898 Z <<= s;
899 Z |= unif01_StripB (gen, r, s);
900 }
901 } else {
902 j++;
903 Mask >>= 1;
904 Key >>= 1;
905 }
906 }
907 j0 = j % s;
908 }
909 #else
910 if (L == s)
911 util_Error ("CATGenereBits: L = s and s > 16");
912 else
913 util_Error ("CATGenereBits: L > s and L + s > 32");
914 #endif
915 }
916
917 *Count = co;
918 }
919
920
921 /*-------------------------------------------------------------------------*/
922
smarsa_CATBits(unif01_Gen * gen,sres_Poisson * res,long N,long n,int r,int s,int L,unsigned long Key)923 void smarsa_CATBits (unif01_Gen *gen, sres_Poisson *res,
924 long N, long n, int r, int s, int L, unsigned long Key)
925 {
926 long Seq;
927 double Lambda; /* Poisson mean */
928 long Sum;
929 long Co;
930 fmass_INFO Mass;
931 char str[LENGTH + 1];
932 lebool localRes = FALSE;
933 chrono_Chrono *Timer;
934 char *TestName = "smarsa_CATBits test";
935
936 Timer = chrono_Create ();
937 Lambda = (n - L + 1) / num_TwoExp[L];
938 if (swrite_Basic)
939 WriteDataCATBits (gen, TestName, N, n, r, s, L, Key, Lambda);
940 util_Assert (L > 1, "smarsa_CATBits: L <= 1");
941
942 TestCATBitsData (L, Key);
943 if (res == NULL) {
944 localRes = TRUE;
945 res = sres_CreatePoisson ();
946 }
947 sres_InitPoisson (res, N, Lambda, "smarsa_CATBits");
948 sprintf (str, "The N statistic values (a Poisson with mean %g):", Lambda);
949 statcoll_SetDesc (res->sVal1, str);
950
951 Sum = 0;
952 for (Seq = 1; Seq <= N; Seq++) {
953 CATGenerBits (gen, n, r, s, L, Key, &Co);
954 statcoll_AddObs (res->sVal1, (double) Co);
955 Sum += Co;
956 }
957
958 res->sVal2 = Sum;
959 Mass = fmass_CreatePoisson (res->Mu);
960 res->pLeft = fdist_Poisson2 (Mass, Sum);
961 res->pRight = fbar_Poisson2 (Mass, Sum);
962 fmass_DeletePoisson (Mass);
963 res->pVal2 = gofw_pDisc (res->pLeft, res->pRight);
964
965 if (swrite_Collectors)
966 statcoll_Write (res->sVal1, 5, 14, 1, 1);
967 if (swrite_Basic) {
968 WriteResultsPoisson (res, N);
969 swrite_Final (gen, Timer);
970 }
971 if (localRes)
972 sres_DeletePoisson (res);
973 chrono_Delete (Timer);
974 }
975
976
977 /*=========================================================================*/
978
WriteDataMatRank(unif01_Gen * gen,char * TestName,long N,long n,int r,int s,int L,int k)979 static void WriteDataMatRank (unif01_Gen * gen, char *TestName,
980 long N, long n, int r, int s, int L, int k)
981 {
982 swrite_Head (gen, TestName, N, n, r);
983 printf (", s = %1d, L = %1d, k = %1d\n\n", s, L, k);
984 }
985
986
987 /*-------------------------------------------------------------------------*/
988 #if 0
989
990 static int RankOfBitMatrix (bitset_BitSet M[], int maxrow)
991 /*
992 * Calculation of the rank of the bit-matrix M
993 */
994 {
995 const int MaxBit = 31; /* number of bits in a word - 1 */
996 bitset_BitSet Swap;
997 int rank = 0;
998 int i;
999 int CL = 1;
1000
1001 while (CL <= MaxBit) {
1002 /* All components of M shift their bits 1 position to the left */
1003 for (i = 0; i < maxrow; i++)
1004 M[i] <<= 1;
1005
1006 /* Search of the first M[i] with 1 as the major bit */
1007 i = rank;
1008 for (;;) {
1009 if ((bitset_TestBit (M[i], MaxBit)) || (i == maxrow - 1))
1010 break;
1011 ++i;
1012 }
1013 /* Diagonalization of matrix M */
1014 if (i < maxrow - 1) {
1015 Swap = M[rank];
1016 M[rank] = M[i];
1017 M[i] = Swap;
1018 for (i = rank + 1; i < maxrow; i++) {
1019 if (bitset_TestBit (M[i], MaxBit))
1020 M[i] ^= M[rank];
1021 }
1022 ++rank;
1023 if (rank == MaxBit)
1024 return rank;
1025 }
1026 ++CL;
1027 }
1028 return rank;
1029 }
1030
1031
1032 /*-------------------------------------------------------------------------*/
1033
1034 #define lmax 64
1035
1036 void smarsa_MatrixRank (unif01_Gen *gen, sres_Chi2 *res,
1037 long N, long n, int r, int s, int l, int k)
1038 {
1039 long Seq;
1040 long Rep;
1041 int j;
1042 int i;
1043 long L; /* One line of bits */
1044 int c; /* Number-1 of U01 used to build a line */
1045 long d; /* Get s bits of a generated U01 */
1046 long a; /* Get b bits of a generated U01 */
1047 int b; /* Number of bits of last U01 of a line */
1048 int Minkl; /* Min (k, l) */
1049 long NbGroups; /* Number of classes for ChiSquare */
1050 long jhigh; /* Index of the highest class */
1051 long jlow; /* Index of the lowest class */
1052 int Rank; /* Rank of matrix */
1053 double X2;
1054 double Prod;
1055 long *Loca; /* Redirections in merging Chi2 classes */
1056 long *Count; /* Observed numbers */
1057 double *NbExp; /* Expected numbers */
1058 bitset_BitSet M[lmax]; /* Matrix */
1059 double V[1]; /* Number of degrees of freedom for Chi2 */
1060 char str[LENGTH + 1];
1061 lebool localRes = FALSE;
1062 chrono_Chrono *Timer;
1063 char *TestName = "smarsa_MatrixRank test";
1064
1065 Timer = chrono_Create ();
1066 /* We shall need c + 1 random numbers to build a line of the matrix */
1067 c = k / s;
1068 b = k % s;
1069 a = num_TwoExp[b];
1070 d = num_TwoExp[s];
1071 if (swrite_Basic)
1072 WriteDataMatRank (gen, TestName, N, n, r, s, l, k);
1073 if (k <= l)
1074 Minkl = k;
1075 else
1076 Minkl = l;
1077 if (res == NULL) {
1078 localRes = TRUE;
1079 res = sres_CreateChi2 ();
1080 }
1081 sres_InitChi2 (res, N, Minkl, "smarsa_MatrixRank");
1082 NbExp = res->NbExp;
1083 Count = res->Count;
1084 Loca = res->Loc;
1085
1086 Prod = n * pow (2.0, -(double) (l * k));
1087 NbExp[0] = Prod;
1088 for (j = 1; j <= Minkl; j++) {
1089 Prod = Prod * pow (2.0, (double) (l + k - 2*j + 1)) *
1090 (1.0 - 1.0 / num_TwoExp[l - j + 1]) *
1091 (1.0 - 1.0 / num_TwoExp[k - j + 1]) /
1092 (1.0 - 1.0 / num_TwoExp[j]);
1093 NbExp[j] = Prod;
1094 }
1095
1096 jlow = 0;
1097 jhigh = Minkl;
1098 if (swrite_Classes)
1099 gofs_WriteClasses (NbExp, Loca, jlow, jhigh, 0);
1100 gofs_MergeClasses (NbExp, Loca, &jlow, &jhigh, &NbGroups);
1101 if (swrite_Classes)
1102 gofs_WriteClasses (NbExp, Loca, jlow, jhigh, NbGroups);
1103 res->jmin = jlow;
1104 res->jmax = jhigh;
1105 res->degFree = NbGroups - 1;
1106
1107 util_Assert (n > 2.0 * gofs_MinExpected,
1108 "smarsa_MatrixRank: n <= 2*gofs_MinExpected");
1109 util_Assert (k <= 31, "smarsa_MatrixRank: k > 31");
1110 util_Assert (l <= lmax, "smarsa_MatrixRank: L > 64");
1111 util_Assert (l * k <= 1020, "smarsa_MatrixRank: L*k > 1020");
1112 util_Assert (NbGroups > 1,
1113 "smarsa_MatrixRank: number of classes = 1."
1114 " Increase n or decrease |L - k|");
1115
1116 sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
1117 " of freedom):", NbGroups - 1);
1118 statcoll_SetDesc (res->sVal1, str);
1119
1120 for (Seq = 1; Seq <= N; Seq++) {
1121 for (i = jlow; i <= jhigh; i++)
1122 Count[i] = 0;
1123 for (Rep = 1; Rep <= n; Rep++) {
1124 /* Generate the l x k matrix and compute its rank */
1125 for (i = 0; i < l; i++) {
1126 /* Build one line of bits L */
1127 L = 0;
1128 for (j = 1; j <= c; j++)
1129 /* Generate s bits */
1130 L = d * L + unif01_StripB (gen, r, s);
1131 /* The last b bits of a line of the matrix */
1132 if (a > 1)
1133 L = a * L + unif01_StripB (gen, r, b);
1134 M[i] = L;
1135 }
1136 /* Set all remaining lines to 0 */
1137 for (i = l; i < lmax; i++)
1138 M[i] = 0;
1139
1140 Rank = RankOfBitMatrix (M, lmax);
1141 ++Count[Loca[Rank]];
1142 }
1143
1144 X2 = gofs_Chi2 (NbExp, Count, jlow, jhigh);
1145 statcoll_AddObs (res->sVal1, X2);
1146 if (swrite_Counters)
1147 tables_WriteTabL (Count, jlow, jhigh, 5, 12, "Observed Numbers");
1148 }
1149
1150 V[0] = NbGroups - 1;
1151 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
1152 res->sVal2, res->pVal2);
1153 res->pVal1->NObs = N;
1154 sres_GetChi2SumStat (res);
1155
1156 if (swrite_Collectors)
1157 statcoll_Write (res->sVal1, 5, 14, 4, 3);
1158
1159 /* !!!! Attention, this Write must use the right pVal */
1160 if (swrite_Basic) {
1161 swrite_AddStrChi (str, LENGTH + 1, res->degFree);
1162 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
1163 swrite_Chi2SumTest (N, res);
1164 swrite_Final (gen, Timer);
1165 }
1166 if (localRes)
1167 sres_DeleteChi2 (res);
1168 chrono_Delete (Timer);
1169 }
1170
1171 #endif
1172
1173 /*=========================================================================*/
1174 #if 0
1175 static void ZeroMat (Matrix * M)
1176 {
1177 int i;
1178
1179 for (i = 0; i < M->nblignes; i++)
1180 PutBVToZero (&(M->lignes[i][0]));
1181 }
1182 #endif
1183
1184 /*-------------------------------------------------------------------------*/
1185
smarsa_MatrixRank(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,int s,int l,int k)1186 void smarsa_MatrixRank (unif01_Gen *gen, sres_Chi2 *res,
1187 long N, long n, int r, int s, int l, int k)
1188 {
1189 long Seq;
1190 long Rep;
1191 int j;
1192 int i;
1193 int c; /* Number-1 of U01 used to build a line */
1194 int b; /* Number of bits of last U01 of a line */
1195 unsigned long bmask; /* b bits mask */
1196 unsigned long smask; /* s bits mask */
1197 int Minkl; /* Min (k, l) */
1198 long NbGroups; /* Number of classes for ChiSquare */
1199 long jhigh; /* Index of the highest class */
1200 long jlow; /* Index of the lowest class */
1201 int Rank; /* Rank of matrix */
1202 double X2;
1203 double temp;
1204 long *Loca; /* Redirections in merging Chi2 classes */
1205 long *Count; /* Observed numbers */
1206 double *NbExp; /* Expected numbers */
1207 double Par[1]; /* Number of degrees of freedom for Chi2 */
1208 char str[LENGTH + 1];
1209 lebool localRes = FALSE;
1210 chrono_Chrono *Timer;
1211 char *TestName = "smarsa_MatrixRank test";
1212 Matrix *M;
1213 BitVect *V;
1214
1215 Timer = chrono_Create ();
1216 /* We shall need ceiling(c) random numbers to build a line of the matrix */
1217 c = k / s;
1218 b = k % s;
1219 bmask = num_TwoExp[b] - 1.0;
1220 /* The b most significant bits are set */
1221 bmask <<= vectorsF2_WL - b;
1222 smask = num_TwoExp[s] - 1.0;
1223 /* The s most significant bits are set */
1224 smask <<= vectorsF2_WL - s;
1225
1226 if (swrite_Basic)
1227 WriteDataMatRank (gen, TestName, N, n, r, s, l, k);
1228 Minkl = util_Min (k, l);
1229 if (res == NULL) {
1230 localRes = TRUE;
1231 res = sres_CreateChi2 ();
1232 }
1233 sres_InitChi2 (res, N, Minkl, "smarsa_MatrixRank");
1234 NbExp = res->NbExp;
1235 Count = res->Count;
1236 Loca = res->Loc;
1237
1238 temp = num_Log2((double) n) - l * k;
1239 NbExp[0] = pow (2.0, temp);
1240 for (j = 1; j <= Minkl; j++) {
1241 temp += l + k - 2*j + 1 +
1242 num_Log2(1.0 - pow (2.0, -(double) (l - j + 1))) +
1243 num_Log2(1.0 - pow (2.0, -(double) (k - j + 1))) -
1244 num_Log2(1.0 - pow (2.0, -(double) j));
1245 NbExp[j] = pow (2.0, temp);
1246 }
1247
1248 jlow = 0;
1249 jhigh = Minkl;
1250 if (swrite_Classes)
1251 gofs_WriteClasses (NbExp, Loca, jlow, jhigh, 0);
1252 gofs_MergeClasses (NbExp, Loca, &jlow, &jhigh, &NbGroups);
1253 if (swrite_Classes)
1254 gofs_WriteClasses (NbExp, Loca, jlow, jhigh, NbGroups);
1255 res->jmin = jlow;
1256 res->jmax = jhigh;
1257 res->degFree = NbGroups - 1;
1258
1259 util_Warning (NbGroups <= 1,
1260 "smarsa_MatrixRank: number of Chi2 classes = 1.\n"
1261 " Increase n or decrease |L - k|.");
1262 if (res->degFree < 1) {
1263 if (localRes)
1264 sres_DeleteChi2 (res);
1265 return;
1266 }
1267 util_Assert (n >= 2.0 * gofs_MinExpected,
1268 "smarsa_MatrixRank: n <= 2*gofs_MinExpected");
1269
1270 sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
1271 " of freedom):", NbGroups - 1);
1272 statcoll_SetDesc (res->sVal1, str);
1273
1274 M = util_Malloc (sizeof (Matrix));
1275 AllocMat (M, l, k, 1);
1276
1277 for (Seq = 1; Seq <= N; Seq++) {
1278 for (i = jlow; i <= jhigh; i++)
1279 Count[i] = 0;
1280
1281 for (Rep = 1; Rep <= n; Rep++) {
1282 /* Generate the l x k matrix and compute its rank */
1283 for (i = 0; i < l; i++) {
1284 V = &(M->lignes[i][0]);
1285 /* Build one line of bits */
1286 for (j = 0; j < c; j++) {
1287 /* Shift by s and generate s new bits */
1288 BVRShiftSelf (V, s);
1289 V->vect[0] |= (smask &
1290 (gen->GetBits (gen->param, gen->state) << r));
1291 }
1292 /* The last b bits of a line of the matrix */
1293 if (b > 0) {
1294 BVRShiftSelf (V, b);
1295 V->vect[0] |= (bmask &
1296 (gen->GetBits (gen->param, gen->state) << r));
1297 }
1298 }
1299 Rank = GaussianElimination (M, l, k, 1);
1300 ++Count[Loca[Rank]];
1301 }
1302
1303 X2 = gofs_Chi2 (NbExp, Count, jlow, jhigh);
1304 statcoll_AddObs (res->sVal1, X2);
1305 if (swrite_Counters)
1306 tables_WriteTabL (Count, jlow, jhigh, 5, 12, "Observed Numbers");
1307 }
1308
1309 FreeMat (M);
1310 util_Free (M);
1311
1312 Par[0] = NbGroups - 1;
1313 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, Par,
1314 res->sVal2, res->pVal2);
1315 res->pVal1->NObs = N;
1316 sres_GetChi2SumStat (res);
1317
1318 if (swrite_Collectors)
1319 statcoll_Write (res->sVal1, 5, 14, 4, 3);
1320
1321 /* !!!! Attention, this Write must use the right pVal */
1322 if (swrite_Basic) {
1323 swrite_AddStrChi (str, LENGTH, res->degFree);
1324 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
1325 swrite_Chi2SumTest (N, res);
1326 swrite_Final (gen, Timer);
1327 }
1328 if (localRes)
1329 sres_DeleteChi2 (res);
1330 chrono_Delete (Timer);
1331 }
1332
1333
1334 /*=========================================================================*/
1335
WriteDataSavir2(unif01_Gen * gen,char * TestName,long N,long n,int r,long m,int t)1336 static void WriteDataSavir2 (unif01_Gen * gen, char *TestName,
1337 long N, long n, int r, long m, int t)
1338 {
1339 swrite_Head (gen, TestName, N, n, r);
1340 printf (", m = %1ld, t = %1d\n\n", m, t);
1341 }
1342
1343
1344 /*-------------------------------------------------------------------------*/
1345
smarsa_Savir2(unif01_Gen * gen,sres_Chi2 * res,long N,long n,int r,long m,int t)1346 void smarsa_Savir2 (unif01_Gen *gen, sres_Chi2 *res,
1347 long N, long n, int r, long m, int t)
1348 {
1349 const double eps = 1.0E-15;
1350 long I;
1351 long msup; /* Dimension - 1 of arrays */
1352 long Seq;
1353 long Rep;
1354 long j;
1355 int i;
1356 long NbGroups; /* Number of classes for ChiSquare */
1357 long jhigh; /* Index of the highest class */
1358 long jlow; /* Index of the lowest class */
1359 double X2; /* ChiSquare Statistic */
1360 double UnSurm = 1.0 / m;
1361 double *Prob; /* Probabilities */
1362 long *Loca;
1363 double V[1]; /* Number degrees of freedom for Chi2 */
1364 char str[LENGTH + 1];
1365 lebool localRes = FALSE;
1366 chrono_Chrono *Timer;
1367 char *TestName = "smarsa_Savir2 test";
1368
1369 Timer = chrono_Create ();
1370 if (swrite_Basic)
1371 WriteDataSavir2 (gen, TestName, N, n, r, m, t);
1372
1373 Prob = util_Calloc ((size_t) m + 2, sizeof (double));
1374 Prob[m + 1] = 0.0;
1375 for (j = 1; j <= m; j++)
1376 Prob[j] = UnSurm;
1377 for (i = 2; i <= t; i++) {
1378 for (j = m; j >= 1; j--)
1379 Prob[j] = Prob[j + 1] + Prob[j] / j;
1380 }
1381 j = 1;
1382 while (Prob[j] > eps)
1383 ++j;
1384 msup = j - 1;
1385
1386 if (res == NULL) {
1387 localRes = TRUE;
1388 res = sres_CreateChi2 ();
1389 }
1390 sres_InitChi2 (res, N, msup, "smarsa_Savir2");
1391
1392 for (j = 1; j <= msup; j++)
1393 res->NbExp[j] = Prob[j] * n;
1394 util_Free (Prob);
1395 Loca = res->Loc;
1396
1397 jlow = 1;
1398 jhigh = msup;
1399 if (swrite_Classes)
1400 gofs_WriteClasses (res->NbExp, Loca, jlow, jhigh, 0);
1401 gofs_MergeClasses (res->NbExp, Loca, &jlow, &jhigh, &NbGroups);
1402 if (swrite_Classes)
1403 gofs_WriteClasses (res->NbExp, Loca, jlow, jhigh, NbGroups);
1404 res->jmin = jlow;
1405 res->jmax = jhigh;
1406 res->degFree = NbGroups - 1;
1407
1408 util_Warning (NbGroups < 2,
1409 "smarsa_Savir2: Number of classes = 1.\n Decrease t or increase n.");
1410 if (res->degFree < 1) {
1411 if (localRes)
1412 sres_DeleteChi2 (res);
1413 return;
1414 }
1415 util_Assert (n >= 2.0 * gofs_MinExpected,
1416 "smarsa_Savir2: n <= 2*gofs_MinExpected");
1417
1418 sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
1419 " of freedom):", NbGroups - 1);
1420 res->sVal1 = statcoll_Create (N, str);
1421
1422 for (Seq = 1; Seq <= N; Seq++) {
1423 for (j = jlow; j <= jhigh; j++)
1424 res->Count[j] = 0;
1425 for (Rep = 1; Rep <= n; Rep++) {
1426 I = m;
1427 for (i = 1; i <= t; i++)
1428 I = 1 + unif01_StripD (gen, r) * I;
1429 if (I > msup)
1430 ++res->Count[Loca[msup]];
1431 else
1432 ++res->Count[Loca[I]];
1433 }
1434
1435 if (swrite_Counters)
1436 tables_WriteTabL (res->Count, jlow, jhigh, 5, 12,
1437 "Observed Numbers");
1438
1439 X2 = gofs_Chi2 (res->NbExp, res->Count, jlow, jhigh);
1440 statcoll_AddObs (res->sVal1, X2);
1441 }
1442
1443 V[0] = NbGroups - 1;
1444 gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
1445 res->sVal2, res->pVal2);
1446 res->pVal1->NObs = N;
1447 sres_GetChi2SumStat (res);
1448
1449 if (swrite_Collectors)
1450 statcoll_Write (res->sVal1, 5, 14, 4, 3);
1451 if (swrite_Basic) {
1452 swrite_AddStrChi (str, LENGTH, res->degFree);
1453 gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
1454 swrite_Chi2SumTest (N, res);
1455 swrite_Final (gen, Timer);
1456 }
1457 if (localRes)
1458 sres_DeleteChi2 (res);
1459 chrono_Delete (Timer);
1460 }
1461
1462
1463 /*=========================================================================*/
1464
WriteDataGCD(unif01_Gen * gen,char * TestName,long N,long n,int r,int s)1465 static void WriteDataGCD (unif01_Gen * gen, char *TestName,
1466 long N, long n, int r, int s)
1467 {
1468 swrite_Head (gen, TestName, N, n, r);
1469 printf (", s = %1d\n\n", s);
1470 }
1471
1472
1473 /*-------------------------------------------------------------------------*/
1474
smarsa_GCD(unif01_Gen * gen,smarsa_Res2 * res,long N,long n,int r,int s)1475 void smarsa_GCD (unif01_Gen *gen, smarsa_Res2 *res,
1476 long N, long n, int r, int s)
1477 {
1478 /*
1479 The theoretical distribution for the number of iterations is unknown.
1480 The binomial is a very rough approximation: thus the printing of the
1481 results is commented out.
1482 */
1483 const double C1 = 6 / (num_Pi * num_Pi);
1484 const double P1 = 0.376;
1485 const int KMAX = 50;
1486 unsigned long U, V, temp;
1487 double X;
1488 double Param[1];
1489 char str[LENGTH + 1];
1490 lebool localRes = FALSE;
1491 chrono_Chrono *Timer;
1492 char *TestName = "smarsa_GCD test";
1493 sres_Chi2 *GCD;
1494 sres_Chi2 *NumIter;
1495 int jmax, j, k;
1496 long Seq, i;
1497 double *NbExp;
1498 long *Loc;
1499 long NbClasses;
1500 fmass_INFO Q;
1501
1502 Timer = chrono_Create ();
1503 if (swrite_Basic)
1504 WriteDataGCD (gen, TestName, N, n, r, s);
1505 if (n < 30) {
1506 util_Warning (TRUE, "n < 30");
1507 return;
1508 }
1509 if (n > pow (2.0, 1.5*s)) {
1510 util_Warning (TRUE, "n > 2^(1.5s)");
1511 return;
1512 }
1513 if (res == NULL) {
1514 localRes = TRUE;
1515 res = smarsa_CreateRes2 ();
1516 }
1517 jmax = 1 + sqrt (C1 * n / gofs_MinExpected);
1518 util_Assert (jmax > 1, "smarsa_GCD: jmax < 2");
1519 InitRes2 (res, N, jmax, KMAX);
1520
1521 GCD = res->GCD;
1522 GCD->jmin = 1;
1523 GCD->jmax = jmax;
1524 GCD->degFree = jmax - 1;
1525 sprintf (str, "GCD; the N statistic values (a ChiSquare with %1d degrees"
1526 " of freedom):", jmax - 1);
1527 statcoll_SetDesc (GCD->sVal1, str);
1528
1529 /* Compute the probabilities for the GCD values */
1530 NbExp = GCD->NbExp;
1531 Loc = GCD->Loc;
1532 X = 0.0;
1533 for (j = 1; j < jmax; j++) {
1534 NbExp[j] = n * C1 / ((double) j * j);
1535 X += NbExp[j];
1536 Loc[j] = j;
1537 }
1538 NbExp[jmax] = n - X;
1539
1540 if (swrite_Classes) {
1541 printf ("Classes for the GCD values:\n");
1542 gofs_WriteClasses (GCD->NbExp, GCD->Count, 1, jmax, 0);
1543 }
1544
1545 NumIter = res->NumIter;
1546 /* Compute expected numbers for number of iterations */
1547 Q = fmass_CreateBinomial (KMAX, P1, 1.0 - P1);
1548 for (i = 0; i <= KMAX; i++)
1549 NumIter->NbExp[i] = n * fmass_BinomialTerm2 (Q, i);
1550 fmass_DeleteBinomial (Q);
1551
1552 NumIter->jmin = 0;
1553 NumIter->jmax = KMAX;
1554 if (swrite_Classes) {
1555 printf ("\nClasses for the number of iterations:\n");
1556 gofs_WriteClasses (NumIter->NbExp, NumIter->Loc, NumIter->jmin,
1557 NumIter->jmax, 0);
1558 }
1559 gofs_MergeClasses (NumIter->NbExp, NumIter->Loc, &NumIter->jmin,
1560 &NumIter->jmax, &NbClasses);
1561
1562 if (swrite_Classes)
1563 gofs_WriteClasses (NumIter->NbExp, NumIter->Loc, NumIter->jmin,
1564 NumIter->jmax, NbClasses);
1565
1566 sprintf (str, "NumIter; the N statistic values (a ChiSquare with %1ld"
1567 " degrees of freedom):", NbClasses - 1);
1568 statcoll_SetDesc (NumIter->sVal1, str);
1569 NumIter->degFree = NbClasses - 1;
1570 util_Assert (NumIter->degFree >= 1, "NumIter->degFree < 1");
1571
1572 for (Seq = 1; Seq <= N; Seq++) {
1573 for (i = 0; i <= KMAX; i++)
1574 NumIter->Count[i] = 0;
1575 for (i = 0; i <= GCD->jmax; i++)
1576 GCD->Count[i] = 0;
1577 for (i = 1; i <= n; i++) {
1578 k = 0;
1579 do {
1580 U = unif01_StripB (gen, r, s);
1581 V = unif01_StripB (gen, r, s);
1582 } while (0 == U || 0 == V);
1583 do {
1584 temp = U % V;
1585 U = V;
1586 V = temp;
1587 k++;
1588 } while (V > 0);
1589 if ((long) U > GCD->jmax)
1590 U = GCD->jmax;
1591 (GCD->Count[U])++;
1592 if (k > KMAX)
1593 k = KMAX;
1594 (NumIter->Count[NumIter->Loc[k]])++;
1595 }
1596 if (swrite_Counters) {
1597 tables_WriteTabL (GCD->Count, GCD->jmin, GCD->jmax, 5, 10,
1598 "Observed numbers for GCD values:");
1599 /* tables_WriteTabL (NumIter->Count, NumIter->jmin, NumIter->jmax, 5,
1600 10, "Observed numbers for number of iterations:");
1601 */
1602 }
1603
1604 X = gofs_Chi2 (GCD->NbExp, GCD->Count, GCD->jmin, GCD->jmax);
1605 statcoll_AddObs (GCD->sVal1, X);
1606 X = gofs_Chi2 (NumIter->NbExp, NumIter->Count, NumIter->jmin,
1607 NumIter->jmax);
1608 statcoll_AddObs (NumIter->sVal1, X);
1609 }
1610
1611 Param[0] = GCD->degFree;
1612 gofw_ActiveTests2 (GCD->sVal1->V, GCD->pVal1->V, N, wdist_ChiSquare,
1613 Param, GCD->sVal2, GCD->pVal2);
1614 GCD->pVal1->NObs = N;
1615 sres_GetChi2SumStat (GCD);
1616 /*
1617 Param[0] = NumIter->degFree;
1618 gofw_ActiveTests2 (NumIter->sVal1->V, NumIter->pVal1->V, N,
1619 wdist_ChiSquare, Param, NumIter->sVal2, NumIter->pVal2);
1620 NumIter->pVal1->NObs = N;
1621 */
1622
1623 if (swrite_Basic) {
1624 if (swrite_Collectors)
1625 statcoll_Write (GCD->sVal1, 5, 14, 4, 3);
1626 printf ("\n-----------------------------------------------\n");
1627 if (N == 1) {
1628 printf ("Number of degrees of freedom : %4ld\n",
1629 GCD->degFree);
1630 printf ("Chi2 statistic for GCD values :");
1631 gofw_Writep2 (GCD->sVal2[gofw_Mean], GCD->pVal2[gofw_Mean]);
1632 } else {
1633 printf ("Test results for GCD values:\n");
1634 gofw_WriteActiveTests0 (N, GCD->sVal2, GCD->pVal2);
1635 swrite_Chi2SumTest (N, GCD);
1636 }
1637 /*
1638 if (swrite_Collectors)
1639 statcoll_Write (NumIter->sVal1, 5, 14, 4, 3);
1640 printf ("\n-----------------------------------------------\n");
1641 if (N == 1) {
1642 printf ("Number of degrees of freedom : %4ld\n",
1643 NumIter->degFree);
1644 printf ("Chi2 statistic for NumIter :");
1645 gofw_Writep2 (NumIter->sVal2[gofw_Mean], NumIter->pVal2[gofw_Mean]);
1646 } else {
1647 printf ("Test results for NumIter:\n");
1648 gofw_WriteActiveTests0 (N, NumIter->sVal2, NumIter->pVal2);
1649 swrite_SumTest (N, NumIter->sVal2[gofw_Sum], NumIter->pVal2[gofw_Sum],
1650 N*NumIter->degFree);
1651 }
1652 */
1653 printf ("\n\n");
1654 swrite_Final (gen, Timer);
1655 }
1656
1657 if (localRes)
1658 smarsa_DeleteRes2 (res);
1659 chrono_Delete (Timer);
1660
1661 }
1662
1663
1664 /*=========================================================================*/
1665