htestu-0.1.0.0: cbits/testu/src/sstring.c
/*************************************************************************\
*
* Package: TestU01
* File: sstring.c
* Environment: ANSI C
*
* Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
* e-mail: lecuyer@iro.umontreal.ca
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted without a fee for private, research,
* academic, or other non-commercial purposes.
* Any use of this software in a commercial environment requires a
* written licence from the copyright owner.
*
* Any changes made to this package must be clearly identified as such.
*
* In scientific publications which used this software, a reference to it
* would be appreciated.
*
* Redistributions of source code must retain this copyright notice
* and the following disclaimer.
*
* THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
\*************************************************************************/
#include "util.h"
#include "chrono.h"
#include "num.h"
#include "tables.h"
#include "bitset.h"
#include "sstring.h"
#include "unif01.h"
#include "wdist.h"
#include "swrite.h"
#include "sres.h"
#include "gofs.h"
#include "gofw.h"
#include "fbar.h"
#include "statcoll.h"
#include <math.h>
#include <float.h>
#include <limits.h>
#include <stdio.h>
#include <string.h>
/*------------------------------ Constants --------------------------------*/
/* Minimal length (number of bits) of a sequence for LongestHeadRun */
#define LMIN 1000
/* Max string length for the correlations in PeriodsInStrings */
#define MAX_CORR 31
/* Max dimension of arrays */
#define DIM 1000
/* Max string lengths */
#define LEN1 200
#define LEN2 200
/*-------------------------------- Types ----------------------------------*/
typedef struct InfoListC *ListC; /* A correlation list */
struct InfoListC {
long Nb; /* Number of bits of a correlation in
the initial computations; then
number of occurences */
bitset_BitSet C; /* A correlation */
long Pop; /* Population related to C (and c) */
ListC Ext; /* The smallest extension of C */
ListC Ext0; /* The smallest extension of D longer
than C, if C is an extension of D.
Initially NULL */
ListC Next; /* Next correlation of same length */
};
/* Corr contains the lists of correlations of each length for s in [0..smax]
*/
typedef struct {
ListC Corr[MAX_CORR + 1];
int smax;
} sstring_Corr;
/*----------------------------- Variables --------------------------------*/
lebool sstring_CorrFlag = FALSE;
lebool sstring_Counters = FALSE;
/*----------------------------- Functions --------------------------------*/
static void InitRes3 (
sstring_Res3 *res, /* Results holder */
long N, /* Number of replications */
int jmax /* Max class index for chi2 */
)
/*
* Initializes the sstring_Res3 structure
*/
{
sres_InitBasic (res->NBits, N, "sstring_Run: Number of Bits");
sres_InitChi2 (res->NRuns, N, jmax, "sstring_Run: Number of Runs");
res->Count0 = util_Realloc (res->Count0, (jmax + 1) * sizeof (long));
res->Count1 = util_Realloc (res->Count1, (jmax + 1) * sizeof (long));
res->NRuns->jmin = 1;
res->NRuns->degFree = jmax - 1;
}
/*-------------------------------------------------------------------------*/
sstring_Res3 * sstring_CreateRes3 (void)
{
sstring_Res3 *res;
res = util_Malloc (sizeof (sstring_Res3));
res->NBits = sres_CreateBasic ();
res->NRuns = sres_CreateChi2 ();
res->Count0 = util_Calloc (1, sizeof (long));
res->Count1 = util_Calloc (1, sizeof (long));
return res;
}
/*-------------------------------------------------------------------------*/
void sstring_DeleteRes3 (sstring_Res3 *res)
{
if (res == NULL)
return;
res->Count0 = util_Free (res->Count0);
res->Count1 = util_Free (res->Count1);
sres_DeleteBasic (res->NBits);
sres_DeleteChi2 (res->NRuns);
util_Free (res);
}
/*=========================================================================*/
static void InitRes2 (
sstring_Res2 *res, /* Results holder */
long N, /* Number of replications */
int jhigh /* Max class index for chi2 */
)
/*
* Initializes the sstring_Res2 structure
*/
{
sres_InitDisc (res->Disc, N,
"sstring_LongestHeadRun: Global longest run of 1's");
sres_InitChi2 (res->Chi, N, jhigh,
"sstring_LongestHeadRun: Block longest runs of 1's");
}
/*-------------------------------------------------------------------------*/
sstring_Res2 * sstring_CreateRes2 (void)
{
sstring_Res2 *res;
res = util_Malloc (sizeof (sstring_Res2));
res->Chi = sres_CreateChi2 ();
res->Disc = sres_CreateDisc ();
return res;
}
/*-------------------------------------------------------------------------*/
void sstring_DeleteRes2 (sstring_Res2 *res)
{
if (res == NULL)
return;
sres_DeleteChi2 (res->Chi);
sres_DeleteDisc (res->Disc);
util_Free (res);
}
/*=========================================================================*/
static void InitRes (
sstring_Res *res, /* Results holder */
long N, /* Number of replications */
int L, /* Size of blocks (number of bits) */
int d, /* Parameter for sub-matrices */
char *nam
)
/*
* Initializes res
*/
{
int i;
sres_InitBasic (res->Bas, N, nam);
if (res->L > 0) {
tables_DeleteMatrixL (&res->Counters);
tables_DeleteMatrixD (&res->ZCounters);
}
res->Counters = tables_CreateMatrixL (L + 2, L + 1);
res->ZCounters = tables_CreateMatrixD (L + 2, L + 1);
if (d < 0)
d = 0;
for (i = d + 1; i <= res->d; i++)
sres_DeleteBasic (res->Block[i]);
for (i = res->d + 1; i <= d; i++)
res->Block[i] = sres_CreateBasic ();
for (i = 1; i <= d; i++)
sres_InitBasic (res->Block[i], N, nam);
res->L = L;
res->d = d;
}
/*-------------------------------------------------------------------------*/
sstring_Res * sstring_CreateRes (void)
{
sstring_Res *res;
res = util_Malloc (sizeof (sstring_Res));
memset (res, 0, sizeof (sstring_Res));
res->Bas = sres_CreateBasic ();
res->Style = tables_Plain;
res->L = -1;
res->d = 0;
return res;
}
/*-------------------------------------------------------------------------*/
void sstring_DeleteRes (sstring_Res *res)
{
if (res == NULL)
return;
if (res->d > 0) {
int i;
for (i = 1; i <= res->d; i++) {
sres_DeleteBasic (res->Block[i]);
}
}
if (res->L > 0) {
tables_DeleteMatrixD (&res->ZCounters);
tables_DeleteMatrixL (&res->Counters);
}
sres_DeleteBasic (res->Bas);
util_Free (res);
}
/*=========================================================================*/
static long Psi (
bitset_BitSet C, /* Correlation making up tail of correlation k */
long j, /* Length of correlation C */
long k /* Correlation made up of 1 followed by 0's
until C, i.e. k = 100...000C */
)
{
/* j <=> c & k <=> k in the article */
if (k > j)
return 0;
if (k <= 0)
return (long) num_TwoExp[-k];
if (bitset_TestBit (C, j - k))
return 1;
else
return 0;
}
/*-------------------------------------------------------------------------*/
static void DeleteCorr (sstring_Corr *corr)
/*
* Delete all correlations.
*/
{
ListC Ci, OldCi;
int i;
if (corr == NULL)
return;
for (i = 0; i <= corr->smax; i++) {
Ci = corr->Corr[i];
while (Ci) {
OldCi = Ci;
Ci = Ci->Next;
util_Free (OldCi);
}
}
util_Free (corr);
}
/*-------------------------------------------------------------------------*/
static sstring_Corr * CreateCorr (int s)
/*
* Compute all possible correlations for strings of length s
*/
{
ListC CjE, Cj, Ci, OldCi, XS;
sstring_Corr *corr;
int j, i, k, Tmax;
long p;
corr = util_Malloc (sizeof (sstring_Corr));
memset (corr, 0, sizeof (sstring_Corr));
corr->smax = s;
XS = corr->Corr[0] = util_Malloc (sizeof (struct InfoListC));
XS->Nb = 0;
XS->Pop = 1;
XS->Ext = NULL;
XS->Ext0 = NULL;
XS->Next = NULL;
XS = corr->Corr[1] = util_Malloc (sizeof (struct InfoListC));
XS->Nb = 1;
bitset_SetBit (XS->C, 0);
XS->Pop = 2;
XS->Ext = NULL;
XS->Ext0 = NULL;
XS->Next = NULL;
for (i = 2; i <= s; i++) { /* i is the string length */
/* Count and build the list of correlations of length i. */
Ci = corr->Corr[i] = util_Malloc (sizeof (struct InfoListC));
for (j = 0; j <= i - 2; j++) {
/* j is the length of correlation C in 100...00C */
Cj = corr->Corr[j];
while (Cj) {
/* Compute the number of strings of length i and with corre- */
/* lation "10...0C", and if > 0, add this corr. to Corr[i] */
p = Cj->Pop * Psi (Cj->C, j, 2*j - i);
/* Check if 1C may be a correlation. Possible only if
C = 111....1, i.e. the last of the list */
if (Cj->Next == NULL)
p -= 2 * Psi (Cj->C, j, 2*j + 2 - i);
CjE = Cj->Ext;
/* Note: j <= j-2, i.e. j+1 <= (i+j) / 2 */
Tmax = (i + j) / 2;
while (CjE && CjE->Nb <= Tmax) {
p -= CjE->Pop * Psi (Cj->C, j, 2*CjE->Nb - i);
CjE = CjE->Ext0;
}
/* p = number of strings looked for */
if (p > 0) {
/* Put this correlation in Ci */
Ci->Nb = i;
Ci->Pop = p;
Ci->Ext = NULL;
Ci->Ext0 = NULL;
/* Ci->C becomes Cj->C shifted right by i-j */
/* positions, with a 1 in first position. */
Ci->C = 0;
bitset_SetBit (Ci->C, 0);
if (j > 0) {
for (k = 0; k < j; k++) {
if (bitset_TestBit (Cj->C, k)) {
bitset_SetBit (Ci->C, k + i - j);
}
}
}
if (Cj->Ext == NULL)
Cj->Ext = Ci;
else {
CjE = Cj->Ext;
while (CjE->Ext0)
CjE = CjE->Ext0;
CjE->Ext0 = Ci;
}
OldCi = Ci;
Ci = util_Malloc (sizeof (struct InfoListC));
OldCi->Next = Ci;
}
Cj = Cj->Next;
}
}
/* For j = i-1, we have the correlation "11...1" */
Ci->C = 0;
for (k = 0; k < i; k++) {
bitset_SetBit (Ci->C, k);
}
Ci->Nb = i;
Ci->Pop = 2;
Ci->Ext = NULL;
Ci->Ext0 = NULL;
Ci->Next = NULL;
}
return corr;
}
/*=========================================================================*/
static void sstring_WriteCorr (sstring_Corr *corr, int s)
{
ListC Cs;
int k;
char str [LEN1 + 1];
if (corr == NULL) {
util_Warning (TRUE,
"sstring_WriteCorr: corr is a NULL pointer");
return;
}
if (corr->smax < s) {
sprintf (str, "sstring_WriteCorr: invalid s = %d", s);
util_Error (str);
}
Cs = corr->Corr[s];
if (Cs == NULL)
return;
printf ("\n-----------------------------------------------------\n"
"List of correlations of length %d and their population\n\n", s);
while (Cs) {
for (k = 0; k < s; k++) {
if (bitset_TestBit (Cs->C, k))
printf ("1");
else
printf ("0");
}
printf ("%12ld\n", Cs->Pop);
Cs = Cs->Next;
}
printf ("\n\n");
}
/*-------------------------------------------------------------------------*/
static bitset_BitSet GenerateC (
unif01_Gen *gen, /* Generator */
int r, /* Drop first r bits of each random number */
int s /* Keep next s bits of each random number */
)
/*
* Generate a string of s bits and return its correlation. To determine the
* correlation of a bit string, compare strings g and d (initially equal).
* For each iteration:
* 1) drop leftmost bit of string g
* 2) drop rightmost bit of string d
* 3) if g = d, bit k of the correlation c is 1, otherwise 0.
*/
{
int k;
unsigned long g, d, lbit = s - 1;
bitset_BitSet c = 0;
/* Generate a random number; drop r most significant bits; keep s next
bits */
g = d = unif01_StripB (gen, r, s);
/* Initialization of correlation, trivial case */
bitset_SetBit (c, 0);
for (k = 1; k < s; k++) {
/* drop leftmost bit of string g */
bitset_ClearBit (g, lbit);
/* drop rightmost bit of string d */
d >>= 1;
/* if g = d, bit k of the correlation c is 1, otherwise 0 */
if (g == d)
bitset_SetBit (c, k);
lbit--;
}
return c;
}
/*=========================================================================*/
static void WriteDataPeriod (
unif01_Gen *gen, /* generator */
char *Test, /* Test name */
long N, /* Number of replications */
long n, /* Sample size */
int r, /* r first bits of each random number dropped */
int s /* s bits of each random number used */
)
{
swrite_Head (gen, Test, N, n, r);
printf (", s = %4d\n\n", s);
}
/*-------------------------------------------------------------------------*/
void sstring_PeriodsInStrings (unif01_Gen *gen, sres_Chi2 *res,
long N, long n, int r, int s)
{
ListC XS, Ci;
sstring_Corr *corr;
long jhigh, /* Highest class for ChiSquare */
jlow, /* Lowest class for ChiSquare */
NbGroups; /* Number of classes for ChiSquare */
long j, i;
long Seq; /* One replication of the test */
double Fraction, X2;
bitset_BitSet D; /* A correlation */
double V[1]; /* Number of ChiSquare degrees of freedom */
char str [LEN1 + 1];
double NbExp [DIM + 1];
long Loca [DIM + 1];
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_PeriodsInStrings test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataPeriod (gen, TestName, N, n, r, s);
util_Assert (r >= 0, "sstring_PeriodsInStrings: r < 0");
util_Assert (r <= 31, "sstring_PeriodsInStrings: r > 31");
util_Assert (r + s <= 31, "sstring_PeriodsInStrings: r + s > 31");
util_Assert (s <= 31, "sstring_PeriodsInStrings: s > 31");
util_Assert (s >= 2, "sstring_PeriodsInStrings: s < 2");
/* util_Assert (n > 2.0 * gofs_MinExpected,
"sstring_PeriodsInStrings: n <= 2*gofs_MinExpected"); */
Fraction = n / num_TwoExp[s];
if (res == NULL) {
localRes = TRUE;
res = sres_CreateChi2 ();
}
corr = CreateCorr (s);
if (sstring_CorrFlag)
sstring_WriteCorr (corr, s);
/* Get the expected numbers of the population count */
XS = Ci = corr->Corr[s];
j = 1;
while (Ci) {
NbExp[j] = Fraction * Ci->Pop;
Ci = Ci->Next;
++j;
util_Assert (j <= DIM, "sstring_PeriodsInStrings: DIM too small");
}
jlow = 1;
jhigh = j - 1;
if (swrite_Classes)
gofs_WriteClasses (NbExp, Loca, jlow, jhigh, 0);
/* Merge classes for the chi-square test */
gofs_MergeClasses (NbExp, Loca, &jlow, &jhigh, &NbGroups);
if (swrite_Classes)
gofs_WriteClasses (NbExp, Loca, jlow, jhigh, NbGroups);
res->degFree = NbGroups - 1;
if (res->degFree < 1) {
if (localRes)
sres_DeleteChi2 (res);
return;
}
sres_InitChi2 (res, N, jhigh, "sstring_PeriodsInStrings");
res->jmin = jlow;
tables_CopyTabD (NbExp, res->NbExp, jlow, jhigh);
tables_CopyTabL (Loca, res->Loc, jlow, jhigh);
sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
" of freedom):", NbGroups - 1);
statcoll_SetDesc (res->sVal1, str);
/* Test begins */
for (Seq = 1; Seq <= N; Seq++) {
/* Zero the population counters */
Ci = XS;
while (Ci) {
Ci->Nb = 0;
Ci = Ci->Next;
}
for (i = 1; i <= n; i++) {
D = GenerateC (gen, r, s);
/* Find the correlation */
Ci = XS;
while (Ci->C != D)
Ci = Ci->Next;
++Ci->Nb;
}
/* Keep the observed numbers in sstring_Count */
for (j = jlow; j <= jhigh; j++)
res->Count[j] = 0;
Ci = XS;
j = 1;
while (Ci) {
if (j >= res->jmax)
res->Count[res->jmax] += Ci->Nb;
else
res->Count[Loca[j]] += Ci->Nb;
Ci = Ci->Next;
++j;
}
X2 = gofs_Chi2 (res->NbExp, res->Count, jlow, jhigh);
statcoll_AddObs (res->sVal1, X2);
if (swrite_Counters)
tables_WriteTabL (res->Count, jlow, jhigh, 5, 10,
"Observed population counts");
}
res->degFree = V[0] = NbGroups - 1;
gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
res->sVal2, res->pVal2);
res->pVal1->NObs = N;
sres_GetChi2SumStat (res);
if (swrite_Collectors)
statcoll_Write (res->sVal1, 5, 14, 4, 3);
if (swrite_Basic) {
swrite_AddStrChi (str, LEN1, NbGroups - 1);
gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
swrite_Chi2SumTest (N, res);
swrite_Final (gen, Timer);
}
DeleteCorr (corr);
if (localRes)
sres_DeleteChi2 (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static double ProbabiliteLHR (long j, double Lnl)
/*
* Returns the probability that the longest series of successive 1 has
* length = j.
*/
{
double x, temp;
temp = (j + 1) * num_Ln2 - Lnl;
x = exp (-exp (-temp));
temp += num_Ln2;
x = exp (-exp (-temp)) - x;
return x;
}
/*-------------------------------------------------------------------------*/
static void WriteDataLongHead (unif01_Gen *gen, char *Test,
long N, long n, int r, int s, long L)
{
swrite_Head (gen, Test, N, n, r);
printf (", s = %1d, L = %1ld\n\n", s, L);
}
/*-------------------------------------------------------------------------*/
void sstring_LongestHeadRun (unif01_Gen *gen, sstring_Res2 *res,
long N, long n, int r, int s, long L)
{
const double eps = DBL_EPSILON;
const long K = L/s; /* Number of iterations */
long Rep; /* Current replication number */
long Seq; /* Current sequence number */
long i;
double LnLen; /* log (L) or log (NnL) */
double X2, temp;
int j;
long longest; /* Longest serie of 1 in a block */
long longest2; /* Longest serie of 1 over all blocks */
long longest3; /* Longest serie of 1 over all Replic. */
long longueur; /* Run length in a block */
long longueur2; /* Run length in a sequence */
long longueur3; /* Run length in a replication */
long jhigh; /* Highest class for Chi2 */
long jhigh2; /* Highest index for CDF[j] */
long NbGroups; /* Number of classes for Chi2 */
bitset_BitSet ensemble; /* Chosen bits in each generated number */
double V[1]; /* Number degrees of freedom for Chi2 */
char str [LEN1 + 1];
double NbExp [DIM + 1]; /* Expected numbers */
double CDF [DIM + 1]; /* Cumulative probabilities */
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_LongestHeadRun test";
sres_Chi2 *Chi;
sres_Disc *Disc;
Timer = chrono_Create ();
L = K * s;
if (swrite_Basic)
WriteDataLongHead (gen, TestName, N, n, r, s, L);
util_Assert (L >= LMIN, "sstring_LongestHeadRun: L < 1000");
if (res == NULL) {
localRes = TRUE;
res = sstring_CreateRes2 ();
}
jhigh = DIM;
/* Get the expected numbers for the chi-square for blocks of L bits */
LnLen = log ((double) L);
CDF[0] = ProbabiliteLHR (0, LnLen);
NbExp[0] = n * CDF[0];
for (j = 1; j < DIM; j++) {
temp = ProbabiliteLHR (j, LnLen);
NbExp[j] = n * temp;
CDF[j] = temp + CDF[j-1];
if ((temp <= eps) && (CDF[j] > 0.5)) {
jhigh = j;
break;
}
}
util_Assert (jhigh > 0, "sstring_LongestHeadRun: jhigh = 0");
NbExp[jhigh] = n * (1.0 - CDF[jhigh - 1]);
/* Get the probabilities for the global run over the N*n*L bits */
LnLen = log (N * (double) n * (double) L); /* Avoid overflow of long */
CDF[0] = ProbabiliteLHR (0, LnLen);
for (j = 1; j < DIM; j++) {
temp = ProbabiliteLHR (j, LnLen);
CDF[j] = temp + CDF[j-1];
if ((temp <= eps) && (CDF[j] > 0.5)) {
jhigh2 = j;
break;
}
}
InitRes2 (res, N, jhigh);
Disc = res->Disc;
Chi = res->Chi;
tables_CopyTabD (NbExp, Chi->NbExp, 0, jhigh);
if (swrite_Classes)
gofs_WriteClasses (Chi->NbExp, Chi->Loc, 0, jhigh, 0);
gofs_MergeClasses (Chi->NbExp, Chi->Loc, &Chi->jmin, &Chi->jmax, &NbGroups);
if (swrite_Classes)
gofs_WriteClasses (Chi->NbExp, Chi->Loc, Chi->jmin, Chi->jmax, NbGroups);
Chi->degFree = NbGroups - 1;
if (Chi->degFree <= 0) {
util_Assert (1, "sstring_LongestHeadRun: Chi->degFree = 0");
if (localRes)
sstring_DeleteRes2 (res);
return;
}
sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
" of freedom):", Chi->degFree);
statcoll_SetDesc (Chi->sVal1, str);
statcoll_SetDesc (Disc->sVal1,
"The longest run of 1 for each replication ");
/* Beginning of test */
longest3 = longueur3 = 0;
for (Rep = 1; Rep <= N; Rep++) {
for (i = Chi->jmin; i <= Chi->jmax; i++)
Chi->Count[i] = 0;
longest2 = -1; /* -1 at the beginning of a new replication */
longueur2 = 0;
for (Seq = 1; Seq <= n; Seq++) {
longest = -1; /* -1 at the beginning of a new sequence */
longueur = 0;
for (i = 1; i <= K; i++) {
/* Now build a block of L bits */
ensemble = unif01_StripB (gen, r, s);
/* Examine each bit of a number */
for (j = s - 1; j >= 0; j--) {
if (bitset_TestBit (ensemble, j))
++longueur;
else {
/* Beginning of a sequence: merge last block of 1's of */
/* last sequence with first block of 1's of new sequence */
if (longest < 0) {
/* Beginning of a replication: merge last sequence of */
/* 1's of last replication with first sequence of 1's */
/* of new replication */
if (longest2 < 0) {
longueur3 += longueur;
if (longueur3 > longest3)
longest3 = longueur3;
}
longueur2 += longueur;
if (longueur2 > longest2)
longest2 = longueur2;
}
if (longueur > longest)
longest = longueur;
longueur = 0;
}
}
}
if (longueur > longest)
longest = longueur;
if (longest >= Chi->jmax)
++Chi->Count[Chi->jmax];
else if (longest <= Chi->jmin)
++Chi->Count[Chi->jmin];
else
++Chi->Count[Chi->Loc[longest]];
if (longest > longest2)
longest2 = longest;
longueur3 = longueur2 = longueur;
}
X2 = gofs_Chi2 (Chi->NbExp, Chi->Count, Chi->jmin, Chi->jmax);
statcoll_AddObs (Chi->sVal1, X2);
statcoll_AddObs (Disc->sVal1, (double) longest2);
if (longest2 > longest3)
longest3 = longest2;
if (swrite_Counters)
tables_WriteTabL (Chi->Count, Chi->jmin, Chi->jmax, 5, 10,
"Observed numbers");
longueur = 0;
for (j = Chi->jmin; j <= Chi->jmax; j++)
longueur += Chi->Count[j];
util_Warning (longueur != n, "Total Count != n");
}
Disc->sVal2 = longest3;
if (longest3 > jhigh2) {
Disc->pLeft = 1.0;
Disc->pRight = 0.0;
} else {
Disc->pLeft = CDF[longest3];
if (longest3 > 0)
Disc->pRight = 1.0 - CDF[longest3 - 1];
else
Disc->pRight = 1.0;
}
Disc->pVal2 = gofw_pDisc (Disc->pLeft, Disc->pRight);
V[0] = Chi->degFree;
gofw_ActiveTests2 (Chi->sVal1->V, Chi->pVal1->V, N, wdist_ChiSquare, V,
Chi->sVal2, Chi->pVal2);
Chi->pVal1->NObs = N;
sres_GetChi2SumStat (Chi);
if (swrite_Collectors) {
statcoll_Write (Chi->sVal1, 5, 14, 4, 3);
statcoll_Write (Disc->sVal1, 5, 14, 0, 0);
}
if (swrite_Basic) {
swrite_AddStrChi (str, LEN1, Chi->degFree);
gofw_WriteActiveTests2 (N, Chi->sVal2, Chi->pVal2, str);
swrite_Chi2SumTest (N, Chi);
printf ("-----------------------------------------------\n");
printf ("Global longest run of 1 :");
gofw_Writep2 (Disc->sVal2, Disc->pVal2);
printf ("\n\n");
swrite_Final (gen, Timer);
}
if (localRes)
sstring_DeleteRes2 (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static void HammingWeight2_L (unif01_Gen * gen, sres_Basic * res,
long N, int r, int s, long L, long K)
/*
* Generate all the n bits for the HammingWeight2 test in the case L > s.
* For the last number generated in a block of L bits, we keep its first
* LMods bits and discard the other bits.
*/
{
const int LDivs = L / s; /* A block uses LDivs numbers ... */
const int LMods = L % s; /* + 1 if LMods > 0 */
const double L2 = L / 2.0;
int co, j;
long i, Seq;
unsigned long Z;
double X2;
for (Seq = 1; Seq <= N; Seq++) {
X2 = 0.0;
for (i = 0; i < K; i++) {
/* Generate a block of L bits */
co = 0;
for (j = 0; j < LDivs; j++) {
Z = unif01_StripB (gen, r, s);
while (Z > 0) { /* Count the number of 1 bits */
Z &= Z - 1; /* Clear lowest 1 bit */
++co;
}
}
/* The last bits of the block */
if (LMods > 0) {
Z = unif01_StripB (gen, r, LMods);
while (Z > 0) {
Z &= Z - 1;
++co;
}
}
X2 += (co - L2)*(co - L2);
}
X2 *= 4.0 / L;
statcoll_AddObs (res->sVal1, X2);
}
}
/*-------------------------------------------------------------------------*/
static void HammingWeight2_S (unif01_Gen * gen, sres_Basic * res,
long N, int r, int s, long L, long K)
/*
* Generate all the n bits for the HammingWeight2 test in the case L <= s.
* A number generates sDivL blocks. If s % L == 0, we use all s bits of the
* number.
*/
{
const int sDivL = s / L; /* A number generates sDivL blocks */
const long Q = K / sDivL + (K % sDivL > 0);
const unsigned long MASK = num_TwoExp[L] - 1.0;
const double L2 = L / 2.0;
int co, j;
long i, Seq;
unsigned long Z, Y;
double X2;
for (Seq = 1; Seq <= N; Seq++) {
X2 = 0.0;
for (i = 0; i < Q; i++) {
Z = unif01_StripB (gen, r, s);
/* Generate sDivL blocks of L bits */
for (j = 0; j < sDivL; j++) {
co = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++co;
}
X2 += (co - L2)*(co - L2);
Z >>= L;
}
}
X2 *= 4.0 / L;
statcoll_AddObs (res->sVal1, X2);
}
}
/*-------------------------------------------------------------------------*/
void sstring_HammingWeight2 (unif01_Gen * gen, sres_Basic * res,
long N, long n, int r, int s, long L)
{
const long K = n / L;
double sum;
double V[1]; /* Number of Chi2 degrees of freedom */
char chaine[LEN1 + 1] = "";
char str[LEN2 + 1] = "";
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_HammingWeight2 test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataLongHead (gen, TestName, N, n, r, s, L);
util_Assert (r + s <= 32, "sstring_HammingWeight2: r + s > 32");
util_Assert (L <= n, "sstring_HammingWeight2: L > n");
util_Assert (L >= 2, "sstring_HammingWeight2: L < 2");
if (res == NULL) {
localRes = TRUE;
res = sres_CreateBasic ();
}
sres_InitBasic (res, N, "sstring_HammingWeight2");
strncpy (chaine, "sVal1: a chi-square with ", (size_t) LEN1);
sprintf (str, "%ld", K);
strncat (chaine, str, (size_t) LEN2);
strncat (chaine, " degrees of freedom", (size_t) LEN1);
statcoll_SetDesc (res->sVal1, chaine);
if (L >= s)
HammingWeight2_L (gen, res, N, r, s, L, K);
else
HammingWeight2_S (gen, res, N, r, s, L, K);
V[0] = K;
gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
res->sVal2, res->pVal2);
res->pVal1->NObs = N;
sum = N * statcoll_Average (res->sVal1);
res->sVal2[gofw_Sum] = sum;
res->pVal2[gofw_Sum] = fbar_ChiSquare2 (N*K, 12, sum);
if (swrite_Collectors)
statcoll_Write (res->sVal1, 5, 14, 2, 1);
if (swrite_Basic) {
swrite_AddStrChi (str, LEN2, K);
gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
swrite_Chi2SumTestb (N, res->sVal2[gofw_Sum], res->pVal2[gofw_Sum], K);
swrite_Final (gen, Timer);
}
if (localRes)
sres_DeleteBasic (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static void HammingWeight_L (unif01_Gen * gen, sres_Chi2 * res,
long N, long n, int r, int s, long L)
/*
* Generate all the n*L bits for the HammingWeight test in the case L > s.
* For the last number generated in a block of L bits, we keep its first
* LMods bits and discard the other bits.
*/
{
const int LDivs = L / s; /* A block uses LDivs numbers ... */
const int LMods = L % s; /* + 1 if LMods > 0 */
int co, j;
long i, Seq;
unsigned long Z;
double X2;
for (Seq = 1; Seq <= N; Seq++) {
for (i = res->jmin; i <= res->jmax; i++)
res->Count[i] = 0;
for (i = 0; i < n; i++) {
/* Generate a block of L bits */
co = 0;
for (j = 0; j < LDivs; j++) {
Z = unif01_StripB (gen, r, s);
while (Z > 0) { /* Count the number of 1 bits */
Z &= Z - 1; /* Clear lowest 1 bit */
++co;
}
}
/* The last bits of the block */
if (LMods > 0) {
Z = unif01_StripB (gen, r, LMods);
while (Z > 0) {
Z &= Z - 1;
++co;
}
}
++res->Count[res->Loc[co]];
}
X2 = gofs_Chi2 (res->NbExp, res->Count, res->jmin, res->jmax);
statcoll_AddObs (res->sVal1, X2);
if (swrite_Counters)
tables_WriteTabL (res->Count, res->jmin, res->jmax, 5, 10,
"Observed numbers of blocks");
}
}
/*-------------------------------------------------------------------------*/
static void HammingWeight_S (unif01_Gen * gen, sres_Chi2 * res,
long N, long n, int r, int s, long L)
/*
* Generate all the n*L bits for the HammingWeight test in the case L <= s.
* A number generates sDivL blocks. If s % L == 0, we use all s bits of the
* number.
*/
{
const int sDivL = s / L; /* A number generates sDivL blocks */
const int s1 = s - s % L;
const long Q = n / sDivL;
const int Q2 = n % sDivL;
const unsigned long MASK = num_TwoExp[L] - 1.0;
int co, j;
long i, Seq;
unsigned long Z, Y;
double X2;
for (Seq = 1; Seq <= N; Seq++) {
for (i = res->jmin; i <= res->jmax; i++)
res->Count[i] = 0;
for (i = 0; i < Q; i++) {
Z = unif01_StripB (gen, r, s1);
/* Generate sDivL blocks of L bits */
for (j = 0; j < sDivL; j++) {
co = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++co;
}
++res->Count[res->Loc[co]];
Z >>= L;
}
}
/* The last bits */
if (Q2 > 0) {
Z = unif01_StripB (gen, r, Q2 * L);
for (j = 0; j < Q2; j++) {
co = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++co;
}
++res->Count[res->Loc[co]];
Z >>= L;
}
}
X2 = gofs_Chi2 (res->NbExp, res->Count, res->jmin, res->jmax);
statcoll_AddObs (res->sVal1, X2);
if (swrite_Counters)
tables_WriteTabL (res->Count, res->jmin, res->jmax, 5, 10,
"Observed numbers of blocks");
}
}
/*-------------------------------------------------------------------------*/
void sstring_HammingWeight (unif01_Gen * gen, sres_Chi2 * res,
long N, long n, int r, int s, long L)
{
long i;
double V[1]; /* Number of Chi2 degrees of freedom */
char str[LEN1 + 1] = "";
fmass_INFO Q;
lebool localRes = FALSE;
chrono_Chrono *Timer;
long jlow, jhigh;
long NbGroups; /* Number of classes */
char *TestName = "sstring_HammingWeight test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataLongHead (gen, TestName, N, n, r, s, L);
util_Assert (r + s <= 32, "sstring_HammingWeight: r + s > 32");
util_Assert (L >= 2, "sstring_HammingWeight: L < 2");
if (res == NULL) {
localRes = TRUE;
res = sres_CreateChi2 ();
}
sres_InitChi2 (res, N, L, "sstring_HammingWeight");
Q = fmass_CreateBinomial (L, 0.5, 0.5);
for (i = 0; i <= L; i++)
res->NbExp[i] = n * fmass_BinomialTerm2 (Q, i);
fmass_DeleteBinomial (Q);
jlow = 0;
jhigh = L;
if (swrite_Classes)
gofs_WriteClasses (res->NbExp, res->Loc, jlow, jhigh, 0);
gofs_MergeClasses (res->NbExp, res->Loc, &jlow, &jhigh, &NbGroups);
if (swrite_Classes)
gofs_WriteClasses (res->NbExp, res->Loc, jlow, jhigh, NbGroups);
res->jmin = jlow;
res->jmax = jhigh;
res->degFree = NbGroups - 1;
if (res->degFree < 1) {
if (localRes)
sres_DeleteChi2 (res);
return;
}
sprintf (str, "The N statistic values (a ChiSquare with %1ld degrees"
" of freedom):", NbGroups - 1);
statcoll_SetDesc (res->sVal1, str);
if (L >= s)
HammingWeight_L (gen, res, N, n, r, s, L);
else
HammingWeight_S (gen, res, N, n, r, s, L);
V[0] = res->degFree;
gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_ChiSquare, V,
res->sVal2, res->pVal2);
res->pVal1->NObs = N;
sres_GetChi2SumStat (res);
if (swrite_Collectors)
statcoll_Write (res->sVal1, 5, 14, 2, 1);
if (swrite_Basic) {
swrite_AddStrChi (str, LEN1, res->degFree);
gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2, str);
swrite_Chi2SumTest (N, res);
swrite_Final (gen, Timer);
}
if (localRes)
sres_DeleteChi2 (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
#if 0
void sstring_Run0 (unif01_Gen * gen, sres_Basic * res,
long N, long n, int r, int s)
{
const long K = n / s; /* If n % s != 0, a string will contain
K * s bits instead of n */
const unsigned long SBIT = 1UL << (s - 1);
unsigned long jBit; /* Position of current bit in Z */
int pBit; /* Previous bit */
long i, Seq;
long co1; /* Counter for number of 1 */
long cor; /* Counter for number of runs */
unsigned long Z;
double X, f1;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_Run test";
Timer = chrono_Create ();
n = K * s;
if (swrite_Basic)
WriteDataPeriod (gen, TestName, N, n, r, s);
util_Assert (r + s <= 32, "sstring_Run: r + s > 32");
/* util_Assert (100 <= n, "sstring_Run: n < 100"); */
if (res == NULL) {
localRes = TRUE;
res = sres_CreateBasic ();
}
sres_InitBasic (res, N, "sstring_Run");
statcoll_SetDesc (res->sVal1, "sVal1: a standard normal");
for (Seq = 1; Seq <= N; Seq++) {
co1 = cor = 0;
/* Be sure to count the first run with pBit != {0, 1} */
pBit = 2;
for (i = 0; i < K; i++) {
Z = unif01_StripB (gen, r, s);
jBit = SBIT;
/* Add the number of 1 bit and number of runs in Z */
while (jBit > 0) {
if (Z & jBit) { /* bit 1 */
co1++;
if (pBit != 1)
cor++;
pBit = 1;
} else { /* bit 0 */
if (pBit != 0)
cor++;
pBit = 0;
}
jBit >>= 1;
}
}
f1 = (double) co1 / (K * s);
X = (cor - n * 2.0 * f1 * (1.0 - f1)) /
(2.0 * sqrt ((double) n) * f1 * (1.0 - f1));
statcoll_AddObs (res->sVal1, X);
}
gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Normal,
(double *) NULL, res->sVal2, res->pVal2);
res->pVal1->NObs = N;
if (swrite_Collectors)
statcoll_Write (res->sVal1, 5, 14, 4, 3);
if (swrite_Basic) {
gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
"Normal statistic :");
swrite_Final (gen, Timer);
}
if (localRes)
sres_DeleteBasic(res);
chrono_Delete (Timer);
}
#endif
/*=========================================================================*/
void sstring_Run (unif01_Gen * gen, sstring_Res3 *res,
long N, long n, int r, int s)
{
const unsigned long SBIT = 1UL << (s - 1);
const double sr = s;
unsigned long jBit; /* Position of current bit in Z */
int pBit; /* Previous bit */
int k, j;
long Seq;
double cob; /* Counter for number of bits */
long cor; /* Counter for number of 1 runs */
int len; /* Length of current run */
unsigned long Z;
double X2, X, temp;
char str[LEN1 + 1];
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_Run test";
sres_Basic *NBits;
sres_Chi2 *NRuns;
long *Count0, *Count1;
double *Prob, *NbExp;
double Param[1];
Timer = chrono_Create ();
k = 1 + num_Log2 (n / gofs_MinExpected);
if (swrite_Basic)
WriteDataPeriod (gen, TestName, N, n, r, s);
util_Assert (r + s <= 32, "sstring_Run: r + s > 32");
/* util_Assert (100 <= n, "sstring_Run: n < 100"); */
if (res == NULL) {
localRes = TRUE;
res = sstring_CreateRes3 ();
}
InitRes3 (res, N, k);
NBits = res->NBits;
NRuns = res->NRuns;
Count0 = res->Count0;
Count1 = res->Count1;
statcoll_SetDesc (NBits->sVal1,
"The N statistic values (a standard normal):");
sprintf (str, "The N statistic values (a ChiSquare with %1d degrees"
" of freedom):", 2*(k - 1));
statcoll_SetDesc (NRuns->sVal1, str);
Prob = util_Calloc (1 + (size_t) k, sizeof (double));
Prob[0] = 1.0;
for (j = 1; j < k; j++) {
Prob[j] = Prob[j - 1] / 2.0;
NRuns->NbExp[j] = n * Prob[j];
}
Prob[k] = Prob[k - 1];
NRuns->NbExp[k] = n * Prob[k];
util_Assert (NRuns->NbExp[k] >= gofs_MinExpected,
"sstring_Run: NRuns->NbExp[k] < gofs_MinExpected");
if (swrite_Classes)
gofs_WriteClasses (NRuns->NbExp, NRuns->Loc, 1, k, 0);
NRuns->jmax = k;
NRuns->jmin = 1;
NRuns->degFree = 2*(k - 1);
if (NRuns->degFree < 1) {
util_Warning (TRUE, "Chi-square with 0 degree of freedom.");
if (localRes)
sstring_DeleteRes3 (res);
return;
}
for (Seq = 1; Seq <= N; Seq++) {
cob = cor = len = 0;
for (j = 1; j <= k; j++) {
Count0[j] = 0;
Count1[j] = 0;
}
/* Make sure to count the first run; set pBit != {0, 1} */
pBit = 2;
while (cor < n) {
Z = unif01_StripB (gen, r, s);
jBit = SBIT;
cob += sr;
if (len >= n) {
util_Warning (TRUE, "sstring_Run: all bits are 0 !");
util_Free (Prob);
if (localRes)
sstring_DeleteRes3 (res);
return;
}
/* Add the number of runs in Z */
while (jBit > 0) {
if (Z & jBit) { /* bit 1 */
if (pBit != 1) {
cor++;
if (len < k)
(Count0[len])++;
else
(Count0[k])++;
len = 1;
} else {
len++;
}
pBit = 1;
} else { /* bit 0 */
if (pBit != 0) {
if (len < k)
(Count1[len])++;
else
(Count1[k])++;
len = 1;
} else {
len++;
}
pBit = 0;
}
jBit >>= 1;
}
}
X2 = 0.0;
NbExp = NRuns->NbExp;
for (j = NRuns->jmin; j <= NRuns->jmax; j++) {
temp = Count0[j] - NbExp[j];
X2 += temp * temp / (NbExp[j] * (1.0 - Prob[j]));
}
X = X2;
X2 = 0.0;
for (j = NRuns->jmin; j <= NRuns->jmax; j++) {
temp = Count1[j] - NbExp[j];
X2 += temp * temp / (NbExp[j] * (1.0 - Prob[j]));
}
statcoll_AddObs (NRuns->sVal1, X2 + X);
if (swrite_Counters) {
tables_WriteTabL (Count0, 1, k, 5, 10,
"Observed number of runs of 0");
tables_WriteTabL (Count1, 1, k, 5, 10,
"Observed number of runs of 1");
}
X = (cob - 4.0 * n) / sqrt (8.0 * n);
statcoll_AddObs (NBits->sVal1, X);
}
Param[0] = 2*(k - 1);
gofw_ActiveTests2 (NRuns->sVal1->V, NRuns->pVal1->V, N, wdist_ChiSquare,
Param, NRuns->sVal2, NRuns->pVal2);
NRuns->pVal1->NObs = N;
sres_GetChi2SumStat (NRuns);
gofw_ActiveTests2 (NBits->sVal1->V, NBits->pVal1->V, N, wdist_Normal,
(double *) NULL, NBits->sVal2, NBits->pVal2);
NBits->pVal1->NObs = N;
sres_GetNormalSumStat (NBits);
if (swrite_Basic) {
printf ("\n-----------------------------------------------\n");
if (N == 1) {
printf ("Total number of 1 runs: %1ld\n\n", cor);
printf ("Number of degrees of freedom : %4ld\n",
NRuns->degFree);
printf ("Chi2 statistic for number of runs :");
gofw_Writep2 (NRuns->sVal2[gofw_Mean], NRuns->pVal2[gofw_Mean]);
} else {
printf ("Test results for the number of runs:\n");
gofw_WriteActiveTests0 (N, NRuns->sVal2, NRuns->pVal2);
swrite_Chi2SumTest (N, NRuns);
}
if (swrite_Collectors)
statcoll_Write (NRuns->sVal1, 5, 14, 4, 3);
printf ("\n-----------------------------------------------\n");
if (N == 1) {
printf ("Total number of bits: %.0f\n\n", cob);
printf ("Normal statistic for number of bits :");
gofw_Writep2 (NBits->sVal2[gofw_Mean], NBits->pVal2[gofw_Mean]);
} else {
printf ("Test results for the number of bits:\n");
gofw_WriteActiveTests0 (N, NBits->sVal2, NBits->pVal2);
swrite_NormalSumTest (N, NBits);
}
if (swrite_Collectors)
statcoll_Write (NBits->sVal1, 5, 14, 4, 3);
printf ("\n\n");
swrite_Final (gen, Timer);
}
util_Free (Prob);
if (localRes)
sstring_DeleteRes3 (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static void WriteDataAutoCor (unif01_Gen *gen, char *Test,
long N, long n, int r, int s, int d)
{
swrite_Head (gen, Test, N, n, r);
printf (", s = %1d, d = %1d\n\n", s, d);
}
/*-------------------------------------------------------------------------*/
void sstring_AutoCor (unif01_Gen * gen, sres_Basic * res,
long N, long n, int r, int s, int d)
{
const long K = (n - d) / s;
const long M = d / s + 2;
unsigned long *Y; /* Circular buffer for random numbers */
unsigned long A; /* Correlation */
unsigned long Z, s1, s2;
unsigned long mask1, mask2; /* Masks of s1, s2 least sig. bits */
double X;
long i, Seq;
int j1, j2;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_AutoCor test";
Timer = chrono_Create ();
/* There are a few bits less than n */
n -= (n - d) % s;
if (swrite_Basic)
WriteDataAutoCor (gen, TestName, N, n, r, s, d);
util_Assert (r + s <= 32, "sstring_AutoCor: r + s > 32");
util_Assert (d <= n / 2, "sstring_AutoCor: d > n/2");
util_Assert (d > 0, "sstring_AutoCor: d < 1");
if (res == NULL) {
localRes = TRUE;
res = sres_CreateBasic ();
}
sres_InitBasic (res, N, "sstring_AutoCor");
Y = util_Calloc ((size_t) M, sizeof (unsigned long));
statcoll_SetDesc (res->sVal1, "sVal1: a standard normal");
s1 = d % s;
s2 = s - s1;
mask1 = num_TwoExp[s1] - 1.0;
mask2 = num_TwoExp[s2] - 1.0;
for (Seq = 1; Seq <= N; Seq++) {
/* Fill circular buffer with first random numbers */
for (i = 0; i < M-1; i++)
Y[i] = unif01_StripB (gen, r, s);
A = 0;
j1 = M - 1;
j2 = M - 2;
for (i = 0; i < K; i++) {
Y[j1] = unif01_StripB (gen, r, s);
j1 = (j1 + 1) % M;
Z = ((Y[j1] >> s1) ^ Y[j2]) & mask2;
while (Z > 0) { /* Count the number of 1 bits in Z */
Z &= Z - 1; /* Clear lowest 1 bit */
++A;
}
j2 = (j2 + 1) % M;
Z = ((Y[j2] >> s2) ^ Y[j1]) & mask1;
while (Z > 0) {
Z &= Z - 1;
++A;
}
}
X = 2.0 * (A - (n - d) / 2.0) / sqrt ((double) (n - d));
statcoll_AddObs (res->sVal1, X);
}
gofw_ActiveTests2 (res->sVal1->V, res->pVal1->V, N, wdist_Normal,
(double *) NULL, res->sVal2, res->pVal2);
res->pVal1->NObs = N;
sres_GetNormalSumStat (res);
if (swrite_Collectors)
statcoll_Write (res->sVal1, 5, 14, 4, 3);
if (swrite_Basic) {
gofw_WriteActiveTests2 (N, res->sVal2, res->pVal2,
"Normal statistic :");
swrite_NormalSumTest (N, res);
swrite_Final (gen, Timer);
}
util_Free (Y);
if (localRes)
sres_DeleteBasic (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static void WriteDataHammingCorr (unif01_Gen *gen, char *TestName,
long N, long n, int r, int s, int L)
{
swrite_Head (gen, TestName, N, n, r);
printf (", s = %1d, L = %1d\n\n\n", s, L);
}
/*-------------------------------------------------------------------------*/
static void HammingCorr_L (unif01_Gen *gen, sstring_Res * res,
long n, int r, int s, int L)
/*
* Generate all the n bits for the HammingCorr test in the case L > s.
* For the last number generated in a block of L bits, we keep its first
* LMods bits and discard the other bits.
*/
{
const int LMods = L % s;
const int LDivs = L / s;
int Pre, X;
int j;
long k;
unsigned long Z;
/* Junk value to avoid a test "if (k == 1)" for every generated
number; it will not be counted. */
Pre = L + 1;
for (k = 1; k <= n; k++) {
/* Generate a sequence of L bits */
X = 0;
for (j = 1; j <= LDivs; j++) {
Z = unif01_StripB (gen, r, s);
/* Count the number of 1 bits */
while (Z > 0) {
Z &= Z - 1; /* Clear lowest 1 bit */
++X;
}
}
/* The last bits of the sequence */
if (LMods > 0) {
Z = unif01_StripB (gen, r, LMods);
while (Z > 0) {
Z &= Z - 1;
++X;
}
}
++res->Counters[Pre][X];
Pre = X;
}
}
/*-------------------------------------------------------------------------*/
static void HammingCorr_S (unif01_Gen *gen, sstring_Res * res,
long n, int r, int s, int L)
/*
* Generate all the n bits for the HammingCorr test in the case L <= s.
* A number generates sDivL blocks. If s % L == 0, we use all s bits of the
* number.
*/
{
const int sDivL = s / L; /* A number generates sDivL blocks */
const long Q = n / sDivL;
const long Q1 = n % sDivL;
const unsigned long MASK = num_TwoExp[L] - 1.0;
int Pre, X;
int j;
long k;
unsigned long Z, Y;
/* Junk value to avoid a test "if (k == 1)" for every generated
number; it will not be counted. */
Pre = L + 1;
for (k = 0; k < Q; k++) {
Z = unif01_StripB (gen, r, s);
for (j = 0; j < sDivL; j++) {
X = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++X;
}
++res->Counters[Pre][X];
Pre = X;
Z >>= L;
}
}
/* The last Q1 blocks */
if (Q1 > 0) {
Z = unif01_StripB (gen, r, s);
for (j = 0; j < Q1; j++) {
X = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++X;
}
++res->Counters[Pre][X];
Pre = X;
Z >>= L;
}
}
}
/*-------------------------------------------------------------------------*/
void sstring_HammingCorr (unif01_Gen * gen, sstring_Res * res,
long N, long n, int r, int s, int L)
{
int i, j;
long Seq;
double Sum;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_HammingCorr test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataHammingCorr (gen, TestName, N, n, r, s, L);
util_Assert (s <= num_MaxTwoExp, "sstring_HammingCorr: s too large");
util_Assert ((unsigned) s <= CHAR_BIT * sizeof (unsigned long),
"sstring_HammingCorr: s too large");
if (res == NULL) {
localRes = TRUE;
res = sstring_CreateRes ();
}
InitRes (res, N, L, -1, "sstring_HammingCorr");
statcoll_SetDesc (res->Bas->sVal1, "HammingCorr sVal1: standard normal");
for (Seq = 1; Seq <= N; Seq++) {
for (i = 0; i <= L; i++) {
for (j = 0; j <= L; j++)
res->Counters[i][j] = 0;
}
if (L >= s)
HammingCorr_L (gen, res, n, r, s, L);
else
HammingCorr_S (gen, res, n, r, s, L);
if (swrite_Counters)
/* Print the matrix of counters */
tables_WriteMatrixL (res->Counters, 0, L, 0, L, 8,
res->Style, "Number of pairs [0..L, 0..L]");
/* Calculate statistic */
Sum = 0.0;
for (i = 0; i <= L; i++) {
for (j = 0; j <= L; j++)
Sum += res->Counters[i][j] * (i - L / 2.0) * (j - L / 2.0);
}
Sum = Sum * 4.0 / (L * sqrt (n - 1.0));
statcoll_AddObs (res->Bas->sVal1, Sum);
}
gofw_ActiveTests2 (res->Bas->sVal1->V, res->Bas->pVal1->V, N,
wdist_Normal, (double *) NULL, res->Bas->sVal2, res->Bas->pVal2);
res->Bas->pVal1->NObs = N;
sres_GetNormalSumStat (res->Bas);
if (swrite_Collectors)
statcoll_Write (res->Bas->sVal1, 5, 14, 4, 3);
if (swrite_Basic) {
gofw_WriteActiveTests2 (N, res->Bas->sVal2, res->Bas->pVal2,
"Normal statistic :");
swrite_NormalSumTest (N, res->Bas);
swrite_Final (gen, Timer);
}
if (localRes)
sstring_DeleteRes (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
static void WriteDataHammingIndep (unif01_Gen * gen, char *TestName,
long N, long n, int r, int s, int L, int d)
{
swrite_Head (gen, TestName, N, n, r);
printf (", s = %1d, L = %1d, d = %1d\n\n\n", s, L, d);
}
/*-------------------------------------------------------------------------*/
static void HammingIndep_L (unif01_Gen *gen, sstring_Res * res,
long n, int r, int s, int L)
/*
* Generate all the n bits for the HammingIndep test in the case L > s.
* For the last number generated in a block of L bits, we keep its first
* LMods bits and discard the other bits.
*/
{
int Pre; /* Previous value of X */
int X;
const int LMods = L % s;
const int LDivs = L / s;
int j;
unsigned long U;
unsigned long TwonUL;
unsigned long ic;
/* For the test with n >= 2^30 */
TwonUL = 2 * (unsigned long) n;
Pre = 0; /* Eliminate a warning from compiler */
for (ic = 1; ic <= TwonUL; ic++) {
/* Generate 1 block of L bits */
X = 0;
for (j = 1; j <= LDivs; j++) {
U = unif01_StripB (gen, r, s);
/* Count the number of 1 bits */
while (U > 0) {
U &= U - 1; /* Clear lowest 1 bit */
++X;
}
}
/* The last bits of the block */
if (LMods > 0) {
U = unif01_StripB (gen, r, LMods);
while (U > 0) {
U &= U - 1;
++X;
}
}
/* Non-overlapping pairs; count only when ic % 2 == 0 */
if (!(ic & 1))
++res->Counters[Pre][X];
Pre = X;
}
}
/*-------------------------------------------------------------------------*/
static void HammingIndep_S (unif01_Gen *gen, sstring_Res * res,
long n, int r, int s, int L)
/*
* Generate all the n bits for the HammingIndep test in the case L <= s.
* A number generates sDivL blocks. If s % L == 0, we use all s bits of the
* number.
*/
{
const int sDivL = s / L; /* A number generates sDivL blocks */
const unsigned long MASK = num_TwoExp[L] - 1.0;
int Pre; /* Previous value of X */
int X;
int j;
unsigned long Q, Q1, i;
unsigned long Z, Y;
unsigned long TwonUL;
unsigned long bloc = 0;
/* For the test with n >= 2^30 */
TwonUL = 2 * (unsigned long) n;
Q = TwonUL / sDivL;
Q1 = TwonUL % sDivL;
Pre = 0; /* Eliminate a warning from compiler */
for (i = 0; i < Q; i++) {
Z = unif01_StripB (gen, r, s);
for (j = 0; j < sDivL; j++) {
X = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++X;
}
/* Non-overlapping pairs; count only when bloc % 2 == 0 */
if (!(++bloc & 1))
++res->Counters[Pre][X];
Pre = X;
Z >>= L;
}
}
/* The last Q1 blocks */
if (Q1 > 0)
Z = unif01_StripB (gen, r, s);
for (i = 0; i < Q1; i++) {
X = 0;
Y = Z & MASK;
while (Y > 0) { /* Count the number of 1 bits */
Y &= Y - 1; /* Clear lowest 1 bit */
++X;
}
if (!(++bloc & 1))
++res->Counters[Pre][X];
Pre = X;
Z >>= L;
}
}
/*-------------------------------------------------------------------------*/
static void CountBlocks (
sstring_Res *res,
int L, /* Length of blocks (num bits) */
int d /* Rank of rows-columns counted */
)
/*
* Add the diagonal blocks for the matrix res->Counters, i.e. count the
* number of values in the 4 corners of the matrix with the center removed.
* We count only rows-columns of rank >= d starting from the center of
* the matrix. We eliminate 2d - 1 rows and columns at the center of the
* matrix when L is even. When L is odd, we eliminate 2d - 2 rows and
* columns at the center (except when d = 1, where we keep all the rows
* and columns).
*/
{
int L2, L1, k, j, i;
L2 = L1 = L / 2;
if ((L & 1))
++L1;
for (k = 1; k <= d; k++) {
res->XD[k][0] = 0;
res->XD[k][1] = 0;
/* Block ++ */
for (i = 0; i <= L1 - k; i++) {
for (j = 0; j <= L1 - k; j++) {
res->XD[k][0] += res->Counters[i][j];
}
}
/* Block -- */
for (i = L2 + k; i <= L; i++) {
for (j = L2 + k; j <= L; j++) {
res->XD[k][0] += res->Counters[i][j];
}
}
/* Block +- */
for (i = 0; i <= L1 - k; i++) {
for (j = L2 + k; j <= L; j++) {
res->XD[k][1] += res->Counters[i][j];
}
}
/* Block -+ */
for (i = L2 + k; i <= L; i++) {
for (j = 0; j <= L1 - k; j++) {
res->XD[k][1] += res->Counters[i][j];
}
}
}
}
/*-------------------------------------------------------------------------*/
static void WriteBlocs (
sstring_Res *res,
int d /* Rank of rows-columns */
)
/*
* Print the sum of counters in the diagonal blocks for different d
*/
{
int i;
printf ("--------------------------------------------------\n");
for (i = 1; i <= d; i++) {
printf ("The number of blocks ++, -- with d >= %1d is %10ld\n",
i, res->XD[i][0]);
printf ("The number of blocks +-, -+ with d >= %1d is %10ld\n\n",
i, res->XD[i][1]);
}
printf ("\n");
}
/*-------------------------------------------------------------------------*/
void sstring_HammingIndep (unif01_Gen * gen, sstring_Res * res,
long N, long n, int r, int s, int L, int d)
{
int Liber = 0; /* Num of degrees of freedom for main test */
int Liberte; /* Num of degrees of freedom for block tests */
int X;
int i, j;
long Seq;
double NbEsp;
double Var;
double NbMoyen;
double X2;
double Sum;
const double nLR = n;
double Z;
double *Prob;
double *NbEsp5;
long *Count5;
double V[1];
fmass_INFO Q;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sstring_HammingIndep test";
char chaine[LEN1 + 1] = "";
char str[LEN2 + 1];
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataHammingIndep (gen, TestName, N, n, r, s, L, d);
if (n < 2.0 * gofs_MinExpected) {
util_Warning (TRUE, "sstring_HammingIndep: n < 20");
return;
}
/* util_Assert (n >= 30, "sstring_HammingIndep: n < 30"); */
util_Assert (d <= sstring_MAXD, "sstring_HammingIndep: d > sstring_MAXD");
util_Assert (((L + 1) / 2) >= d, "sstring_HammingIndep: d > (L + 1) / 2");
util_Assert (s <= num_MaxTwoExp, "sstring_HammingIndep: s too large");
util_Assert ((unsigned) s <= CHAR_BIT * sizeof (unsigned long),
"sstring_HammingIndep: s too large");
if (res == NULL) {
localRes = TRUE;
res = sstring_CreateRes ();
}
InitRes (res, N, L, d, "sstring_HammingIndep");
for (i = 1; i <= d; i++) {
strncpy (chaine, "HammingIndep Block[", (size_t) LEN1);
sprintf (str, "%1d", i);
strncat (chaine, str, (size_t) LEN2);
strncat (chaine, "]", (size_t) 2);
statcoll_SetDesc (res->Block[i]->sVal1, chaine);
}
strncpy (chaine, "\nCounters with expected numbers >= ", (size_t) LEN1);
sprintf (str, "%g", gofs_MinExpected);
strncat (chaine, str, (size_t) LEN2);
statcoll_SetDesc (res->Bas->sVal1, chaine);
Prob = util_Calloc ((size_t) L + 1, sizeof (double));
NbEsp5 = util_Calloc ((size_t) (L + 1)*(L + 1) + 1, sizeof (double));
Count5 = util_Calloc ((size_t) (L + 1)*(L + 1) + 1, sizeof (long));
Q = fmass_CreateBinomial (L, 0.5, 0.5);
for (i = 0; i <= L; i++)
Prob[i] = fmass_BinomialTerm2 (Q, i);
fmass_DeleteBinomial (Q);
for (Seq = 1; Seq <= N; Seq++) {
for (i = 0; i <= L; i++) {
for (j = 0; j <= L; j++) {
res->Counters[i][j] = 0;
}
}
if (L >= s)
HammingIndep_L (gen, res, n, r, s, L);
else
HammingIndep_S (gen, res, n, r, s, L);
/* The cells for which the expected number >= gofs_MinExpected will */
/* correspond to one class each. Merge all other cells such that */
/* their expected number is < gofs_MinExpected into 1 big class: Z */
/* will contain the sum of all their expected numbers, and the */
/* counter X will contain the sum of all their observed numbers. We */
/* shall apply a chi-square test with Liber degrees of freedom. We */
/* have (Liber + 1) classes. */
Liber = 0;
Z = 0.0;
X = 0;
for (i = 0; i <= L; i++) {
for (j = 0; j <= L; j++) {
res->ZCounters[i][j] = nLR * Prob[i] * Prob[j];
if (res->ZCounters[i][j] >= gofs_MinExpected) {
NbEsp5[Liber] = res->ZCounters[i][j];
Count5[Liber] = res->Counters[i][j];
++Liber;
} else {
Z += res->ZCounters[i][j];
X += res->Counters[i][j];
}
}
}
if (Z >= gofs_MinExpected) {
/* We have one more class */
/* if (swrite_Classes && Seq == 1) {
printf ("All cells with NbExp < ");
printf ("%s", str);
printf (" are merged into one big class with NbExp = %f\n\n", Z);
} */
NbEsp5[Liber] = Z;
Count5[Liber] = X;
} else if (Liber > 0) {
/* We add them to the last class instead */
--Liber;
NbEsp5[Liber] += Z;
Count5[Liber] += X;
}
/* Everything has been put in a single class; separate all in two
classes */
if (Liber == 0) {
Z = 0.0;
X = 0;
for (i = 0; i <= L; i++) {
for (j = 0; j <= L / 2; j++) {
res->ZCounters[i][j] = nLR * Prob[i] * Prob[j];
Z += res->ZCounters[i][j];
X += res->Counters[i][j];
}
}
NbEsp5[0] = Z;
Count5[0] = X;
Z = 0.0;
X = 0;
for (i = 0; i <= L; i++) {
for (j = 1 + L / 2; j <= L; j++) {
res->ZCounters[i][j] = nLR * Prob[i] * Prob[j];
Z += res->ZCounters[i][j];
X += res->Counters[i][j];
}
}
NbEsp5[1] = Z;
Count5[1] = X;
Liber = 1;
}
if (Liber > 0) {
X2 = gofs_Chi2 (NbEsp5, Count5, 0, Liber);
statcoll_AddObs (res->Bas->sVal1, X2);
}
/* Compute the normalized observed number for each pair */
for (i = 0; i <= L; i++) {
for (j = 0; j <= L; j++) {
NbEsp = nLR * Prob[i] * Prob[j];
Var = NbEsp * (1.0 - Prob[i] * Prob[j]);
if (Var <= 0.0) {
/* This case will occur when NbEsp = 0; if so, fail the test */
/* when res->Counters[i][j] != 0 */
res->ZCounters[i][j] = (res->Counters[i][j] - NbEsp) * 1.E100;
} else {
res->ZCounters[i][j] = (res->Counters[i][j] - NbEsp) /
sqrt (Var);
}
}
}
if (sstring_Counters)
/* Print the matrix of counters */
tables_WriteMatrixL (res->Counters, 0, L, 0, L, 8, res->Style,
"res->Counters, the number of pairs [0..L, 0..L]");
if (swrite_Counters)
/* Print the matrix of normalized counters */
tables_WriteMatrixD (res->ZCounters, 0, L, 0, L, 12, 4,
res->Style, "res->ZCounters, the normalized counters");
/* These blocks are sub-matrices symmetrically placed with respect to */
/* the diagonals in the matrix of the number of pairs [Xi, X(i+1)]. */
/* For those, we shall apply a chi-square test for the total number */
/* in the diagonal sub-matrices. */
CountBlocks (res, L, d);
if (swrite_Counters)
WriteBlocs (res, d);
for (i = 1; i <= d; i++) {
double NumExp[3];
long Count[3];
Sum = 0.0;
for (j = 0; j <= (L + 1) / 2 - i; j++)
Sum += Prob[j];
/* Probability of a sub-matrix block */
Sum *= Sum;
/* Average number for each of the 2 blocks */
NbMoyen = Sum * nLR * 2.0;
NumExp[0] = NumExp[1] = NbMoyen;
if (2.0*NbMoyen < gofs_MinExpected)
printf ("******* sample too small for chi-square for d = %d\n", i);
NumExp[2] = nLR - 2.0*NbMoyen;
Count[0] = res->XD[i][0];
Count[1] = res->XD[i][1];
Count[2] = n - Count[0] - Count[1];
X2 = gofs_Chi2 (NumExp, Count, 0, 2);
statcoll_AddObs (res->Block[i]->sVal1, X2);
}
}
for (i = 1; i <= d; i++) {
/* Degrees of freedom */
if ((L & 1) && i == 1)
Liberte = 1;
else
Liberte = 2;
V[0] = Liberte;
gofw_ActiveTests2 (res->Block[i]->sVal1->V, res->Block[i]->pVal1->V, N,
wdist_ChiSquare, V, res->Block[i]->sVal2, res->Block[i]->pVal2);
res->Block[i]->pVal1->NObs = N;
Sum = N * statcoll_Average (res->Block[i]->sVal1);
res->Block[i]->sVal2[gofw_Sum] = Sum;
res->Block[i]->pVal2[gofw_Sum] = fbar_ChiSquare2 (N*Liberte, 12, Sum);
if (swrite_Basic) {
printf ("\nDiagonal blocks with d = %2d", i);
swrite_AddStrChi (str, LEN2, Liberte);
gofw_WriteActiveTests2 (N, res->Block[i]->sVal2,
res->Block[i]->pVal2, str);
swrite_Chi2SumTestb (N, res->Block[i]->sVal2[gofw_Sum],
res->Block[i]->pVal2[gofw_Sum], Liberte);
if (swrite_Collectors) {
strncpy (chaine, res->Block[i]->sVal1->Desc, (size_t) LEN1);
strncat (chaine, ": a chi2 with ", (size_t) LEN1);
sprintf (str, "%1d", Liberte);
strncat (chaine, str, (size_t) LEN2);
strncat (chaine, " degrees of freedom", (size_t) LEN1);
statcoll_SetDesc (res->Block[i]->sVal1, chaine);
statcoll_Write (res->Block[i]->sVal1, 5, 14, 4, 3);
}
}
}
if (Liber > 0) {
V[0] = Liber;
gofw_ActiveTests2 (res->Bas->sVal1->V, res->Bas->pVal1->V, N,
wdist_ChiSquare, V, res->Bas->sVal2, res->Bas->pVal2);
res->Bas->pVal1->NObs = N;
Sum = N * statcoll_Average (res->Bas->sVal1);
res->Bas->sVal2[gofw_Sum] = Sum;
res->Bas->pVal2[gofw_Sum] = fbar_ChiSquare2 (N*Liber, 12, Sum);
if (swrite_Basic) {
printf ("%s", res->Bas->sVal1->Desc);
swrite_AddStrChi (str, LEN2, Liber);
gofw_WriteActiveTests2 (N, res->Bas->sVal2, res->Bas->pVal2, str);
swrite_Chi2SumTestb (N, res->Bas->sVal2[gofw_Sum],
res->Bas->pVal2[gofw_Sum], Liber);
if (swrite_Collectors) {
strncpy (chaine, res->Bas->sVal1->Desc, (size_t) LEN1);
strncat (chaine, ": a ChiSquare with ", (size_t) LEN1);
sprintf (str, "%1d", Liber);
strncat (chaine, str, (size_t) LEN2);
strncat (chaine, " degrees of freedom", (size_t) LEN1);
statcoll_SetDesc (res->Bas->sVal1, chaine);
statcoll_Write (res->Bas->sVal1, 5, 14, 4, 3);
}
}
} else {
/* for the module tvaria */
res->Bas->pVal2[gofw_Mean] = -1.0;
if (d < 1)
util_Warning (1,
"n is too small: ChiSquare with 0 degree of freedom");
}
if (swrite_Basic)
swrite_Final (gen, Timer);
util_Free (Prob);
util_Free (NbEsp5);
util_Free (Count5);
if (localRes)
sstring_DeleteRes (res);
chrono_Delete (Timer);
}