/**
 *  @file   test2.cpp
 *  @brief  C言語によるアルゴリズム辞典、のrndtest.c による測定
 *  @author tenk*
 *  @date   2010-1
 */

#define XORSHIFTRAND_TEST
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "XorShiftRand.h"
#include "MtRand.h"

// ===========================================================================

/// 適当な線形合同法の乱数. (http://www001.upp.so-net.ne.jp/isaku/rand.html の rand( )その１)
class Test_Rand1 {
    unsigned    s_;
public:
    Test_Rand1() : s_(1) {}

    double operator()() {
        return (s_ = s_ * 1103515245+12345) / (0xFFFFFFFF + 1.0);
    }
};


/// C標準ライブラリの rand() を呼び出す.
struct Test_Clib_Rand {
    double operator()() { return rand() / (RAND_MAX+1.0) ; }
};


template<class R>
struct Test_Rnd : public R {
    double operator()() { return R::get() / (0xFFFFFFFF + 1.0); }
};



// ===========================================================================
// 計測＆チェック

template<class RND>
void test(RND rnd, size_t count, unsigned flags, const char* ttl, double& r1, double& r2)
{
    unsigned char   verboss     = flags & 1;
    if (verboss) {
        fprintf(stderr, "\t%s:start\n", ttl);
    }
    double x0    = rnd() - 0.5;
    double xprev = x0;
    double s1    = x0;
    double s2    = x0 * x0;
    double r     = 0;
    double x     = 0;
    unsigned  n  = count;
    for (unsigned i = 1; i < n; ++i) {
        x     = rnd() - 0.5;
        s1   += x;
        s2   += x * x;
        r    += xprev * x;
        xprev = x;
    }
    r  = (n * (r + x * x0) - s1 * s1) /(n * s2 - s1 * s1);
    r1 = s1 * sqrt( 12.0 / n );
    r2 = ((n-1) * r + 1) * sqrt((n + 1.0) / (n * (n - 3.0)));
    return;
}


static void test_all(unsigned count, unsigned flags)
{
    enum { N = 64 };
    const char* nm[N]       = {0};
    double      r1[N] = {0}, r2[N] = {0};

    printf("C言語によるアルゴリズム辞典、のrndtest.c による測定\n", count);
    printf(" -0.5〜0.5の乱数を %9d 回実行\n", count);

    // いくつかのテスト呼び出し.
    {
        unsigned k     = 0;
        nm[k] = "XorShift  32"  ; test(Test_Rnd<XorShiftRand<32> >() , count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift  64"  ; test(Test_Rnd<XorShiftRand<64> >() , count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift  96"  ; test(Test_Rnd<XorShiftRand<96> >() , count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift 128"  ; test(Test_Rnd<XorShiftRand<128> >(), count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift 160"  ; test(Test_Rnd<XorShiftRand<160> >(), count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift 192"  ; test(Test_Rnd<XorShiftRand<192> >(), count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "XorShift 256"  ; test(Test_Rnd<XorShiftRand<256> >(), count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "LCG         "  ; test(Test_Rand1()                  , count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "clib rand   "  ; test(Test_Clib_Rand()              , count, flags, nm[k], r1[k], r2[k]); ++k;
        nm[k] = "MT Rand     "  ; test(Test_Rnd<MtRand>()            , count, flags, nm[k], r1[k], r2[k]); ++k;
    }

    printf(" 以下は期待値との差を標準誤差で割ったもの.\n");
    printf(" 20回中19回は ±2以内に入るはず.\n");
    printf("[1] 平均値              \n");
    printf("[2] 隣どうしの相関係数:  \n");
    printf(",%-21s ,%6s,%6s\n", "方法", "[1]", "[2]");
    for (unsigned j = 0; j < N; ++j) {
        if (nm[j] == 0)
            continue;
        printf(",%-21s ,%6.3f,%6.3f\n", nm[j], r1[j], r2[j]);
    }
}



// ===========================================================================

/** 使い方.
 */
static int usage()
{
    fprintf(stderr, "usage>test [-opts] file(s)\n"
           "  xor shift 乱数チェック.\n"
           "  -cN   試行回数の指定\n"
           "  -v    経過メッセージを増やす\n"
          );
    return 1;
}



/** メイン.
 */
int main(int argc, const char* argv[])
{
    int      i;
    unsigned count = 1;
    unsigned flags = 0;

    if (argc < 2)
        return usage();

    // オプション取得.
    for (i = 1; i < argc; ++i) {
        const char* a = argv[i];
        if (a[0] == '-') {
            if (a[1] == 'v') {
                flags = 1;      // verboss;
            } else if (a[1] == 'c') {
                count   = strtol(a+2,0,0);
                if (count < 1)
                    count = 1;
            } else {
                return usage();
            }
        }
    }

    test_all(count, flags);

    return 0;
}