commit 0fe369bc5ae5dcf1bc9f6e1739fe109cbc0bcb4c Author: randogoth Date: Sat Mar 2 17:58:28 2024 +0200 removed QWQNG, added stdin diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100755 index 0000000..3f740bf --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,36 @@ +option(QNGMETER "Build qngMeter program" ON) + +if (QNGMETER) + + FIND_PACKAGE( OpenMP REQUIRED) + if(OPENMP_FOUND) + message("OPENMP FOUND") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + endif() + + ADD_DEFINITIONS("-std=c++0x") + + # Includes + include_directories( ${CMAKE_CURRENT_BINARY_DIR} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + + message(STATUS "Building qngMeter program.") + + # Source includes + include_directories(${CMAKE_SOURCE_DIR}/src) + include_directories(${CMAKE_CURRENT_SOURCE_DIR}/RandTest) + + # Target + add_executable(qngmeter QNGmeter.cpp RandTest/BiasAndAC.cpp RandTest/Entropy.cpp RandTest/Serial.cpp RandTest/Monkey.cpp + RandTest/MonkeyBitmap.cpp RandTest/Bias.cpp RandTest/AutoCorrelation.cpp RandTest/Gamma.cpp + RandTest/KolmogorovSmirnov.cpp RandTest/Stat.cpp RandTest/BitCount.cpp) + + # Dependencies + target_link_libraries(qngmeter) + +else(QNGMETER) + message(STATUS "Not building QngMeter.") +endif(QNGMETER) diff --git a/QNGmeter.cpp b/QNGmeter.cpp new file mode 100755 index 0000000..57350b1 --- /dev/null +++ b/QNGmeter.cpp @@ -0,0 +1,388 @@ +#ifdef _WIN32 +#include +#include +#endif + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "BiasAndAC.h" +#include "Monkey.h" +#include "Serial.h" +#include "Entropy.h" +#include "Stat.h" +#include "Gamma.h" +#include "KolmogorovSmirnov.h" + +int deviceType_; +int rxBytes; +int blockIntCount; +bool doExit = false; + +using namespace std; + +#ifdef __linux + #include + // Terminal Signal Handler + void sigintevent(int) + { + doExit = true; + } +#endif + +#define _BSD_SOURCE +#include + +int main() +{ + // try this to resize terminal +// resizeterm(42, 80); + + #ifdef _WIN32 + // hide cursor + CONSOLE_CURSOR_INFO info; + info.dwSize = 100; + info.bVisible = FALSE; + SetConsoleCursorInfo(GetStdHandle(STD_OUTPUT_HANDLE), &info); + // clear screen + system("cls"); + #elif __linux + // clear screen + cout << "\E[?25l\E[H\E[2J"; + // Signal handler to catch CTRL-C in terminal + signal(SIGINT, (__sighandler_t)&sigintevent); + #elif MACOSX + #endif + + // randtest classes + CBiasAndAC biasAndAc; + biasAndAc.ResetAll(); + CMonkey oqso; + CSerial serial; + CEntropy entropy; + CGamma gamma; + CKolmogorovSmirnov ks; + + double serialP = 0.5; + double serialZ = 0.0; + double KSP; + double KSN; + + time_t startTime; + time(&startTime); + string sTime = ctime(&startTime); + char sStartTime[16]; + sTime.copy(sStartTime, 15, 4); + + time_t timeNow; + time_t timePrev; + time(&timePrev); + timePrev -= 8; + + queue > > dataQueue; + + double bitsThroughputCount = 0; + double prevBitsThroughputCount = 0; + double bitsTestedCount = 0; + double bitsTestedRatio = 1; + double prevBitsTestedCount = 0; + + #ifdef _WIN32 + LARGE_INTEGER countFreq; + QueryPerformanceFrequency(&countFreq); + LARGE_INTEGER prevCount; + QueryPerformanceCounter(&prevCount); + LARGE_INTEGER nowCount; + #elif __linux + timeval start, stop; + double newTimeInterval; + gettimeofday(&start, NULL); + #elif MACOSX + #endif + + double throughput = 0; + + vector metaPs; + vector meterZs; + vector meterFlags(36, 0); + double meterScore = 0; + bool meterFreeze = false; + + omp_set_nested(true); + + #pragma omp parallel sections num_threads(2) + { + + #pragma omp section + { + // This section reads from stdin and enqueues data for testing + while (!doExit) { + shared_ptr> newBuffer(new vector()); + uint32_t input; + // Assuming the input stream provides data in a format that can be directly read into uint32_t + while (cin.read(reinterpret_cast(&input), sizeof(uint32_t))) { + newBuffer->push_back(input); + if (newBuffer->size() == 2048) { // Arbitrary buffer size, adjust based on your needs + #pragma omp critical + dataQueue.push(newBuffer); + newBuffer = shared_ptr>(new vector()); + } + } + if (!newBuffer->empty()) { + #pragma omp critical + dataQueue.push(newBuffer); + } + doExit = true; // Stop if stdin closes or reaches EOF + } + } + + + #pragma omp section + { + // This section processes the enqueued data + while (!doExit || !dataQueue.empty()) { + shared_ptr> testBuffer = nullptr; + #pragma omp critical + { + if (!dataQueue.empty()) { + testBuffer = dataQueue.front(); + dataQueue.pop(); + } + } + + if (testBuffer != nullptr) { + // Parallel processing of the testBuffer + #pragma omp parallel sections num_threads(4) + { + #pragma omp section + { for (uint32_t value : *testBuffer) biasAndAc.InsertWord32(value); } + #pragma omp section + { for (uint32_t value : *testBuffer) oqso.InsertWord32(value); } + #pragma omp section + { for (uint32_t value : *testBuffer) serial.InsertWord32(value); } + #pragma omp section + { for (uint32_t value : *testBuffer) entropy.InsertWord32(value); } + } + + // Display results + time(&timeNow); + if (difftime(timeNow, timePrev) >= 10) + { + // calc rates + #ifdef _WIN32 + QueryPerformanceCounter(&nowCount); + double newTimeInterval = ((double)nowCount.QuadPart - prevCount.QuadPart) / countFreq.QuadPart; + prevCount = nowCount; + #elif __linux + gettimeofday(&stop, NULL); + newTimeInterval = (stop.tv_sec - start.tv_sec); // sec + newTimeInterval += (stop.tv_usec - start.tv_usec) /1000000.0; // us to sec + start = stop; + #elif MACOSX + #endif + + double newBitInterval; + double newRate; + #pragma omp critical + { + newBitInterval = (bitsThroughputCount - prevBitsThroughputCount); + prevBitsThroughputCount = bitsThroughputCount; + } + + newRate = newBitInterval / newTimeInterval; + if (throughput == 0) + throughput = newRate; + else + throughput = (2*throughput + newRate) / 3; + + double newBitsTestedRatio = (bitsTestedCount-prevBitsTestedCount) / newBitInterval; + prevBitsTestedCount = bitsTestedCount; + bitsTestedRatio = (2*bitsTestedRatio + newBitsTestedRatio) / 3; + if (bitsTestedRatio > 1) + bitsTestedRatio = 1; + + // meta test and meter + metaPs.clear(); + meterZs.clear(); + + if (bitsTestedCount >= 65536) + { + // Autocorrelation KS test + for (int i=0; i<32; i++) + { + metaPs.push_back(biasAndAc.AC.P_Chi2[i]); + meterZs.push_back(biasAndAc.AC.cumulativeACZScore[i]); + } + double AcKSP; + double AcKSN; + ks.KSUP(&AcKSP, &AcKSN, &metaPs[0], metaPs.size()); + + // Combined KS test + metaPs.push_back(AcKSP); + + metaPs.push_back(biasAndAc.Bias.P_Chi2); + meterZs.push_back(biasAndAc.Bias.cumulativeBiasZScore); + } + + if (bitsTestedCount >= 4194304) + { + metaPs.push_back(serial.P_Chi2); + serialP = gamma.Gamma(128., serial.cumulativeSerialChi2); + serialZ = ks.PtoZ(serialP); + meterZs.push_back(serialZ); + + metaPs.push_back(entropy.P_Chi2); + meterZs.push_back(entropy.cumulativeZScore); + } + + if (bitsTestedCount >= 10485775) + { + metaPs.push_back(oqso.P_Chi2); + meterZs.push_back(oqso.cumulativeZScore); + } + + if (bitsTestedCount >= 65536) + { + // This KS is combined AC KSP plus with other tests + ks.KSUP(&KSP, &KSN, &metaPs[32], metaPs.size()-32); + + meterFreeze = false; + for (int i=0; i4.264897 || (metaPs[i]<0.00001 || metaPs[i]>0.99999)) + meterFlags[i] = -1; + // unfreeze condition + if (meterFlags[i] == -1) + { + if (fabs(meterZs[i])<2.326348 && (metaPs[i]>0.01 && metaPs[i]<0.99)) + meterFlags[i] = 0; + } + + if (meterFlags[i] == -1) + meterFreeze = true; + } + + // meter calc + if (meterFreeze == false) + meterScore = log(bitsTestedCount)/log(2.); + } + + + cout << endl; + cout << " QNGmeter Console 1.0 Test Type z-score p[z<=x] p[chi2<=x] " << endl; + cout << " +---------------------------+------------------------------------------------+" << endl; + cout << " | | 1/0 Balance " << setiosflags(ios::fixed) << setprecision(3) << showpos << biasAndAc.Bias.cumulativeBiasZScore << " " << setprecision(4) << noshowpos << CStat::ZtoP(biasAndAc.Bias.cumulativeBiasZScore) << " " << biasAndAc.Bias.P_Chi2 << " |" << endl; + cout << " | | Serial Test " << setiosflags(ios::fixed) << setprecision(3) << showpos << serialZ << " " << setprecision(4) << noshowpos << serialP << " " << serial.P_Chi2 << " |" << endl; + cout << " | | OQSO Test " << setiosflags(ios::fixed) << setprecision(3) << showpos << oqso.cumulativeZScore << " " << setprecision(4) << noshowpos << CStat::ZtoP(oqso.cumulativeZScore) << " " << oqso.P_Chi2 << " |" << endl; + cout << " | | Entropy Test " << setiosflags(ios::fixed) << setprecision(3) << showpos << entropy.cumulativeZScore << " " << setprecision(4) << noshowpos << CStat::ZtoP(entropy.cumulativeZScore) << " " << serial.P_Chi2 << " |" << endl; + cout << " | | H: " << setiosflags(ios::fixed) << setprecision(9) << entropy.E << " |" << endl; + + cout << " | | |" << endl; + + for (int i=1; i<=32; i++) + { + switch(i) + { + case 2: + cout << " | Start Time |"; + break; + case 3: + cout << " | " << sStartTime << " |"; + break; + case 5: + cout << " | Total Bits Tested |"; + break; + case 6: + cout << " | " << scientific << setw(9) << setprecision(2) << bitsTestedCount << fixed << " |"; + break; + case 8: + cout << " | Throughput |"; + break; + case 9: + cout << " | " << setiosflags(ios::fixed) << setw(4) << setprecision(1) << (double)(throughput/1000000.0) << " Mbps |"; + break; + case 11: + cout << " | Bits Tested Percent |"; + break; + case 12: + cout << " | " << setiosflags(ios::fixed) << setw(5) << setprecision(1) << (100*bitsTestedRatio) << "% |"; + break; + case 18: + cout << " | Meta KS+ Test |"; + break; + case 19: + cout << " | " << setiosflags(ios::fixed) << setw(5) << setprecision(3) << KSP << " |"; + break; + case 22: + cout << " | Meta KS- Test |"; + break; + case 23: + cout << " | " << setiosflags(ios::fixed) << setw(5) << setprecision(3) << KSN << " |"; + break; + case 29: + cout << " | QNGmeter Score |"; + break; + case 30: + cout << " | " << setiosflags(ios::fixed) << setw(4) << setprecision(1) << abs(meterScore) << ((meterScore<0)? "-" : (meterFreeze==false)? "+" : " ") << " |"; + break; + default: + cout << " | |"; + } + cout << " " << setw(2) << i << "st AutoCorr " << setiosflags(ios::fixed) << setprecision(3) << showpos << biasAndAc.AC.cumulativeACZScore[i-1] << " " << setprecision(4) << noshowpos << CStat::ZtoP(biasAndAc.AC.cumulativeACZScore[i-1]) << " " << biasAndAc.AC.P_Chi2[i-1] << " |" << endl; + } + cout << " +---------------------------+------------------------------------------------+" << endl; + + #ifdef _WIN32 + // put cursor in top corner + COORD coord; + coord.X = 0; + coord.Y = 0; + SetConsoleCursorPosition(GetStdHandle(STD_OUTPUT_HANDLE), coord); + #elif __linux + // clear screen + printf("\E[H"); + #elif MACOSX + #endif + + timePrev = timeNow; + } + + // End on an 'x' keypress + #ifdef _WIN32 + if (kbhit()) + { + char c = getch_(); + + if ( tolower(c) == 'x' ) + { + doExit = true; + break; + } + } + #elif __linux + // CRTL-C + #elif MACOSX + #endif + } + else + { + // give this thread a break from tight loop - waiting for data + usleep(1000); + } + } + } + } +} + diff --git a/RandTest/ACBinomialChi2.cpp b/RandTest/ACBinomialChi2.cpp new file mode 100755 index 0000000..f54ff71 --- /dev/null +++ b/RandTest/ACBinomialChi2.cpp @@ -0,0 +1,106 @@ +//#include "StdAfx.h" +#include "ACBinomialChi2.h" +#include "math.h" + +CACBinomialChi2::CACBinomialChi2(void) +{ + InitPTable(); +} + +CACBinomialChi2::~CACBinomialChi2(void) +{ +} + +// Bins incomming values for future Chi2 test +void CACBinomialChi2::BinIt(double BinVal) +{ + int BinIndex = -1; + if (BinVal <= 16090) + BinIndex = 0; + else if (BinVal >= 16679) + BinIndex = 50; + else { + BinIndex = (int)BinVal - 16091; + BinIndex /= 12; + BinIndex += 1; + } + + pBin[BinIndex] += 1.0; + + // update Total count + Total += 1.0; +} + +// Calculate cumulative chi^2 +double CACBinomialChi2::Calc() +{ + double Sum = 0; + for (int i=0; i<51; i++) + Sum += (pBin[i] * pBin[i])/(Total * pTable[i]); + + return Sum-Total; +} + +// Resets the cumulative test completely +void CACBinomialChi2::ResetAll(void) +{ + ZeroMemory(pBin, 51*sizeof(double)); + + Total = 0; +} + +// Initializes pTable +void CACBinomialChi2::InitPTable(void) +{ + pTable[0] = 0.0199136; + pTable[1] = 0.00443703; + pTable[2] = 0.00523732; + pTable[3] = 0.00613872; + pTable[4] = 0.00714366; + pTable[5] = 0.00825334; + pTable[6] = 0.00946896; + pTable[7] = 0.0107857; + pTable[8] = 0.0121988; + pTable[9] = 0.0136985; + pTable[10] = 0.0152738; + pTable[11] = 0.0169106; + pTable[12] = 0.0185897; + pTable[13] = 0.0202911; + pTable[14] = 0.021992; + pTable[15] = 0.0236668; + pTable[16] = 0.0252903; + pTable[17] = 0.0268349; + pTable[18] = 0.0282722; + pTable[19] = 0.0295779; + pTable[20] = 0.0307263; + pTable[21] = 0.0316944; + pTable[22] = 0.0324642; + pTable[23] = 0.0330194; + pTable[24] = 0.0333472; + pTable[25] = 0.0334425; + pTable[26] = 0.0333044; + pTable[27] = 0.0329349; + pTable[28] = 0.0323407; + pTable[29] = 0.031536; + pTable[30] = 0.0305358; + pTable[31] = 0.0293616; + pTable[32] = 0.0280353; + pTable[33] = 0.026582; + pTable[34] = 0.0250301; + pTable[35] = 0.0234031; + pTable[36] = 0.0217306; + pTable[37] = 0.020037; + pTable[38] = 0.018347; + pTable[39] = 0.0166828; + pTable[40] = 0.0150652; + pTable[41] = 0.0135096; + pTable[42] = 0.0120311; + pTable[43] = 0.0106397; + pTable[44] = 0.0093446; + pTable[45] = 0.00814988; + pTable[46] = 0.00705929; + pTable[47] = 0.00607215; + pTable[48] = 0.00518732; + pTable[49] = 0.00440057; + pTable[50] = 0.0200105; +} diff --git a/RandTest/ACBinomialChi2.h b/RandTest/ACBinomialChi2.h new file mode 100755 index 0000000..eca2ce9 --- /dev/null +++ b/RandTest/ACBinomialChi2.h @@ -0,0 +1,21 @@ +#pragma once +#include "BinomialChi2.h" + +class CACBinomialChi2 : + public CBinomialChi2 +{ +public: + CACBinomialChi2(void); + ~CACBinomialChi2(void); + // Overridden to initialze pTable for auto-correlation + virtual void InitPTable(void); + // Oberridden to calc index for ac binning + virtual void BinIt(double BinVal); + + virtual double Calc(); + virtual void ResetAll(); + +private: + double pBin[54]; + double pTable[54]; +}; diff --git a/RandTest/ACBinomialChi2.old.cpp b/RandTest/ACBinomialChi2.old.cpp new file mode 100755 index 0000000..6ad1ea2 --- /dev/null +++ b/RandTest/ACBinomialChi2.old.cpp @@ -0,0 +1,74 @@ +#include "StdAfx.h" +#include "acbinomialchi2.h" +#include "math.h" + +CACBinomialChi2::CACBinomialChi2(void) +{ + InitPTable(); +} + +CACBinomialChi2::~CACBinomialChi2(void) +{ +} + +// Overridden to initialze pTable for auto-correlation +VOID CACBinomialChi2::InitPTable(void) +{ + // 2^16 auto-correlation + pTable[0] = 0.0436123510244955321; + pTable[1] = 0.0494945236349887629; + pTable[2] = 0.0483966799920995590; + pTable[3] = 0.0465904838376547746; + pTable[4] = 0.0441572550221125783; + pTable[5] = 0.0412031166704026490; + pTable[6] = 0.0378513315493016954; + pTable[7] = 0.0342338133652554904; + pTable[8] = 0.0304826229422669920; + pTable[9] = 0.0267222012692252800; + pTable[10] = 0.0230629500401081457; + pTable[11] = 0.0195965713435074457; + pTable[12] = 0.0163933532864769463; + pTable[13] = 0.0135013695454263963; + pTable[14] = 0.0109473753452615719; + pTable[15] = 0.0087390490108419405; + pTable[16] = 0.0068681557366127721; + pTable[17] = 0.0053141974233424296; + pTable[18] = 0.0040481501288843967; + pTable[19] = 0.0030359644605488801; + pTable[20] = 0.0022415976474860848; + pTable[21] = 0.0016294434474942837; + pTable[22] = 0.0011661147954655219; + pTable[23] = 0.0008216056587517001; + pTable[24] = 0.0016958983342359457; +} + +// Oberridden to calc index for ac binning +VOID CACBinomialChi2::BinIt(double BinVal) +{ + // center bin (7 slots for ac) + int BinIndex = 0; + int Side = 1; // middle bin and right side + + // calc for autocorrelation + if ((BinVal<-3.) || (BinVal>3.)) + BinIndex = 1 + (abs((int)BinVal)-4)/8; // 8 slots per bin + + if (BinVal<-3.) + Side = -1; // left side + + // tails of distribution fall in one bin + if (BinIndex>24) + BinIndex = 24; + + if (Side>0) + { // fill in positive side + pBinPositive[BinIndex]+=1.; + } + else + { // fill in negative side + pBinNegative[BinIndex]+=1.; + } + + // update Total count + Total += 1.; +} diff --git a/RandTest/AutoCorrelation.cpp b/RandTest/AutoCorrelation.cpp new file mode 100755 index 0000000..c2fa982 --- /dev/null +++ b/RandTest/AutoCorrelation.cpp @@ -0,0 +1,125 @@ +//#include "StdAfx.h" +#include "AutoCorrelation.h" +#include "math.h" + +CAutoCorrelation::CAutoCorrelation(void) +: maxOrder(32) +, bitStream(0) +{ + CreateChi2Tests(); +// ZeroMemory(cumulativeACZScore, 32*sizeof(double)); + memset(cumulativeACZScore, 0, 32*sizeof(double)); + for (int i=0; i<32; i++) + P_Chi2[i] = .5; + ResetTest(); + // Initializing app must call ResetAll()!!! +} + +CAutoCorrelation::~CAutoCorrelation(void) +{ + for (int i=0; i<32; i++) + delete MetaChi2[i]; +} + +// Inserts a 32 bit word into the unit test +// Must insert word into Bias test first, AC is completely dependent upon Bias +void CAutoCorrelation::InsertWord32(uint32_t InWord32) +{ + bitStream >>= 32; + bitStream |= (uint64_t)InWord32 << 32; + + int Order; + for (Order=1; Order<=maxOrder; Order++) + { + totalAndCount[Order-1] += BitCount.GetCount32((uint32_t)bitStream & (uint32_t)(bitStream>>Order)); + blockXorCount[Order-1] += BitCount.GetCount32((uint32_t)bitStream ^ (uint32_t)(bitStream>>Order)); + } + + // Every 2048 32 bit words do a unit and cumulative calculation + if ((++blockWordCount)>=2048) + { + for (Order=1; Order<=maxOrder; Order++) + { + totalXorCount[Order-1] += blockXorCount[Order-1]; + double NN = (65536.*(totalBlockCount+1)); +// double XX = totalXorCount[Order-1]; + +// cumulativeACZScore[Order-1] = -(2*XX-NN) / sqrt(NN); + + // fractional ANDs + double fAnd = totalAndCount[Order-1] / NN; + // fractional bias + double fBias = *pTotalBiasCount / NN; + // full AC calculation (a-b^2) / ((1-b) * b) + if ((fBias*fBias)<=0 || fBias>=1) + cumulativeACZScore[Order-1] = 0; + else + cumulativeACZScore[Order-1] = sqrt(NN) * (fAnd - fBias*fBias)/(fBias - fBias*fBias); + + MetaChi2[Order-1]->Insert(blockXorCount[Order-1]); + P_Chi2[Order-1] = MetaChi2[Order-1]->GetPvalue();//Gamma.Gamma(50, ACChi2[Order-1].Calc()); + } + ResetTest(); + totalBlockCount++; + } +} + +// Resets current unit testing +void CAutoCorrelation::ResetTest(void) +{ + blockWordCount = 0; + for (int i=0; i<32; i++) + blockXorCount[i] = 0; +} + +// Resets cumulative scores +void CAutoCorrelation::ResetAll(double* pTotalBiasCount, double* pBlockBiasCount) +{ + this->pBlockBiasCount = pBlockBiasCount; + this->pTotalBiasCount = pTotalBiasCount; + + totalBlockCount = 0; +// ZeroMemory(totalAndCount, 32*sizeof(double)); + memset(totalAndCount, 0, 32*sizeof(double)); +// ZeroMemory(totalXorCount, 32*sizeof(double)); + memset(totalXorCount, 0, 32*sizeof(double)); + for (int i=0; i<32; i++) + P_Chi2[i] = .5; + + for (int i=0; i<32; i++) + { + MetaChi2[i]->Reset(); +// ACChi2[i].ResetAll(); + cumulativeACZScore[i] = 0.; + prevCumulativeACZScore[i] = 0.; + P_Chi2[i] = .5; + } + + bitStream = 0; + ResetTest(); +} + +void CAutoCorrelation::CreateChi2Tests(void) { + double pTable[] = { + 0.0016751591157826984, 0.0008125704223141780, 0.0011538175281394167, 0.0016130041110183372, 0.0022200160266961827, + 0.0030081473713346360, 0.0040129572987985991, 0.0052705078207208266, 0.0068149541705240074, 0.0086755302334023398, + 0.0108730598645528967, 0.0134162222471676280, 0.0162978932289932796, 0.0194919592594017396, 0.0229510396552503468, + 0.0266055419211793493, 0.0303644043131267252, 0.0341177483099464753, 0.0377414795440626524, 0.0411036575708872835, + 0.0440722296391575050, 0.0465235232774633227, 0.0483507488411394182, 0.0494717022548410683, 0.0467242519481981447, + 0.0494717022548410683, 0.0483507488411394182, 0.0465235232774633227, 0.0440722296391575050, 0.0411036575708872835, + 0.0377414795440626524, 0.0341177483099464753, 0.0303644043131267252, 0.0266055419211793493, 0.0229510396552503468, + 0.0194919592594017396, 0.0162978932289932796, 0.0134162222471676280, 0.0108730598645528967, 0.0086755302334023398, + 0.0068149541705240074, 0.0052705078207208266, 0.0040129572987985991, 0.0030081473713346360, 0.0022200160266961827, + 0.0016130041110183372, 0.0011538175281394167, 0.0008125704223141780, 0.0016751591157826984 + }; + double boundryTable[] = { + 32392, 32408, 32424, 32440, 32456, 32472, 32488, 32504, 32520, + 32536, 32552, 32568, 32584, 32600, 32616, 32632, 32648, 32664, 32680, + 32696, 32712, 32728, 32744, 32760, 32775, 32791, 32807, 32823, 32839, + 32855, 32871, 32887, 32903, 32919, 32935, 32951, 32967, 32983, 32999, + 33015, 33031, 33047, 33063, 33079, 33095, 33111, 33127, 33143, 1e100 + }; + + for (int i=0; i<32; i++) + MetaChi2[i] = new Chi2(49, pTable, false, boundryTable, false); +} diff --git a/RandTest/AutoCorrelation.h b/RandTest/AutoCorrelation.h new file mode 100755 index 0000000..cf74fc3 --- /dev/null +++ b/RandTest/AutoCorrelation.h @@ -0,0 +1,69 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "BitCount.h" +#include "ACBinomialChi2.h" +#include "Gamma.h" +#include "Chi2.hpp" + +class CAutoCorrelation +{ +public: + CAutoCorrelation(void); + ~CAutoCorrelation(void); + // Inserts a 32 bit word into the unit test + void InsertWord32(uint32_t InWord32); + + // Resets current unit testing + void ResetTest(void); + // Resets cumulative scores + void ResetAll(double* pTotalBiasCount, double* pBlockBiasCount); + + // Cumulative z-score + double cumulativeACZScore[32]; + double prevCumulativeACZScore[32]; + + // Keeps track of all blocks in cumulative testing + double totalBlockCount; + + // P-values of cumulative chi^2 tests + double P_Chi2[32]; + + // Highest order to be tested + int maxOrder; + +protected: + // History bit stream of last 32 plus current 32 bits (=64bits) + uint64_t bitStream; + + // First word in block to calc bias difference +// uint32_t firstBlockWord32; + + // Counts up "ands" for a block = multiply in cross correlation +// int blockAndCount[32]; + double blockXorCount[32]; + + // Counts up cumulative "ands" for each AC order + double totalXorCount[32]; + double totalAndCount[32]; + + // Chi^2 test cumulative for each order + void CreateChi2Tests(); + Chi2* MetaChi2[32]; + + // Incomplete Gamma function + CGamma Gamma; + + // Counts 32 bits words in current block test + int blockWordCount; + + // To quickly count bits within 32 bit word + CBitCount BitCount; + + // Pointer to the just calculated block bias one count + double* pBlockBiasCount; + // Pointer to the cumulative (including current block) one count + double* pTotalBiasCount; +}; diff --git a/RandTest/Bias.cpp b/RandTest/Bias.cpp new file mode 100755 index 0000000..72cec85 --- /dev/null +++ b/RandTest/Bias.cpp @@ -0,0 +1,85 @@ +//#include "StdAfx.h" +#include "Bias.h" +#include "math.h" + +CBias::CBias(void) +: cumulativeBiasZScore(0) +, P_Chi2(.5) +{ + CreateChi2Test(); + + ResetAll(); +} + +CBias::~CBias(void) +{ + delete MetaChi2; +} + +// Tests a 32 bit word +void CBias::InsertWord32(uint32_t InWord32) +{ + blockBiasCount += BitCount.GetCount32(InWord32); + + // Every 2^16 bits calc cumulative bias and chi^2 + if ( (++blockWordCount)>=2048 ) + { + // calculate unit bias + double BlockZScore = (2.*blockBiasCount-(32*blockWordCount))/(sqrt(32.0*blockWordCount)); + + // Chi^2 for bias using binomial chi^2 test Chi2Binomial + MetaChi2->Insert(blockBiasCount); + P_Chi2 = MetaChi2->GetPvalue(); + + totalBiasCount += blockBiasCount; + formerBlockBiasCount = blockBiasCount; + totalBlockCount++; + + cumulativeBiasZScore = ((2.*totalBiasCount) - (65536.*totalBlockCount)) / sqrt(65536.*totalBlockCount); + + ResetTest(); + } +} + +// Resets current test +void CBias::ResetTest() +{ + blockWordCount = 0; + blockBiasCount = 0; +} + +// Resets all +void CBias::ResetAll() +{ + P_Chi2 = .5; + cumulativeBiasZScore = 0.; + totalBiasCount = 0; + totalBlockCount = 0; +// BiasChi2.ResetAll(); + MetaChi2->Reset(); + ResetTest(); +} + +void CBias::CreateChi2Test(void) { + double pTable[] = { + 0.0016751591157826984, 0.0008125704223141780, 0.0011538175281394167, 0.0016130041110183372, 0.0022200160266961827, + 0.0030081473713346360, 0.0040129572987985991, 0.0052705078207208266, 0.0068149541705240074, 0.0086755302334023398, + 0.0108730598645528967, 0.0134162222471676280, 0.0162978932289932796, 0.0194919592594017396, 0.0229510396552503468, + 0.0266055419211793493, 0.0303644043131267252, 0.0341177483099464753, 0.0377414795440626524, 0.0411036575708872835, + 0.0440722296391575050, 0.0465235232774633227, 0.0483507488411394182, 0.0494717022548410683, 0.0467242519481981447, + 0.0494717022548410683, 0.0483507488411394182, 0.0465235232774633227, 0.0440722296391575050, 0.0411036575708872835, + 0.0377414795440626524, 0.0341177483099464753, 0.0303644043131267252, 0.0266055419211793493, 0.0229510396552503468, + 0.0194919592594017396, 0.0162978932289932796, 0.0134162222471676280, 0.0108730598645528967, 0.0086755302334023398, + 0.0068149541705240074, 0.0052705078207208266, 0.0040129572987985991, 0.0030081473713346360, 0.0022200160266961827, + 0.0016130041110183372, 0.0011538175281394167, 0.0008125704223141780, 0.0016751591157826984 + }; + double boundryTable[] = { + 32392, 32408, 32424, 32440, 32456, 32472, 32488, 32504, 32520, + 32536, 32552, 32568, 32584, 32600, 32616, 32632, 32648, 32664, 32680, + 32696, 32712, 32728, 32744, 32760, 32775, 32791, 32807, 32823, 32839, + 32855, 32871, 32887, 32903, 32919, 32935, 32951, 32967, 32983, 32999, + 33015, 33031, 33047, 33063, 33079, 33095, 33111, 33127, 33143, 1e100 + }; + + MetaChi2 = new Chi2(49, pTable, false, boundryTable, false); +} diff --git a/RandTest/Bias.h b/RandTest/Bias.h new file mode 100755 index 0000000..8cbd2ea --- /dev/null +++ b/RandTest/Bias.h @@ -0,0 +1,57 @@ +// 1/0 Balance - expected value is p(0) = p(1) = 0.5 + +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "BitCount.h" +#include "Chi2.hpp" +#include "Gamma.h" + +class CBias +{ + friend class CBiasAndAC; + +public: + CBias(void); + ~CBias(void); + + // Tests a 32 bit word + void InsertWord32(uint32_t InWord32); + + // Resets all + void ResetAll(); + + // Cumulative z score + double cumulativeBiasZScore; + + // Cumulative chi^2 results + double P_Chi2; + + // Keeps track of all blocks in cumulative testing + double totalBlockCount; + +protected: + // Counts 1-bits per block + double blockBiasCount; + // Cumulative bias 1-bit count + double totalBiasCount; + double formerBlockBiasCount; + + // Resets current test + void ResetTest(); + + // Quick 32 bit word bit count table lookup + CBitCount BitCount; + + // Keeps track of words tested in block + int blockWordCount; + + // Calculates the chi^2 + void CreateChi2Test(); + Chi2* MetaChi2; + + // Incomplete Gamma function + CGamma Gamma; +}; diff --git a/RandTest/BiasAndAC.cpp b/RandTest/BiasAndAC.cpp new file mode 100755 index 0000000..4dc51dc --- /dev/null +++ b/RandTest/BiasAndAC.cpp @@ -0,0 +1,25 @@ +//#include "StdAfx.h" +#include "BiasAndAC.h" + +CBiasAndAC::CBiasAndAC(void) +{ + ResetAll(); +} + +CBiasAndAC::~CBiasAndAC(void) +{ +} + +// Inserts a 32 bit word into both Bias an AC tests +void CBiasAndAC::InsertWord32(uint32_t InWord32) +{ + Bias.InsertWord32(InWord32); + AC.InsertWord32(InWord32); +} + +// Resets all in Bias and AutoCorrelation +void CBiasAndAC::ResetAll(void) +{ + Bias.ResetAll(); + AC.ResetAll(&Bias.totalBiasCount, &Bias.formerBlockBiasCount); +} diff --git a/RandTest/BiasAndAC.h b/RandTest/BiasAndAC.h new file mode 100755 index 0000000..8dbf743 --- /dev/null +++ b/RandTest/BiasAndAC.h @@ -0,0 +1,24 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "Bias.h" +#include "AutoCorrelation.h" + +class CBiasAndAC +{ +public: + CBiasAndAC(void); + ~CBiasAndAC(void); + // Inserts a 32 bit word into both Bias an AC tests + void InsertWord32(uint32_t InWord32); + + // Bias test + CBias Bias; + // Auto-correlation test + CAutoCorrelation AC; + + // Resets all in Bias and AutoCorrelation + void ResetAll(void); +}; diff --git a/RandTest/BinomialChi2.cpp b/RandTest/BinomialChi2.cpp new file mode 100755 index 0000000..3fccebb --- /dev/null +++ b/RandTest/BinomialChi2.cpp @@ -0,0 +1,97 @@ +#include "StdAfx.h" +#include "binomialchi2.h" +#include "math.h" + +CBinomialChi2::CBinomialChi2(void) +{ + InitPTable(); +} + +CBinomialChi2::~CBinomialChi2(void) +{ +} + +// Bins incomming values for future Chi2 test +void CBinomialChi2::BinIt(double BinVal) +{ + // center bin (15 slots for 1/0) + int BinIndex = 0; + int Side = 1; // middle bin and right side + + // calc Index for 1/0 balance + if ((BinVal>32775.) || (BinVal<32761.)) + BinIndex = 1 + (abs((int)BinVal-32768)-8)/16; // 16 slots per bin + + if (BinVal<32761) + Side = -1; // left side + + // tails of distribution fall in one bin + if (BinIndex>24) + BinIndex = 24; + + if (Side>0) + { // fill in positive side + pBinPositive[BinIndex]+=1.; + } + else + { // fill in negative side + pBinNegative[BinIndex]+=1.; + } + + // update Total count + Total += 1.; +} + +// Calculate cumulative chi^2 +double CBinomialChi2::Calc() +{ + double Sum = (pBinPositive[0] * pBinPositive[0])/(Total * pTable[0]); + + for (int i=1; i<25; i++) + { + Sum += (pBinPositive[i] * pBinPositive[i])/(Total * pTable[i]); + Sum += (pBinNegative[i] * pBinNegative[i])/(Total * pTable[i]); + } + + return Sum-Total; +} + +// Resets the cumulative test completely +void CBinomialChi2::ResetAll(void) +{ + ZeroMemory(pBinPositive, 30*sizeof(double)); + ZeroMemory(pBinNegative, 30*sizeof(double)); + + Total = 0; +} + +// Initializes pTable +void CBinomialChi2::InitPTable(void) +{ + // 2^16 1/0 bias distribution table + pTable[0] = 0.0467242519481981447; + pTable[1] = 0.0494717022548410683; + pTable[2] = 0.0483507488411394182; + pTable[3] = 0.0465235232774633227; + pTable[4] = 0.0440722296391575050; + pTable[5] = 0.0411036575708872835; + pTable[6] = 0.0377414795440626524; + pTable[7] = 0.0341177483099464753; + pTable[8] = 0.0303644043131267252; + pTable[9] = 0.0266055419211793493; + pTable[10] = 0.0229510396552503468; + pTable[11] = 0.0194919592594017396; + pTable[12] = 0.0162978932289932796; + pTable[13] = 0.0134162222471676280; + pTable[14] = 0.0108730598645528967; + pTable[15] = 0.0086755302334023398; + pTable[16] = 0.0068149541705240074; + pTable[17] = 0.0052705078207208266; + pTable[18] = 0.0040129572987985991; + pTable[19] = 0.0030081473713346360; + pTable[20] = 0.0022200160266961827; + pTable[21] = 0.0016130041110183372; + pTable[22] = 0.0011538175281394167; + pTable[23] = 0.0008125704223141780; + pTable[24] = 0.0016751591157826984; +} diff --git a/RandTest/BinomialChi2.h b/RandTest/BinomialChi2.h new file mode 100755 index 0000000..c7302c8 --- /dev/null +++ b/RandTest/BinomialChi2.h @@ -0,0 +1,29 @@ +#pragma once + +class CBinomialChi2 +{ +public: + CBinomialChi2(void); + ~CBinomialChi2(void); + + // Bins incomming values for future Chi2 test + virtual void BinIt(double BinVal); + // Calculate cumulative chi^2 + double Calc(void); + + // Resets the cumulative test completely + void ResetAll(void); + +protected: + // Positive side bins + double pBinPositive[30]; + // Negative side bins + double pBinNegative[30]; + // Total values binned + double Total; + + // Table with bin probabilities + double pTable[25]; + // Initializes pTable + virtual void InitPTable(void); +}; diff --git a/RandTest/BitCount.cpp b/RandTest/BitCount.cpp new file mode 100755 index 0000000..eca4e87 --- /dev/null +++ b/RandTest/BitCount.cpp @@ -0,0 +1,29 @@ +//#include "StdAfx.h" +#include "BitCount.h" + +uint8_t CBitCount::Table16[65536]; +bool CBitCount::IsInitialized = false; + +CBitCount::CBitCount(void) +{ + // Initialize, if not already done so + if (!IsInitialized) + { + for (int i=0; i<=65535; i++) + { + // Fill up table + uint8_t BitCount = 0; + for (int b=0; b<16; b++) + BitCount += ((i>>b)&0x1); + + Table16[i] = BitCount; + } + + IsInitialized = true; + } +} + +CBitCount::~CBitCount(void) +{ +} + diff --git a/RandTest/BitCount.h b/RandTest/BitCount.h new file mode 100755 index 0000000..b20ae5a --- /dev/null +++ b/RandTest/BitCount.h @@ -0,0 +1,24 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) + +class CBitCount +{ +public: + CBitCount(void); + ~CBitCount(void); + + // Get bitcount for this 32 bit word +// __forceinline uint8_t GetCount32(uint32_t InWord32) { + uint8_t GetCount32(uint32_t InWord32) { +// PUSHORT pWordDiv16 = (PUSHORT)&InWord32; + return (Table16[(uint16_t)InWord32] + Table16[(uint16_t)(InWord32>>16)]); + } + +private: + // Table to keep bit counts of 16 bit values + static uint8_t Table16[65536]; + static bool IsInitialized; +}; diff --git a/RandTest/Chi2.cpp b/RandTest/Chi2.cpp new file mode 100755 index 0000000..d70e265 --- /dev/null +++ b/RandTest/Chi2.cpp @@ -0,0 +1,39 @@ +#include "StdAfx.h" +#include "chi2.h" + +CChi2::CChi2(void) +{ + Init(); +} + +CChi2::~CChi2(void) +{ +} + +// Clear Chi2 test and set number of bins to use +void CChi2::Init(double nBins) +{ + nCount = 0; // nCount + this->nBins = nBins; // nBins + pBin = 1./nBins; // pBin = 1/nBins + SumBinsSquared = 0; // SumBinsSquared + memset( Bin, 0, 65536*4 ); // Zero all Bins +} + +// Inserts a probability and calculates new chi^2 +double CChi2::CalcChi2(double pValue) +{ + double np; + double BinAdd; + UINT BinIndex; + + nCount++; + np = nCount*pBin; + BinIndex = (UINT)(nBins * pValue); + + BinAdd = 2 * Bin[BinIndex] + 1; + Bin[BinIndex]++; + SumBinsSquared += BinAdd; + + return (SumBinsSquared/np - nCount); +} diff --git a/RandTest/Chi2.h b/RandTest/Chi2.h new file mode 100755 index 0000000..75496b0 --- /dev/null +++ b/RandTest/Chi2.h @@ -0,0 +1,20 @@ +#pragma once + +class CChi2 +{ +public: + CChi2(void); + ~CChi2(void); + + // Clear Chi2 test and set number of bins to use + void Init(double nBins=32.); + // Inserts a probability and calculates new chi^2 + double CalcChi2(double pValue); + +protected: + double Bin[65536]; + double SumBinsSquared; + double pBin; // Probability of falling into any given bin + double nBins; // Total number of bins + double nCount; // Number of pValues binned +}; diff --git a/RandTest/Chi2.hpp b/RandTest/Chi2.hpp new file mode 100755 index 0000000..06986ac --- /dev/null +++ b/RandTest/Chi2.hpp @@ -0,0 +1,104 @@ +#pragma once + +#include +#include "Gamma.h" + + +class Chi2 { +public: + Chi2(int binCount=10, double* pTable=0, bool pTableCumulative=false, double* boundryTable=0, bool reverseBoundrySense=false) { + this->binCount = binCount; + this->binned = new double[binCount]; + this->boundryTable = new double[binCount]; + this->reverseBoundrySense = reverseBoundrySense; + this->pTable = new double[binCount]; + + if (pTable == 0) { + for (int i=0; ipTable[i] = 1.0 / binCount; + this->boundryTable[i] = (1.0+i) / binCount; + } + } + else { + memcpy(this->pTable, pTable, sizeof(double)*binCount); + memcpy(this->boundryTable, boundryTable, sizeof(double)*binCount); + } + + if (pTableCumulative == true) { + for (int i=(binCount-1); i>=1; i--) { + this->pTable[i] = pTable[i] - pTable[i-1]; + } + } + + Reset(); + } + + ~Chi2() { + delete[] pTable; + delete[] boundryTable; + delete[] binned; + } + + void Reset() { + doRecalc = true; + for (int i=0; i boundryTable[i]) + i++; + } + else { + while (inValue < boundryTable[i]) + i++; + } + binned[i]++; + doRecalc = true; + } + + double GetChi2() { + if (doRecalc == true) + Recalc(); + + return chi2; + } + + double GetPvalue() { + if (doRecalc == true) + Recalc(); + + return (Gamma.Gamma(binCount-1, chi2)); + } + +private: + void Recalc() { + double sum = 0; + double total = 0; + + for (int i=0; iBitCount = BitCount; + + // get initialization blocks used + Q = qFactor * 256; // 20*256 = 5120 + + // get total number of m-bit blocks + nb = (unsigned long)floor( BitCount / 8.0 ); // 524288 + + // get number of calculation blocks used + K = nb - Q; // 519168 + + // calc sd + sd = CalcStandardDeviation( K, 8 ); + + // zero out tab + memset( tab, 0, 256*4 ); + + blockWordCount = 0; + Sum = 0; +} + +void CEntropy::InsertWord32(uint32_t InWord32) +{ + int i; + uint8_t RndByte; + + for ( i=0; i<=3; i++ ) + { + RndByte = (uint8_t)(InWord32>>(8*i)); + + blockWordCount++; + + if ( blockWordCount<=Q ) + { + tab[RndByte] = blockWordCount; + } + else + { + Sum += fcoef( blockWordCount - tab[RndByte] ); + tab[RndByte] = blockWordCount; + } + } + + if ( (blockWordCount*8)>=BitCount ) + { + H = Sum/((double)K); + + ETotal += H/8.; + totalBlockCount++; + E = ETotal / (totalBlockCount); + double BlockZScore = ( (H-8.)/sd ); + ZScoreTotal += BlockZScore; + MetaChi2->Insert(Sum); + P_Chi2 = MetaChi2->GetPvalue(); + if (totalBlockCount!=0) + cumulativeZScore = ZScoreTotal / sqrt(totalBlockCount); + Initialize( 4194304, 20 ); + } +} + +double CEntropy::fcoef(uint32_t i) +{ + // set constants + const double l2 = log(2.0); + const double c = -0.8327462; + const int limit = 23; + + double retval; // return value + unsigned long kk; // universal index + int j; // intermediate calc value + + retval = 0; + if ( iReset(); +// Chi2.Init(); + Initialize( 4194304, 20 ); +} + +void CEntropy::CreateChi2Test(void) { + +double pTable[] = {0.02012108166,0.04008956584,0.06018692852,0.08018501372,0.1001359343,0.120192745,0.1401963459,0.16005801,0.1802017087,0.200031418,0.2200301029,0.2399762859,0.2599755987,0.2800970982,0.3000327584,0.3202120847,0.340234622,0.360070732,0.3800778627,0.3998619416,0.4199041271,0.4398490175,0.4601300439,0.4798811487,0.5000167861,0.5201045044,0.5399738275,0.5600664317,0.5798227432,0.5998761946,0.619978835,0.6402984742,0.6599255874,0.6799283734,0.7001158234,0.7199709181,0.7401361424,0.7599607824,0.7802489806,0.8001210836,0.8200605293,0.8400672073,0.8601944042,0.880037642,0.900350094,0.9203249058,0.9401774272,0.9600776027,0.9807874005,1.0}; + + double boundryTable[] = { + 4151655.83724745500, 4151904.48938296180, 4152066.26934954900, 4152189.24478157190, 4152290.52925698830, + 4152378.39198942950, 4152456.23714622110, 4152526.37086037970, 4152591.89903046660, 4152652.04845737200, + 4152709.17705694870, 4152763.25623371170, 4152815.06782631900, 4152865.15788602550, 4152913.08992392620, + 4152960.15932061100, 4153005.64889538610, 4153049.71294774950, 4153093.31400005010, 4153135.74272596740, + 4153178.16347156510, 4153219.94024266860, 4153262.09135455310, 4153302.92978435200, 4153344.45511117900, + 4153385.87801960110, 4153426.95058751110, 4153468.69217352150, 4153510.04595550760, 4153552.45207967800, + 4153595.52232766520, 4153639.76729502670, 4153683.32874056320, 4153728.72452108000, 4153775.75883651480, + 4153823.43728935770, 4153873.56984770070, 4153924.85234509460, 4153979.78749536580, 4154036.49777039420, + 4154096.92328405190, 4154161.94629989700, 4154232.98319796660, 4154310.19911737930, 4154399.23572011480, + 4154500.86437934400, 4154623.25103316270, 4154783.62485346620, 4155045.02386120380, 1e100 + }; + + MetaChi2 = new Chi2(50, pTable, true, boundryTable, false); +} diff --git a/RandTest/Entropy.h b/RandTest/Entropy.h new file mode 100755 index 0000000..2f390e4 --- /dev/null +++ b/RandTest/Entropy.h @@ -0,0 +1,45 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "Gamma.h" +#include "Chi2.hpp" +#include "Stat.h" + +class CEntropy +{ +public: + double P_Chi2; + inline double fcoef(uint32_t i); + double sd; + double E; + double H; + double cumulativeZScore; + double totalBlockCount; + uint32_t blockWordCount; + void InsertWord32(uint32_t InWord32); + void Initialize(uint32_t BitCount, uint32_t qFactor); + CEntropy(); + virtual ~CEntropy(); + +protected: + CStat Stat; + double pow2_32; + CGamma Gamma; + double ETotal; + double ZScoreTotal; + double Sum; + uint32_t K; + uint32_t nb; + uint32_t Q; + uint32_t BitCount; + uint32_t tab[256]; + double CalcStandardDeviation(uint32_t K, uint32_t m); + Chi2* MetaChi2; + void CreateChi2Test(void); + +public: + // Reset cumulative test + void ResetAll(void); +}; diff --git a/RandTest/Gamma.cpp b/RandTest/Gamma.cpp new file mode 100755 index 0000000..ff0101a --- /dev/null +++ b/RandTest/Gamma.cpp @@ -0,0 +1,218 @@ +//#include "StdAfx.h" +#include "Gamma.h" +#include "math.h" +#include + +CGamma::CGamma() +{ + +} + +CGamma::~CGamma() +{ + +} + +double CGamma::gammln(double xx) +{ + double RetVal; + + double y; + double x; + double tmp; + double ser; + int j; + + double cof[6]; + + cof[0] = 76.18009172947146; cof[1] = -86.50532032941677; cof[2] = 24.01409824083091; + cof[3] = -1.231739572450155; cof[4] = 1.208650973866179e-3; cof[5] = -5.395239384953e-6; + + y = x = xx; + + tmp = (x+5.5) - (x+0.5)*(log(x+5.5)); + ser = 1.000000000190015; + + for ( j=0; j<=5; j++ ) + { + y += 1.0; + ser += cof[j]/y; + } + + RetVal = log( 2.506628274631*ser/x ) - tmp; + + return RetVal; +} + +double CGamma::gamser(double a, double x) +{ + double RetVal; + + double gln; + double ap; + double del; + double sm; + int ITMAX; + int m; + + if ( x!=0 ) + { + ITMAX = 31; + gln = gammln( a ); + ap = a; + del = sm = 1.0/a; + + for( m=1; m<=ITMAX; m++ ) + { + ap += 1; + del *= x/ap; + sm += del; + } + + RetVal = sm * exp( a*log(x) - x - gln ); + } + else + RetVal = 0; + + return RetVal; +} + +double CGamma::gcf(double a, double x) +{ + double RetVal; + + double gln; + double b; + double c; + double d; + double h; + double an; + double del; + double FPMIN; + int ITMAX; + int i; + + gln = gammln( a ); + b = x + 1.0 - a; + FPMIN = pow( 10., -30 ); + c = 1.0 / FPMIN; + d = 1.0 / b; + h = d; + + ITMAX = 30; + for ( i=1; i<=ITMAX; i++ ) + { + an = (-i) * (i-a); + b += 2.0; + d = an*d + b; + if ( abs((int)d) < FPMIN ) d = FPMIN; + c = b + an/c; + if ( abs((int)c) < FPMIN ) c = FPMIN; + d = 1.0/d; + del = d * c; + h *= del; + } + + RetVal = 1.0 - (exp( a*log(x) - x - gln ) * h); + + return RetVal; +} + +double CGamma::Gamma(double a, double chi2) +{ + double RetVal; + + double x; + + a = a / 2; + x = chi2 / 2; + + if ( x<=a ) + RetVal = gamser( a, x ); + else + RetVal = gcf( a, x ); + + return 1-RetVal; +} + +double CGamma::Beta(double z, double a, double b) +{ + double bt; + + if (z < 0.0) + z = 0.0; + if (z > 1.0) + z = 1.0; + + if (z == 0.0 || z == 1.0) + bt=0.0; + else + // Factors in front of the continued fraction. + bt = exp( gammln(a+b) - gammln(a) - gammln(b) + a*log(z) + b*log(1.0-z) ); + + if (z < (a+1.0)/(a+b+2.0)) // Use continued fraction directly. + return bt*betacf(a, b, z)/a; + else // Use continued fraction after making the symmetry transformation. + return 1.0-bt*betacf(b, a, 1.0-z)/b; +} + +double CGamma::betacf(double a, double b, double x) +{ + int m,m2; + double aa,c,d,del,h,qab,qam,qap; + + double MAXIT = 100; + double EPS = 3.0e-7; + double FPMIN = 1.0e-30; + + qab = a+b; + qap = a+1.0; + qam = a-1.0; + + c = 1.0; // First step of Lentz's method. + d = 1.0 - qab*x/qap; + + if (fabs(d) < FPMIN) + d = FPMIN; + + d = 1.0/d; + h = d; + + for (m=1; m<=MAXIT; m++) +// while (true) + { + m2 = 2*m; + aa = m*(b-m)*x / ((qam+m2)*(a+m2)); + d = 1.0 + aa*d; // One step (the even one) of the recurrence. + + if (fabs(d) < FPMIN) + d=FPMIN; + + c = 1.0 + aa/c; + + if (fabs(c) < FPMIN) + c = FPMIN; + + d = 1.0/d; + h *= d*c; + aa = -(a+m)*(qab+m)*x / ((a+m2)*(qap+m2)); + d = 1.0 + aa*d; // Next step of the recurrence (the odd one). + + if (fabs(d) < FPMIN) + d = FPMIN; + + c = 1.0 + aa/c; + + if (fabs(c) < FPMIN) + c = FPMIN; + + d = 1.0/d; + del = d*c; + h *= del; + + if (fabs(del-1.0) < EPS) // Are we done? + break; + } + + return h; +} diff --git a/RandTest/Gamma.h b/RandTest/Gamma.h new file mode 100755 index 0000000..136239c --- /dev/null +++ b/RandTest/Gamma.h @@ -0,0 +1,15 @@ +#pragma once + +class CGamma +{ +public: + double betacf(double a, double b, double x); + double Beta(double z, double a, double b); + double Gamma( double a, double chi2 ); + CGamma(); + virtual ~CGamma(); +protected: + double gcf( double a, double x ); + double gamser( double a, double x ); + double gammln( double xx ); +}; diff --git a/RandTest/KolmogorovSmirnov.cpp b/RandTest/KolmogorovSmirnov.cpp new file mode 100755 index 0000000..e1629f7 --- /dev/null +++ b/RandTest/KolmogorovSmirnov.cpp @@ -0,0 +1,249 @@ +//#include "StdAfx.h" +#include "KolmogorovSmirnov.h" +#include "memory.h" +#include "stdlib.h" +#include "math.h" +#include + +CKolmogorovSmirnov::CKolmogorovSmirnov() +{ +} + + +/////////////////////////////////////////////////////////////////////////////// +// double CKolmogorovSmirnov::ZtoP +// +// Returns: cumulative normal distribution value based on input z-score +// Accuracy better than 1% for zscore<=±7.5; better than 0.05% for zscore=±4 +double CKolmogorovSmirnov::ZtoP( double zscore ) +{ + double retval; + + // calculation variables + double w; + double y; + double t; + double num; + double denom; + + // calculation constants + double c[8]; + c[1] = 2.506628275; c[2] = 0.31938153; c[3] = -0.356563782; c[4] = 1.781477937; + c[5] = -1.821255978; c[6] = 1.330274429; c[7] = 0.2316419; + + w = (zscore>=0)? 1 : -1; + y = 1./( 1. + (c[7]*w*zscore) ); + t = 1. + (c[7]*w*zscore); + + num = w * ( 0.5 - c[2] + ((c[6] + (c[5]*t) + (c[4]*t*t) + (c[3]*t*t*t))/(t*t*t*t)) ); + denom = c[1] * pow( 10, pow(0.5*zscore, 2) ); + + retval = 0.5 + (num/denom); + + return retval; +} + + +/////////////////////////////////////////////////////////////////////////////// +// double CKolmogorovSmirnov::PtoZ +// +// Returns: z-score based on input cumulative normal distribution value +// Accuracy better than 0.1% for zscore<=±7.5; 6 digits for zscore=±6 +double CKolmogorovSmirnov::PtoZ( double p ) +{ + double retval = -8.2; + if (p <= DBL_EPSILON) + return retval; + + // calculation variables + double pp; + double y; + double num; + double denom; + + // calculation constants + double P[5]; + P[0] = -0.322232431088; P[1] = -1.0; P[2] = -0.342242088547; + P[3] = -0.0204231210245; P[4] = -0.453642210148e-4; + double q[5]; + q[0] = 0.099348462606; q[1] = 0.588581570495; q[2] = 0.531103462366; + q[3] = 0.10353775285; q[4] = 0.38560700634e-2; + + pp = (p<0.5)? p : (1.-p); + y = sqrt( log(1./(pp*pp)) ); + num = y*(y*(y*(y*P[4]+P[3]) + P[2]) + P[1]) + P[0]; + denom = y*(y*(y*(y*q[4]+q[3]) + q[2]) + q[1]) + q[0]; + + retval = y + (num/denom); + retval = (p<0.5)? -retval : retval; + + return retval; +} + + +void CKolmogorovSmirnov::KSUP( double* pKn_pos, double* pKn_neg, double* U, unsigned long n ) +{ + unsigned long i; // universal index + double val; // calc value + double* Us = new double[n]; + double Kn_pos; + double Kn_neg; + double* cump1 = new double[n+2]; + double* cump2 = new double[n+2]; + double* kn = new double[n+2]; + double max; + double min; + + memcpy( (void*)Us, (void*)U, n * sizeof(double) ); + qsort( Us, n, sizeof(double), KSUPcompare ); + + cump1[1] = Us[0]; + cump2[1] = 0.; + for ( i=1; i<=n; i++ ) + { + cump1[i+1] = Us[i-1]; + val = (double)i/n; + cump2[i+1] = val; + } + + min = cump2[2] - cump1[2]; + max = min; + for ( i=2; i<=(n+1); i++ ) + { + kn[i-1] = cump2[i] - cump1[i]; + if ( kn[i-1] > max ) + { + max = kn[i-1]; + } + if ( kn[i-1] < min ) + { + min = kn[i-1]; + } + } + + Kn_pos = sqrt((double)n) * (double)max; + Kn_neg = (-sqrt((double)n)) * ((double)min - (1/(double)n)); + + *pKn_pos = KSProb( n, Kn_pos ); + *pKn_neg = KSProb( n, Kn_neg ); + + if (Us!=NULL) + { + delete Us; + Us = NULL; + } + + if (kn!=NULL) + { + delete kn; + kn = NULL; + } + + if (cump1!=NULL) + { + delete cump1; + cump1 = NULL; + } + + if (cump2!=NULL) + { + delete cump2; + cump2 = NULL; + } +} + + +double CKolmogorovSmirnov::KSProb( unsigned long n, double Kn ) +{ + double retval; + double e; + double las; + unsigned long j; + double cof[7]; + + cof[1] = 76.18009172947146; cof[2] = -86.50532032941677; cof[3] = 24.01409824083091; + cof[4] = -1.231739572450155; cof[5] = 1.208650973866179e-3; cof[6] = -5.395239384953e-6; + + if ( Kn<=0 ) + { + retval = 0; + } + else if ( Kn>=sqrt((double)n) ) + { + retval = 1; + } + else + { + e = Kn / sqrt((double)n); + las = floor( (double)n - ((double)n*e) ); + + retval = 0; + for ( j=0; j<=las; j++ ) + { + double ee = exp(lnbin(n, j)); + double p1 = pow((e+(double)j/(double)n),((double)j-1.)); + double p2 = pow((1.0-e-(double)j/(double)n),(double)(n-j)); + retval += ee * p1 * p2; + } + + retval = 1.-e*retval; + } + + return retval; +} + +double CKolmogorovSmirnov::lnbin( double n, double k ) +{ + double retval = 0; + if (k==0) + retval = 0; + else + { + retval = lnf(n) - lnf(k) - lnf(n-k); + } + + return retval; +} + + +double CKolmogorovSmirnov::lnf( double xx ) +{ + double retval = 0; + double x1 = xx + 1.; + + double cof[7]; + + cof[1] = 76.18009172947146; cof[2] = -86.50532032941677; cof[3] = 24.01409824083091; + cof[4] = -1.231739572450155; cof[5] = 1.208650973866179e-3; cof[6] = -5.395239384953e-6; + + if (x1<=1.) + retval = 0; + else + { + double x = x1; + double y = x1; + + double tmp = x + 5.5 - (x+.5)*log(x+5.5); + double ser = 1.000000000190015; + + for (int j=0; j<=5; j++) + { + y = y+1.; + ser += cof[j+1]/y; + } + + retval = log(2.506628274631 * ser/x) - tmp; + } + + return retval; +} + +int CKolmogorovSmirnov::KSUPcompare( const void* elem1, const void* elem2 ) +{ + if ((*(double*)elem1)==(*(double*)elem2)) + return 0; + if ((*(double*)elem1)>(*(double*)elem2)) + return 1; + else + return -1; +} diff --git a/RandTest/KolmogorovSmirnov.h b/RandTest/KolmogorovSmirnov.h new file mode 100755 index 0000000..dfd386b --- /dev/null +++ b/RandTest/KolmogorovSmirnov.h @@ -0,0 +1,20 @@ +#pragma once + +#include + +class CKolmogorovSmirnov +{ +public: + CKolmogorovSmirnov(); + +public: + double ZtoP( double zscore ); + double PtoZ( double p ); + void KSUP( double* pKn_pos, double* pKn_neg, double* U, unsigned long n ); + double lnf( double xx ); + double lnbin( double n, double k ); + double KSProb( unsigned long n, double Kn ); + +private: + static int KSUPcompare( const void* elem1, const void* elem2 ); +}; \ No newline at end of file diff --git a/RandTest/MT19937.cpp b/RandTest/MT19937.cpp new file mode 100755 index 0000000..f8a3f8e --- /dev/null +++ b/RandTest/MT19937.cpp @@ -0,0 +1,73 @@ +#include "StdAfx.h" +#include "mt19937.h" + +CMT19937::CMT19937() +{ + mti = N+1; + mag01[0] = 0; + + /* mag01[x] = x * MATRIX_A for x=0,1 */ + mag01[1] = MATRIX_A; +} + +CMT19937::~CMT19937() +{ + +} + +void CMT19937::SGenRand(uint32_t seed) +{ + /* setting initial seeds to mt[N] using */ + /* the generator Line 25 of Table 1 in */ + /* [KNUTH 1981, The Art of Computer Programming */ + /* Vol. 2 (2nd Ed.), pp102] */ +// mt[0]= seed & 0xffffffff; +// for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } + +} + +uint32_t CMT19937::GenRand() +{ + uint32_t y; + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if sgenrand() has not been called, */ + SGenRand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1]; + } + for (;kk> 1) ^ mag01[y & 0x1]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]; + + mti = 0; + } + + y = mt[mti++]; + y ^= TEMPERING_SHIFT_U(y); + y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B; + y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C; + y ^= TEMPERING_SHIFT_L(y); + + return y; +} diff --git a/RandTest/MT19937.h b/RandTest/MT19937.h new file mode 100755 index 0000000..18da894 --- /dev/null +++ b/RandTest/MT19937.h @@ -0,0 +1,64 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0df /* constant vector a */ +#define UPPER_MASK 0x80000000 /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffff /* least significant r bits */ + +/* Tempering parameters */ +#define TEMPERING_MASK_B 0x9d2c5680 +#define TEMPERING_MASK_C 0xefc60000 +#define TEMPERING_SHIFT_U(y) (y >> 11) +#define TEMPERING_SHIFT_S(y) (y << 7) +#define TEMPERING_SHIFT_T(y) (y << 15) +#define TEMPERING_SHIFT_L(y) (y >> 18) + +class CMT19937 +{ +public: + uint32_t GenRand(); + CMT19937(); + virtual ~CMT19937(); + void SGenRand(uint32_t seed); + +private: + uint32_t mag01[2]; + int mti; + uint32_t mt[N]; +}; + +// C++ encapsulated modified from: +// +/* A C-program for MT19937: Integer version */ +/* genrand() generates one pseudorandom unsigned integer (32bit) */ +/* which is uniformly distributed among 0 to 2^32-1 for each */ +/* call. sgenrand(seed) set initial values to the working area */ +/* of 624 words. Before genrand(), sgenrand(seed) must be */ +/* called once. (seed is any 32-bit integer except for 0). */ +/* Coded by Takuji Nishimura, considering the suggestions by */ +/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */ + +/* This library is free software; you can redistribute it and/or */ +/* modify it under the terms of the GNU Library General Public */ +/* License as published by the Free Software Foundation; either */ +/* version 2 of the License, or (at your option) any later */ +/* version. */ +/* This library is distributed in the hope that it will be useful, */ +/* but WITHOUT ANY WARRANTY; without even the implied warranty of */ +/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */ +/* See the GNU Library General Public License for more details. */ +/* You should have received a copy of the GNU Library General */ +/* Public License along with this library; if not, write to the */ +/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */ +/* 02111-1307 USA */ + +/* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */ +/* Any feedback is very welcome. For any question, comments, */ +/* see http://www.math.keio.ac.jp/matumoto/emt.html or email */ +/* matumoto@math.keio.ac.jp */ diff --git a/RandTest/Monkey.cpp b/RandTest/Monkey.cpp new file mode 100755 index 0000000..b4ecd3a --- /dev/null +++ b/RandTest/Monkey.cpp @@ -0,0 +1,102 @@ +//#include "StdAfx.h" +#include "Monkey.h" +#include "math.h" + +CMonkey::CMonkey(void) +{ + CreateChi2Test(); + + pow2_32 = pow( 2.0, 32 ); + ResetAll(); +} + +CMonkey::~CMonkey(void) +{ + delete MetaChi2; +} + + +// Submits a 32 bit word for OQSO testing +void CMonkey::InsertWord32(uint32_t InWord32) +{ + WordStream >>= 32; + WordStream |= (uint64_t)InWord32 << 32; + + // Submit as many overlapped words as possible from wordstream + while ( LetterPointer >= 0 ) + { + // Word is 20 bits long, overlapped every 5 bits + CurrentWord = (uint32_t)(WordStream>>LetterPointer); + CurrentWord &= 0x000fffff; + + if ( !(MonkeyBitmap.CheckWord( CurrentWord )) ) + MissingWords--; + + // Test 2097152 words (or 10485775 bits) + WordCount++; + if ( WordCount>=2097152 ) + { + // Calc current and cumulative z-scores + totalBlockCount++; + + //double UnitZScore = -((double)MissingWords-141909.1945)/294.656; + + MetaChi2->Insert(MissingWords); + + //ZScoreTotal += UnitZScore; + MissingWordsTotal += (uint64_t)MissingWords; + + P_Chi2 = MetaChi2->GetPvalue(); + + if (totalBlockCount!=0) { + //cumulativeZScore = ZScoreTotal / sqrt(totalBlockCount); + cumulativeZScore = -((double)MissingWordsTotal - (totalBlockCount * 141909.104)) / (sqrt(totalBlockCount) * 294.656); + } + + ResetTest(); + } + + LetterPointer -= 5; + } + + LetterPointer += 32; +} + +// Resets cumulative and current test +void CMonkey::ResetAll(void) +{ + cumulativeZScore = 0; + WordStream = 0; + ZScoreTotal = 0; + MissingWordsTotal = 0; + totalBlockCount = 0.; + LetterPointer = 12; + P_Chi2 = .5; + MetaChi2->Reset(); + + ResetTest(); +} + +// Resets current test +void CMonkey::ResetTest(void) +{ + MissingWords = 1048576; + MonkeyBitmap.Clearmap(); + WordCount = 0; +} + +void CMonkey::CreateChi2Test(void) { + double pTable[] = { + 0.020035263116770285, 0.04007592896841139, 0.059900259673160094, 0.07990555454845777, 0.09973585094717441, 0.12019013336507156, 0.14037340713890295, 0.16024040737242062, 0.18006661762415743, 0.20040959652731188, 0.22005106263679536, 0.2407529898402695, 0.26025548231499596, 0.2805415527829255, 0.300377215413428, 0.3208231113598671, 0.3405774462690355, 0.36077911411228614, 0.38137696414272143, 0.4009988146079183, 0.4208736799354211, 0.4409523283066323, 0.4611839548890032, 0.4815165669138008, 0.5018973841149283, 0.522273249393845, 0.5425910443898824, 0.5614554567534166, 0.5815124046447061, 0.6026740748901427, 0.6222435463400903, 0.6427760999481655, 0.6616581917860309, 0.6825745983586168, 0.7017433191528237, 0.722654551841068, 0.7428220165027213, 0.7621991815446298, 0.7817526213870567, 0.8022460422705294, 0.8224190682008065, 0.8420653268505389, 0.8617367254552067, 0.8817041720232961, 0.9019201354926609, 0.9214978872400847, 0.9412230272875216, 0.9610265302820049, 0.981060464990248, 1.0 + }; + + double boundaryTable[] = { + 141306, 141395, 141452, 141496, 141532, 141564, 141592, 141617, 141640, 141662, + 141682, 141702, 141720, 141738, 141755, 141772, 141788, 141804, 141820, 141835, + 141850, 141865, 141880, 141895, 141910, 141925, 141940, 141954, 141969, 141985, + 142000, 142016, 142031, 142048, 142064, 142082, 142100, 142118, 142137, 142158, + 142180, 142203, 142228, 142256, 142288, 142324, 142368, 142426, 142518, 1e100 + }; + + MetaChi2 = new Chi2(50, pTable, true, boundaryTable, false); +} diff --git a/RandTest/Monkey.h b/RandTest/Monkey.h new file mode 100755 index 0000000..fb9abb9 --- /dev/null +++ b/RandTest/Monkey.h @@ -0,0 +1,70 @@ +// OQSO - Overlapping-Quadruples-Sparse-Occupancy test +// Expected values are mean = 141909.47 and standard deviation = 295 +// (G. Marsaglia and A. Zaman, Computers Math. Applic., +// Vol. 26, No. 9, pp 1-10, 1993) + +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "MonkeyBitmap.h" +#include "Chi2.hpp" +#include "Gamma.h" +#include "Stat.h" + +class CMonkey +{ +public: + CMonkey(void); + ~CMonkey(void); + + // Cumulative z-score + double cumulativeZScore; + + // Total number of unit z-scores + double totalBlockCount; + + // Cumulative chi^2 results + double P_Chi2; + + // Submits a 32 bit word for OQSO testing + void InsertWord32(uint32_t InWord32); + // Resets cumulative and current test + void ResetAll(void); + +protected: + // Resets current test + void ResetTest(void); + + void CreateChi2Test(); + + // Keeps track of missing words + CMonkeyBitmap MonkeyBitmap; + + // Incomplete Gamma function + CGamma Gamma; + // Chi2 Test + Chi2* MetaChi2; + // Adds ZtoP transformation + CStat Stat; + + // Missing word count + uint32_t MissingWords; + // Words (20 bit overlapped) tested + uint32_t WordCount; + + // Actual wordstream + uint64_t WordStream; + // Wordstream 5-bit letter pointer + int LetterPointer; + // Current 20-bit word under investigation + uint32_t CurrentWord; + + // Sum of unit z-scores + double ZScoreTotal; + uint64_t MissingWordsTotal; + + // Stored value of 2^32 + double pow2_32; +}; diff --git a/RandTest/MonkeyBitmap.cpp b/RandTest/MonkeyBitmap.cpp new file mode 100755 index 0000000..fb1d1ec --- /dev/null +++ b/RandTest/MonkeyBitmap.cpp @@ -0,0 +1,40 @@ +//#include "StdAfx.h" +#include "MonkeyBitmap.h" +#include + +CMonkeyBitmap::CMonkeyBitmap(void) +{ + Clearmap(); +} + +CMonkeyBitmap::~CMonkeyBitmap(void) +{ +} + +// Sets the entire bit memory map to 0's +void CMonkeyBitmap::Clearmap(void) +{ + memset( Word, 0, 131072 ); +} + +// Check in Bitmap if this word has already been tested, true if tested previously, false otherwise; +bool CMonkeyBitmap::CheckWord(uint32_t Word20Bit) +{ + bool bRet = false; + + uint32_t Index; + uint32_t BitMask; + + // Find the index in an array of 32 bit words + Index = Word20Bit / 32; + // Now find the slot of the remainder + BitMask = 1<<(Word20Bit%32); + + // Check if bit-slot already filled + if ( BitMask & Word[Index] ) + bRet = true; + // Fill bit-slot for future + Word[Index] |= BitMask; + + return bRet; +} diff --git a/RandTest/MonkeyBitmap.h b/RandTest/MonkeyBitmap.h new file mode 100755 index 0000000..0d0af11 --- /dev/null +++ b/RandTest/MonkeyBitmap.h @@ -0,0 +1,21 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) + +class CMonkeyBitmap +{ +public: + CMonkeyBitmap(void); + ~CMonkeyBitmap(void); + + // Sets the entire bit memory map to 0's + void Clearmap(void); + // Check in Bitmap if this word has already been used, true if tested previously, false otherwise; + bool CheckWord(uint32_t Word20Bit); + +protected: + // Bitmap + uint32_t Word[32768]; +}; diff --git a/RandTest/QuickLfsrCorrector.hpp b/RandTest/QuickLfsrCorrector.hpp new file mode 100755 index 0000000..a81ad91 --- /dev/null +++ b/RandTest/QuickLfsrCorrector.hpp @@ -0,0 +1,49 @@ +#pragma once + +template +class QuickLfsrCorrector { +public: + QuickLfsrCorrector() { + T initVal = -1; + T deltaVal = initVal / 11; + for (int i=0; i<15; i++) { + lfsr[i] = initVal; + initVal -= deltaVal; + } + + inPointer = 0; + tab3 = initVal; + initVal -= deltaVal; + tab5 = initVal; + initVal -= deltaVal; + tab7 = initVal; + } + + T Correct(T inVal) { + inVal ^= tab3 ^ tab5 ^ tab7; + + inPointer--; + if (inPointer<0) + inPointer = 14; + lfsr[inPointer] = inVal; + + int tabPointer = inPointer + 3; + tabPointer %= 15; + tab3 ^= lfsr[tabPointer]; + tabPointer += 5; + tabPointer %= 15; + tab5 ^= lfsr[tabPointer]; + tabPointer += 7; + tabPointer %= 15; + tab7 ^= lfsr[tabPointer]; + + return inVal; + } + +private: + T lfsr[15]; + T tab3; + T tab5; + T tab7; + int inPointer; +}; \ No newline at end of file diff --git a/RandTest/RandTest.h b/RandTest/RandTest.h new file mode 100755 index 0000000..dcd88ec --- /dev/null +++ b/RandTest/RandTest.h @@ -0,0 +1,21 @@ +#include "Windows.h" + +#include "Bias.h" +#include "AutoCorrelation.h" +#include "BiasAndAC.h" +#include "BitCount.h" + +#include "Monkey.h" +#include "MonkeyBitMap.h" +#include "Entropy.h" +#include "Serial.h" + +#include "Chi2.h" +#include "BinomialChi2.h" +#include "ACBinomialChi2.h" + +#include "Stat.h" +#include "Gamma.h" + + + diff --git a/RandTest/Serial.cpp b/RandTest/Serial.cpp new file mode 100755 index 0000000..9aa181a --- /dev/null +++ b/RandTest/Serial.cpp @@ -0,0 +1,258 @@ +//#include "StdAfx.h" +#include "Serial.h" +#include "math.h" + +//CRITICAL_SECTION cs; +//FILE* pFile; +//char filename[200]; + +CSerial::CSerial(void) +{ +/*InitializeCriticalSection(&cs); +EnterCriticalSection(&cs); +unsigned short num = GetTickCount(); +sprintf(filename, "STable_%i.txt", num); +LeaveCriticalSection(&cs);*/ + +for (int i=0; i<50; i++) + CreateChi2Test(); + ResetAll(); + P_Chi2 = 0.5; + cumulativeSerialChi2 = 127.3339; +} + +CSerial::~CSerial(void) +{ + delete MetaChi2; +} + +// Inserts a 32 bit word into the serial test +void CSerial::InsertWord32(uint32_t InWord32) +{ + wordStream >>= 32; + wordStream |= (uint64_t)InWord32 << 32; + + // Shift through each bit to create a sub-word + while ( bitPointer >= 1 ) + { + // Take an overlapping 8 and 7 bit word from stream + uint32_t Bit8Word = ((uint32_t)(wordStream>>bitPointer)) & 0x000000ff; + uint32_t Bit7Word = Bit8Word>>1; + + // Stuff these words into independent bins for this test block + Bin8[Bit8Word]++; + Bin7[Bit7Word]++; + + // Every 8192 words binned do a block and cumulative calc + if ( (++blockWordCount)>=(16*262144) ) + { + // Calc chi^2 value for this block + double BlockSerialChi2 = BlockSumBinsSquared8() - BlockSumBinsSquared7(); + + MetaChi2->Insert(BlockSerialChi2); + + // Insert into meta chi^2 test + P_Chi2 = MetaChi2->GetPvalue(); + + // Calc cumulative chi^2 test + totalBlockCount++; + + cumulativeSerialChi2 = QuickCompare8() - QuickCompare7(); + ResetTest(); + } + + bitPointer --; + } + + bitPointer = 32; +} + +// Resets cumulative testing +void CSerial::ResetAll(void) +{ + wordStream = 0; + bitPointer = 24; + cumulativeSerialChi2 = 127.3339; + P_Chi2 = 0.5; + totalBlockCount = 0; +// ZeroMemory(CumulativeBin8, 256*sizeof(double)); + memset(CumulativeBin8, 0, 256*sizeof(double)); +// ZeroMemory(CumulativeBin7, 128*sizeof(double)); + memset(CumulativeBin7, 0, 128*sizeof(double)); + MetaChi2->Reset(); + ResetTest(); +} + +// Resets current block test +void CSerial::ResetTest(void) +{ +// ZeroMemory(Bin8, 256*sizeof(double)); + memset(Bin8, 0, 256*sizeof(double)); +// ZeroMemory(Bin7, 128*sizeof(double)); + memset(Bin7, 0, 128*sizeof(double)); + blockWordCount = 0; +} + +// Returns sum of the bis squared for the block of 8bit words +double CSerial::BlockSumBinsSquared8(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<256; i++) + { + val = Bin8[i] - ((blockWordCount/256.)); + Sum += val*val; + } + + Sum /= ((blockWordCount/256.)); + + return Sum; +} + +// Returns sum of the bis squared for the block of 7bit words +double CSerial::BlockSumBinsSquared7(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<128; i++) + { + val = Bin7[i] - ((blockWordCount/128.)); + Sum += val*val; + } + + Sum /= ((blockWordCount/128.)); + + return Sum; +} + +// Adds new block to cumulative Squared 8bit words +double CSerial::CumulativeBins8(void) +{ + double MSBS8 = 0.; + + for (int i=0; i<256; i++) + { + CumulativeBin8[i] += Bin8[i]; + MSBS8 += CumulativeBin8[i]*CumulativeBin8[i]; + } + + return MSBS8; +} + +// Adds new block to cumulative Squared 7bit words +double CSerial::CumulativeBins7(void) +{ + double MSBS7 = 0.; + + for (int i=0; i<128; i++) + { + CumulativeBin7[i] += Bin7[i]; + MSBS7 += CumulativeBin7[i]*CumulativeBin7[i]; + } + + return MSBS7; +} + +double CSerial::QuickCompare8(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<256; i++) + { + CumulativeBin8[i] += Bin8[i]; + val = (CumulativeBin8[i] / (totalBlockCount * (blockWordCount/256.))) - 2.; + Sum += val * CumulativeBin8[i]; + } + Sum += (blockWordCount*totalBlockCount); + + Sum = 0; + + for (int i=0; i<256; i++) + { + val = CumulativeBin8[i] - ((blockWordCount/256.)*totalBlockCount); + Sum += val*val; + } + + Sum /= ((blockWordCount/256.)*totalBlockCount); + + return Sum; +} + +double CSerial::QuickCompare7(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<128; i++) + { + CumulativeBin7[i] += Bin7[i]; + val = (CumulativeBin7[i] / (totalBlockCount * (blockWordCount/128.))) - 2.; + Sum += val * CumulativeBin7[i]; + } + Sum += (blockWordCount*totalBlockCount); + + Sum = 0; + + for (int i=0; i<128; i++) + { + val = CumulativeBin7[i] - ((blockWordCount/128.)*totalBlockCount); + Sum += val*val; + } + + Sum /= ((blockWordCount/128.)*totalBlockCount); + + return Sum; +} + +double CSerial::QuickCompare88(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<256; i++) + { + val = (CumulativeBin8[i] / (totalBlockCount * (blockWordCount/256.))) - 2.; + Sum += val * CumulativeBin8[i]; + } + Sum += (blockWordCount*totalBlockCount); + + return Sum; +} + +double CSerial::QuickCompare77(void) +{ + double Sum = 0; + + double val; + for (int i=0; i<128; i++) + { + val = (CumulativeBin7[i] / (totalBlockCount * (blockWordCount/128.))) - 2.; + Sum += val * CumulativeBin7[i]; + } + Sum += (blockWordCount*totalBlockCount); + + + return Sum; +} + +void CSerial::CreateChi2Test(void) { + double pTable[] = {0.019993076,0.019996996,0.020004305,0.020003831,0.019996602,0.020006223,0.020008612,0.019993278,0.020002222,0.020011793,0.019995081,0.020002287,0.019998775,0.020004405,0.019999817,0.020010438,0.019994155,0.020003369,0.020008374,0.019995879,0.020008228,0.019993876,0.020004749,0.020003208,0.019893452,0.020001708,0.020002342,0.020000863,0.020008678,0.020005665,0.020001275,0.02000936,0.020000663,0.019986486,0.020003196,0.020004123,0.020002657,0.020016139,0.019989378,0.020004423,0.019999113,0.019998812,0.020002742,0.020004996,0.020003909,0.019997932,0.020000403,0.020004695,0.019999863,0.020017545}; + + double boundryTable[] = { + 162.960449218750, 157.338134765625, 153.770507812500, 151.080566406250, 148.885253906250, + 147.007324218750, 145.352050781250, 143.862548828125, 142.498779296875, 141.234375000000, + 140.051757812500, 138.935058593750, 137.873779296875, 136.858642578125, 135.883056640625, + 134.940429687500, 134.027099609375, 133.137695312500, 132.268798828125, 131.417968750000, + 130.581298828125, 129.757324218750, 128.942626953125, 128.135498046875, 127.337646484375, + 126.538574218750, 125.740966796875, 124.942871093750, 124.141845703125, 123.336181640625, + 122.523681640625, 121.701416015625, 120.867431640625, 120.019042968750, 119.151611328125, + 118.261962890625, 117.345703125000, 116.396728515625, 115.410400390625, 114.376464843750, + 113.285400390625, 112.123046875000, 110.870117187500, 109.499511718750, 107.970458984375, + 106.216796875000, 104.117187500000, 101.411621093750, 97.322021484375, 0.0 + }; + + MetaChi2 = new Chi2(50, pTable, false, boundryTable, true); +} diff --git a/RandTest/Serial.h b/RandTest/Serial.h new file mode 100755 index 0000000..e792894 --- /dev/null +++ b/RandTest/Serial.h @@ -0,0 +1,73 @@ +#pragma once + +#pragma warning( disable : 4005 ) +#include +#pragma warning( default : 4005 ) +#include "Gamma.h" +#include "Chi2.hpp" +#include "stdio.h" + +class CSerial +{ +public: + CSerial(void); + ~CSerial(void); + // Inserts a 32 bit word into the serial test + void InsertWord32(uint32_t InWord32); + + // Resets cumulative testing + void ResetAll(void); + // Resets current block test + void ResetTest(void); + + // Serial test cumulative result; + double cumulativeSerialChi2; + + // Keeps track of all block in cumulative results + double totalBlockCount; + + // Cumulative chi^2 results + double P_Chi2; + +protected: + // Wordstream + uint64_t wordStream; + // Pointer to current bit in wordStream + int bitPointer; + + // Counts the number of 8/7 bit words binned in this block + uint32_t blockWordCount; + + // Binning for 8-bit words + double Bin8[256]; + // Binning for 7-bit words + double Bin7[128]; + + // Cumulative binning for 8-bit words + double CumulativeBin8[256]; + // Cumulative binning for 7-bit words + double CumulativeBin7[128]; + + // Returns sum of the bis squared for the block 8bit words + double BlockSumBinsSquared8(void); + // Returns sum of the bis squared for the block 7bit words + double BlockSumBinsSquared7(void); + + // Meta chi^2 test + Chi2* MetaChi2; + void CreateChi2Test(); + + // Incomplete gamma function + CGamma Gamma; + + // Adds new block to cumulative bins 8 bit words + double CumulativeBins8(void); + + // Adds new block to cumulative bins 7 bit words + double CumulativeBins7(void); + + double QuickCompare8(void); + double QuickCompare7(void); + double QuickCompare88(void); + double QuickCompare77(void); +}; diff --git a/RandTest/Stat.cpp b/RandTest/Stat.cpp new file mode 100755 index 0000000..76742e3 --- /dev/null +++ b/RandTest/Stat.cpp @@ -0,0 +1,53 @@ +//#include "StdAfx.h" +#include "Stat.h" +#include "math.h" + +CStat::CStat(void) +{ +} + +CStat::~CStat(void) +{ +} + +// Input z-score, returns cumulative normal distribution value +// Accuracy better than 1% to z=+/-7.5; .05% to z=+/-4. +double CStat::ZtoP( double zScore ) +{ + double retval; + + // calculation variables + double w; + double y; + double t; + double num; + double denom; + + // check + if ( zScore > 8. ) + { + zScore = 8.; + } + else + { + if ( zScore < -8. ) + zScore = -8.; + } + + + // calculation constants + double c[8]; + c[1] = 2.506628275; c[2] = 0.31938153; c[3] = -0.356563782; c[4] = 1.781477937; + c[5] = -1.821255978; c[6] = 1.330274429; c[7] = 0.2316419; + + w = (zScore>=0)? 1 : -1; + t = 1. + (c[7]*w*zScore); + y = 1./t; + + num = c[2] + (c[6] + (c[5]*t) + (c[4]*t*t) + (c[3]*t*t*t)) / (t*t*t*t) ; + denom = c[1] * exp( .5*zScore*zScore ) * t; + + retval = 0.5 + w * ( .5 - (num/denom) ); + + return retval; +} diff --git a/RandTest/Stat.h b/RandTest/Stat.h new file mode 100755 index 0000000..808ff4f --- /dev/null +++ b/RandTest/Stat.h @@ -0,0 +1,11 @@ +#pragma once + +class CStat +{ +public: + CStat(void); + ~CStat(void); + + // Calculate a p-value from a z-score + static double ZtoP(double zScore); +}; diff --git a/RandTest/Well44497.h b/RandTest/Well44497.h new file mode 100755 index 0000000..f98202e --- /dev/null +++ b/RandTest/Well44497.h @@ -0,0 +1,168 @@ +/* ***************************************************************************** */ +/* Copyright: Francois Panneton and Pierre L'Ecuyer, University of Montreal */ +/* Makoto Matsumoto, Hiroshima University */ +/* Notice: This code can be used freely for personal, academic, */ +/* or non-commercial purposes. For commercial purposes, */ +/* please contact P. L'Ecuyer at: lecuyer@iro.UMontreal.ca */ +/* A modified "maximally equidistributed" implementations */ +/* by Shin Harase, Hiroshima University. */ +/* ***************************************************************************** */ + +#define W 32 +#define R 1391 +#define DISCARD 15 +#define MASKU (0xffffffffU>>(W-DISCARD)) +#define MASKL (~MASKU) + +#define M1 23 +#define M2 481 +#define M3 229 + +#define MAT0POS(t,v) (v^(v>>t)) +#define MAT0NEG(t,v) (v^(v<<(-(t)))) +#define MAT1(v) v +#define MAT2(a,v) ((v & 1U)?((v>>1)^a):(v>>1)) +#define MAT3POS(t,v) (v>>t) +#define MAT3NEG(t,v) (v<<(-(t))) +#define MAT4POS(t,b,v) (v ^ ((v>> t ) & b)) +#define MAT4NEG(t,b,v) (v ^ ((v<<(-(t))) & b)) +#define MAT5(r,a,ds,dt,v) ((v & dt)?((((v<>(W-r)))&ds)^a):(((v<>(W-r)))&ds)) +#define MAT7(v) 0 + +#define V0 STATE[state_i] +#define VM1Over STATE[state_i+M1-R] +#define VM1 STATE[state_i+M1] +#define VM2Over STATE[state_i+M2-R] +#define VM2 STATE[state_i+M2] +#define VM3Over STATE[state_i+M3-R] +#define VM3 STATE[state_i+M3] +#define Vrm1 STATE[state_i-1] +#define Vrm1Under STATE[state_i+R-1] +#define Vrm2 STATE[state_i-2] +#define Vrm2Under STATE[state_i+R-2] + +#define newV0 STATE[state_i-1] +#define newV0Under STATE[state_i-1+R] +#define newV1 STATE[state_i] +#define newVRm1 STATE[state_i-2] +#define newVRm1Under STATE[state_i-2+R] + +/*output transformation parameter*/ +#define newVM2Over STATE[state_i+M2-R+1] +#define newVM2 STATE[state_i+M2+1] +#define BITMASK 0x48000000 + +static unsigned int STATE[R]; +static unsigned int z0,z1,z2; +static int state_i=0; + +static unsigned int case_1(void); +static unsigned int case_2(void); +static unsigned int case_3(void); +static unsigned int case_4(void); +static unsigned int case_5(void); +static unsigned int case_6(void); + +unsigned int (*WELLRNG44497)(void); + +void SeedWELL(unsigned int seed) { + state_i=0; + WELLRNG44497 = case_1; + + if (seed == 0U) + seed = 5489U; + + STATE[0] = seed & 0xffffffffUL; + + // Same generator used to seed Mersenne twister + for (int i=1; i> 30)) + i); + + // mix it up to avoid bias + for (int i=0; i<10000*R; i++) + WELLRNG44497(); +} + +void InitWELLRNG44497(unsigned int *init ){ + int j; + state_i=0; + WELLRNG44497 = case_1; + for(j=0;j= R + z0 = (Vrm1 & MASKL) | (Vrm2 & MASKU); + z1 = MAT0NEG(-24,V0) ^ MAT0POS(30,VM1Over); + z2 = MAT0NEG(-10,VM2Over) ^ MAT3NEG(-26,VM3Over); + newV1 = z1 ^ z2; + newV0 = MAT1(z0) ^ MAT0POS(20,z1) ^ MAT5(9,0xb729fcecU,0xfbffffffU,0x00020000U,z2) ^ MAT1(newV1); + state_i--; + if(state_i+M1= R + z0 = (Vrm1 & MASKL) | (Vrm2 & MASKU); + z1 = MAT0NEG(-24,V0) ^ MAT0POS(30,VM1); + z2 = MAT0NEG(-10,VM2Over) ^ MAT3NEG(-26,VM3Over); + newV1 = z1 ^ z2; + newV0 = MAT1(z0) ^ MAT0POS(20,z1) ^ MAT5(9,0xb729fcecU,0xfbffffffU,0x00020000U,z2) ^ MAT1(newV1); + state_i--; + if (state_i+M3 < R) + WELLRNG44497 = case_5; + return (STATE[state_i] ^ (newVM2Over & BITMASK)); +} + +static unsigned int case_5(void){ + //state_i+M2 >= R + z0 = (Vrm1 & MASKL) | (Vrm2 & MASKU); + z1 = MAT0NEG(-24,V0) ^ MAT0POS(30,VM1); + z2 = MAT0NEG(-10,VM2Over) ^ MAT3NEG(-26,VM3); + newV1 = z1 ^ z2; + newV0 = MAT1(z0) ^ MAT0POS(20,z1) ^ MAT5(9,0xb729fcecU,0xfbffffffU,0x00020000U,z2) ^ MAT1(newV1); + state_i--; + if(state_i+M2 < R) + WELLRNG44497 = case_6; + return (STATE[state_i] ^ (newVM2Over & BITMASK)); +} + +static unsigned int case_6(void){ + // 2 <= state_i <= R-M2-1 + z0 = (Vrm1 & MASKL) | (Vrm2 & MASKU); + z1 = MAT0NEG(-24,V0) ^ MAT0POS(30,VM1); + z2 = MAT0NEG(-10,VM2) ^ MAT3NEG(-26,VM3); + newV1 = z1 ^ z2; + newV0 = MAT1(z0) ^ MAT0POS(20,z1) ^ MAT5(9,0xb729fcecU,0xfbffffffU,0x00020000U,z2) ^ MAT1(newV1); + state_i--; + if(state_i == 1 ) + WELLRNG44497 = case_2; + return (STATE[state_i] ^ (newVM2 & BITMASK)); +} diff --git a/RandTest/stdafx.cpp b/RandTest/stdafx.cpp new file mode 100755 index 0000000..dd1aedb --- /dev/null +++ b/RandTest/stdafx.cpp @@ -0,0 +1,8 @@ +// stdafx.cpp : source file that includes just the standard includes +// RandTest.pch will be the pre-compiled header +// stdafx.obj will contain the pre-compiled type information + +#include "stdafx.h" + +// TODO: reference any additional headers you need in STDAFX.H +// and not in this file diff --git a/RandTest/stdafx.h b/RandTest/stdafx.h new file mode 100755 index 0000000..da438e7 --- /dev/null +++ b/RandTest/stdafx.h @@ -0,0 +1,11 @@ +// stdafx.h : include file for standard system include files, +// or project specific include files that are used frequently, but +// are changed infrequently +// + +#pragma once + +#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers + +// TODO: reference additional headers your program requires here +#include "RandTest.h"