/*
 *  psychlops_m_util.cpp
 *  Psychlops Standard Library (Universal)
 *
 *  Last Modified 2006/01/04 by Kenchi HOSOKAWA
 *  (C) 2006 Kenchi HOSOKAWA, Kazushi MARUYA and Takao SATO
 */


#include <stdlib.h>
#include <Math.h>
#include <time.h>


#if defined(__GNUC__)
#  include "dSFMT2.0/dSFMT.c"
#elif defined(_MSC_VER)
#  include "dSFMT2.0/dSFMT.c"
#elif defined(__BORLANDC__)
#  include "dSFMT2.0/dSFMT.c"
#endif

#include "dSFMT2.0/dSFMT.h"
#include "okumura/igamma.h"



#include "../ApplicationInterfaces/psychlops_code.h"
#include "psychlops_m_util.h"
#include "psychlops_m_matrix.h"


namespace Psychlops {


	const double PI	= 3.14159265358979;
	const double Math::E = 2.718281828459045;
	const double Math::LOG2E = 1.44269504088896340736;


	/*		Functions		*/

	double random() {
		return dsfmt_gv_genrand_close_open();
	}
	double random(double factor) {
		return dsfmt_gv_genrand_close_open()*factor;
	}
	int random(int factor) {
		return (int)(dsfmt_gv_genrand_close_open()*factor);
	}
	double random(double min, double max) {
		double dif = (max-min);
		if(dif<0) throw Exception("min exceeds max in random(min, max)");
		return random(dif)+min;
	}
	void random(double array[], int size) {
		for(int i=0; i<size; i++) array[i]=random();
	}
	void random(double array[], int size, double factor) {
		for(int i=0; i<size; i++) array[i]=random(factor);
	}
	void random(double array[], int size, double min, double max) {
		double factor = (max-min);
		if(factor<0) throw Exception("min exceeds max in random(min, max)");
		for(int i=0; i<size; i++) array[i]=random(factor)+min;
	}
	void random(Matrix &mtx) {
		double *array, *str; int step; mtx.get_element_ptr(array, str, step);
		int size = mtx.getRows()*mtx.getCols();
		for(int i=0; i<size; i++) array[i]=random();
	}
	void random(Matrix &mtx, double factor) {
		double *array, *str; int step; mtx.get_element_ptr(array, str, step);
		int size = mtx.getRows()*mtx.getCols();
		for(int i=0; i<size; i++) array[i]=random(factor);
	}
	void random(Matrix &mtx, double min, double max) {
		double *array, *str; int step; mtx.get_element_ptr(array, str, step);
		int size = mtx.getRows()*mtx.getCols();
		double factor = (max-min);
		if(factor<0) throw Exception("min exceeds max in random(min, max)");
		for(int i=0; i<size; i++) array[i]=random(factor)+min;
	}
	void randomize(unsigned long seed) {
		srand( (int)seed );
		dsfmt_gv_init_gen_rand( seed );
//		for(int i=0; i<1000000; i++) random();
	}
	void randomize() {
		srand( (int)time(NULL) );
		dsfmt_gv_init_gen_rand( (unsigned long)time(NULL) );
//		for(int i=0; i<1000000; i++) random();
	}

	void Math::setscramble(/*unsigned*/ int array[], unsigned int length) {
		unsigned int rnd;
		for(unsigned int i=0;i<length;i++) { array[i] = -1; }
		for(unsigned int i=0;i<length;i++) {
			rnd = rand()%length;
			while(array[rnd]!=-1) {
				rnd++;
				if(rnd>=length) { rnd=0; }
			}
			array[rnd] = i;
		}
	}


	int Math::round(double val) {
		double integer_part, particle = modf(val, &integer_part);
		return ( ( particle<0.5 | (particle==0.5 && (int)integer_part%2==0) ) ? (int)integer_part : (int)integer_part+1 );
	}
	double Math::log2(double val) {
		return log(val)*LOG2E;
	}
	/*
	double Math::max(double val1, double val2) {
		return val1>val2 ? val1 : val2;
	}
	double Math::min(double val1, double val2) {
		return val1<val2 ? val1 : val2;
	}
	*/


	double Math::gaussian(double x, double sigma) {
		return exp( -(x*x) / (2.0*sigma*sigma) );
	}

	double Math::normalDistibution(double x, double mu, double sigma) {
		return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);
	}

	double Math::cumulativeNormalDistibution(double x, double mu, double sigma) {
		return .5 + .5*erf( (x-mu)/(sigma*sqrt(2.0) ) );
	}

/*
	double Math::mod(double lhs, double rhs) {
		return lhs - floor(lhs/rhs)*rhs;
	}

	double Math::abs(double val) {
		return (val<0 ? -val : val);
	}
*/


}	/*	<- namespace Psycholops 	*/
