htestu-0.1.0.0: cbits/testu/src/sspectral.c
/*************************************************************************\
*
* Package: TestU01
* File: sspectral.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 "num.h"
#include "chrono.h"
#include "tables.h"
#include "sspectral.h"
#include "swrite.h"
#include "wdist.h"
#include "unif01.h"
#include "gofw.h"
#include "statcoll.h"
#include "fftc.c"
#include <math.h>
#include <stdlib.h>
#include <string.h>
/*-------------------------------- Functions ------------------------------*/
static void InitRes (
sspectral_Res *res,
long N,
long jmin,
long jmax,
char *nam
)
/*
* Initializes the sspectral_Res structure
*/
{
long j;
sres_InitBasic (res->Bas, N, nam);
if (jmax > res->jmax)
res->Coef = util_Realloc (res->Coef, (jmax + 200) * sizeof (double));
for (j = 0; j <= jmax; j++)
res->Coef[j] = 0.0;
res->jmin = jmin;
res->jmax = jmax;
res->Bas->name = util_Realloc (res->Bas->name,
1 + strlen (nam) * sizeof (char));
strcpy (res->Bas->name, nam);
}
/*-------------------------------------------------------------------------*/
sspectral_Res * sspectral_CreateRes (void)
{
sspectral_Res *res;
res = util_Malloc (sizeof (sspectral_Res));
res->Bas = sres_CreateBasic ();
res->Coef = util_Calloc (1, sizeof (double));
res->jmax = 0;
return res;
}
/*-------------------------------------------------------------------------*/
void sspectral_DeleteRes (sspectral_Res *res)
{
if (res == NULL)
return;
sres_DeleteBasic (res->Bas);
util_Free (res->Coef);
util_Free (res);
}
/*=========================================================================*/
static void WriteDataFour (
unif01_Gen *gen, /* generator */
char *Test, /* Test name */
long N, /* Number of replications */
int k, /* Sample size n = 2^k */
int r, /* r first bits of each random number dropped */
int s /* s bits of each random number used */
)
{
long n;
n = num_TwoExp[k];
swrite_Head (gen, Test, N, n, r);
printf (", s = %4d, k = %4d\n\n", s, k);
}
/*-------------------------------------------------------------------------*/
void sspectral_Fourier1 (unif01_Gen *gen, sspectral_Res *res,
long N, int t, int r, int s)
{
const unsigned long SBIT = 1UL << (s - 1);
unsigned long jBit;
unsigned long Z;
long k, KALL, Seq, n, i;
double x, NbExp, h, per;
long co;
double *A;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sspectral_Fourier1 test";
Timer = chrono_Create ();
util_Assert (t <= 20, "sspectral_Fourier1: k > 20");
util_Assert (t > 1, "sspectral_Fourier1: k < 2");
if (swrite_Basic)
WriteDataFour (gen, TestName, N, t, r, s);
if (res == NULL) {
localRes = TRUE;
res = sspectral_CreateRes ();
}
n = num_TwoExp[t];
KALL = n / s;
if (n % s > 0)
KALL++;
per = 0.95;
NbExp = per * (n / 2 + 1);
/* h = 3.0 * n; */
h = 2.995732274 * n;
InitRes (res, N, 0, n, "sspectral_Fourier1");
statcoll_SetDesc (res->Bas->sVal1, "sVal1: a standard normal");
A = res->Coef;
for (Seq = 1; Seq <= N; Seq++) {
/* Fill array A: 1 for bit 1, -1 for bit 0 */
i = 0;
for (k = 0; k < KALL; k++) {
Z = unif01_StripB (gen, r, s);
jBit = SBIT;
while (jBit) {
if (jBit & Z)
A[i] = 1.0;
else
A[i] = -1.0;
jBit >>= 1;
i++;
}
}
/*
* Compute the Fourier transform of A and return the result in A. The
* first half of the array, (from 0 to n/2) is filled with the real
* components of the FFT. The second half of the array (from n/2+1 to
* n-1) is filled with the imaginary components of the FFT.
* The n new elements of A are thus:
* [Re(0), Re(1), ...., Re(n/2), Im(n/2-1), ..., Im(1)]
* The procedure is due to H.V. Sorensen, University of Pennsylvania
* and is found in file fftc.c.
*/
rsrfft (A, t);
/* Count the number of Fourier coefficients smaller than h */
co = 0;
for (i = 1; i < n / 2; i++) {
x = A[i] * A[i] + A[n - i] * A[n - i];
if (x < h)
co++;
}
if (A[0] * A[0] < h)
co++;
/* Compute the NIST statistic */
x = (co - NbExp) / sqrt (NbExp * (1.0 - per));
statcoll_AddObs (res->Bas->sVal1, x);
if (swrite_Counters) {
tables_WriteTabD (res->Coef, 0, n - 1, 5, 14, 5, 5,
"Fourier coefficients");
}
}
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_Basic) {
gofw_WriteActiveTests2 (N, res->Bas->sVal2, res->Bas->pVal2,
"Normal statistic :");
swrite_NormalSumTest (N, res->Bas);
if (swrite_Collectors)
statcoll_Write (res->Bas->sVal1, 5, 14, 4, 3);
swrite_Final (gen, Timer);
}
if (localRes)
sspectral_DeleteRes (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
void sspectral_Fourier3 (unif01_Gen *gen, sspectral_Res *res,
long N, int t, int r, int s)
{
const unsigned long SBIT = 1UL << (s - 1);
unsigned long jBit;
unsigned long Z;
long k, KALL, Seq, n, i;
double *A, *B;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sspectral_Fourier3 test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataFour (gen, TestName, N, t, r, s);
util_Assert (r + s <= 32, "sspectral_Fourier3: r + s > 32");
util_Assert (t <= 26, "sspectral_Fourier3: k > 26");
util_Assert (t >= 2, "sspectral_Fourier3: k < 2");
if (res == NULL) {
localRes = TRUE;
res = sspectral_CreateRes ();
}
n = num_TwoExp[t];
KALL = n / s + 1;
InitRes (res, n/4 + 1, 0, n, "sspectral_Fourier3");
statcoll_SetDesc (res->Bas->sVal1, "sVal1: a standard normal");
B = res->Bas->sVal1->V;
A = res->Coef;
for (i = 0; i <= n / 4; i++)
B[i] = 0.0;
for (Seq = 1; Seq <= N; Seq++) {
/* Fill array A: 1 for bit 1, -1 for bit 0 */
i = 0;
for (k = 0; k < KALL; k++) {
Z = unif01_StripB (gen, r, s);
jBit = SBIT;
while (jBit) {
if (jBit & Z)
A[i] = 1.0;
else
A[i] = -1.0;
jBit >>= 1;
i++;
}
}
/*
* Compute the Fourier transform of A and return the result in A. The
* first half of the array, (from 0 to n/2) is filled with the real
* components of the FFT. The second half of the array (from n/2+1 to
* n-1) is filled with the imaginary components of the FFT.
* The n new elements of A are thus:
* [Re(0), Re(1), ...., Re(n/2), Im(n/2-1), ..., Im(1)]
* The procedure is due to H.V. Sorensen, University of Pennsylvania
* and is found in file fftc.c.
*/
rsrfft (A, t);
/* Add the squares of the Fourier coefficients over the N replications
for each i = [1, ..., n/4], and keep them in B[i] */
for (i = 1; i <= n / 4; i++)
B[i] += A[i] * A[i] + A[n - i] * A[n - i];
if (0 && swrite_Counters)
tables_WriteTabD (B, 1, n / 4, 5, 14, 5, 5,
"Sums of square of Fourier coefficients");
}
/* There is an extra sqrt (n) factor between the Fourier coefficients
of Sorensen and those of Erdmann */
for (i = 1; i <= n / 4; i++)
B[i] /= n;
/* The N random variables have been added for each i and kept in B[i].
Their mean (1) and variance (~1) is known from Diane Erdmann. Now
consider the B[i] as n/4 normal random variables. */
for (i = 1; i <= n / 4; i++) {
B[i] = (B[i] - N) / sqrt (N * (1.0 - 2.0 / n));
statcoll_AddObs (res->Bas->sVal1, B[i]);
}
gofw_ActiveTests2 (res->Bas->sVal1->V, res->Bas->pVal1->V, n/4, wdist_Normal,
(double *) NULL, res->Bas->sVal2, res->Bas->pVal2);
res->Bas->pVal1->NObs = n/4;
if (swrite_Basic) {
gofw_WriteActiveTests2 (n/4, res->Bas->sVal2, res->Bas->pVal2,
"Normal statistic :");
if (swrite_Collectors)
statcoll_Write (res->Bas->sVal1, 5, 14, 4, 3);
swrite_Final (gen, Timer);
}
if (localRes)
sspectral_DeleteRes (res);
chrono_Delete (Timer);
}
/*=========================================================================*/
void sspectral_Fourier2 (unif01_Gen *gen, sspectral_Res *res,
long N, int t, int r, int s)
{
const unsigned long SBIT = 1UL << (s - 1);
unsigned long jBit;
unsigned long Z;
long k, KALL, Seq, n, i;
double *A;
double x, sum;
lebool localRes = FALSE;
chrono_Chrono *Timer;
char *TestName = "sspectral_Fourier2 test";
Timer = chrono_Create ();
if (swrite_Basic)
WriteDataFour (gen, TestName, N, t, r, s);
util_Assert (r + s <= 32, "sspectral_Fourier2: r + s > 32");
util_Assert (t <= 26, "sspectral_Fourier2: k > 26");
util_Assert (t >= 2, "sspectral_Fourier2: k < 2");
if (res == NULL) {
localRes = TRUE;
res = sspectral_CreateRes ();
}
n = num_TwoExp[t];
KALL = n / s + 1;
InitRes (res, N, 0, n, "sspectral_Fourier2");
statcoll_SetDesc (res->Bas->sVal1, "sVal1: a standard normal");
A = res->Coef;
for (Seq = 1; Seq <= N; Seq++) {
/* Fill array A: 1 for bit 1, -1 for bit 0 */
i = 0;
for (k = 0; k < KALL; k++) {
Z = unif01_StripB (gen, r, s);
jBit = SBIT;
while (jBit) {
if (jBit & Z)
A[i] = 1.0;
else
A[i] = -1.0;
jBit >>= 1;
i++;
}
}
/*
* Compute the Fourier transform of A and return the result in A. The
* first half of the array, (from 0 to n/2) is filled with the real
* components of the FFT. The second half of the array (from n/2+1 to
* n-1) is filled with the imaginary components of the FFT.
* The n new elements of A are thus:
* [Re(0), Re(1), ...., Re(n/2), Im(n/2-1), ..., Im(1)]
* The procedure is due to H.V. Sorensen, University of Pennsylvania
* and is found in file fftc.c.
*/
rsrfft (A, t);
/* Sum the square of the Fourier coefficients (only half of them) */
sum = 0.0;
for (i = 1; i <= n / 4; i++)
sum += A[i] * A[i] + A[n - i] * A[n - i];
/* There is an extra sqrt (n) factor between the Fourier coefficients
of Sorensen and those of Erdmann */
sum /= n;
/* Standardize the statistic */
x = 2.0*(sum - n / 4.0) / sqrt (n - 2.0);
statcoll_AddObs (res->Bas->sVal1, x);
if (swrite_Counters) {
tables_WriteTabD (res->Coef, 0, n - 1, 5, 14, 5, 5,
"Fourier coefficients");
}
}
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_Basic) {
gofw_WriteActiveTests2 (N, res->Bas->sVal2, res->Bas->pVal2,
"Normal statistic :");
swrite_NormalSumTest (N, res->Bas);
if (swrite_Collectors)
statcoll_Write (res->Bas->sVal1, 5, 14, 4, 3);
swrite_Final (gen, Timer);
}
if (localRes)
sspectral_DeleteRes (res);
chrono_Delete (Timer);
}