/*10:*/
#line 628 "./wiener.w"

/*11:*/
#line 653 "./wiener.w"

#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 
#include <ctype.h> 

/*:11*/
#line 629 "./wiener.w"

/*12:*/
#line 679 "./wiener.w"

#define VERSION "1.0"
#define COPYRIGHT "Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>"
#define SUCCESS (0)
#define FAILURE (1)
#define QUALITY (1009)
#define KK (100)
#define LL  (37)
#define DEFAULT_SEED (310952)
#define cmd_match(s,l,c) ((!strcmp((c),(s)))||(!strcmp((c),(l))))
#define MODE_WIENER_PROCESS (0)
#define MODE_LOCKED_UNIFORM_DISTRIBUTION (1)
#define MODE_LOCKED_NORMAL_DISTRIBUTION (2)

/*:12*/
#line 630 "./wiener.w"

/*13:*/
#line 701 "./wiener.w"

extern char*optarg;
char*progname;
double ranf_arr_buf[QUALITY];
double ranf_arr_dummy= -1.0,ranf_arr_started= -1.0;
double*ranf_arr_ptr= &ranf_arr_dummy;

/*:13*/
#line 631 "./wiener.w"

/*14:*/
#line 712 "./wiener.w"

/*15:*/
#line 760 "./wiener.w"

#define mod_sum(x,y) (((x)+(y))-(int)((x)+(y)))
double ran_u[KK];

void ranf_array(double aa[],int n)
{
register int i,j;

for(j= 0;j<KK;j++)aa[j]= ran_u[j];
for(;j<n;j++)aa[j]= mod_sum(aa[j-KK],aa[j-LL]);
for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK],aa[j-LL]);
for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK],ran_u[i-LL]);
}

void ranf_matrix(double**aa,int mm,int dd)
{
register int i,j,col;

for(col= 0;col<dd;col++){
for(j= 0;j<KK;j++)aa[j][col]= ran_u[j];
for(;j<mm;j++)aa[j][col]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],ran_u[i-LL]);
}
}

/*:15*/
#line 713 "./wiener.w"

/*16:*/
#line 798 "./wiener.w"

#define TT  70   
#define is_odd(s) ((s)&1)

void ranf_start(long seed){
register int t,s,j;
double u[KK+KK-1];
double ulp= (1.0/(1L<<30))/(1L<<22);
double ss= 2.0*ulp*((seed&0x3fffffff)+2);

for(j= 0;j<KK;j++){
u[j]= ss;
ss+= ss;
if(ss>=1.0)ss-= 1.0-2*ulp;
}
u[1]+= ulp;
for(s= seed&0x3fffffff,t= TT-1;t;){
for(j= KK-1;j> 0;j--)
u[j+j]= u[j],u[j+j-1]= 0.0;
for(j= KK+KK-2;j>=KK;j--){
u[j-(KK-LL)]= mod_sum(u[j-(KK-LL)],u[j]);
u[j-KK]= mod_sum(u[j-KK],u[j]);
}
if(is_odd(s)){
for(j= KK;j> 0;j--)u[j]= u[j-1];
u[0]= u[KK];
u[LL]= mod_sum(u[LL],u[KK]);
}
if(s)s>>= 1;else t--;
}
for(j= 0;j<LL;j++)ran_u[j+KK-LL]= u[j];
for(;j<KK;j++)ran_u[j-LL]= u[j];
for(j= 0;j<10;j++)ranf_array(u,KK+KK-1);
ranf_arr_ptr= &ranf_arr_started;
}

#define ranf_arr_next() (*ranf_arr_ptr>=0? *ranf_arr_ptr++: ranf_arr_cycle())
double ranf_arr_cycle()
{
if(ranf_arr_ptr==&ranf_arr_dummy)
ranf_start(314159L);
ranf_array(ranf_arr_buf,QUALITY);
ranf_arr_buf[KK]= -1;
ranf_arr_ptr= ranf_arr_buf+1;
return ranf_arr_buf[0];
}

/*:16*/
#line 714 "./wiener.w"

/*17:*/
#line 874 "./wiener.w"

#define PI (3.14159265358979323846264338)
void normdist(double**aa,int mm,int dd){
register int j,k;
register double f,z;

for(j= 0;j<dd;j++){
for(k= 0;k<mm-1;k+= 2){
if(aa[k][j]> 0.0){
f= sqrt(-2*log(aa[k][j]));
z= 2.0*PI*aa[k+1][j];
aa[k][j]= f*cos(z);
aa[k+1][j]= f*sin(z);
}else{
fprintf(stderr,"%s: Error: Zero element detected!\n",progname);
fprintf(stderr,"%s: (row %d, column %d)\n",progname,k,j);
}
}
}
return;
}
#undef PI

/*:17*/
#line 715 "./wiener.w"

/*18:*/
#line 942 "./wiener.w"

void wiener(double**aa,int mm,int dd,int seed,short mode)
{
register int j,k;

ranf_start(seed);
ranf_matrix(aa,mm,dd);
if(mode==MODE_LOCKED_UNIFORM_DISTRIBUTION)return;
normdist(aa,mm,dd);
if(mode==MODE_LOCKED_NORMAL_DISTRIBUTION)return;
for(j= 0;j<dd;j++){
aa[0][j]= 0.0;
for(k= 1;k<mm;k++){
aa[k][j]+= aa[k-1][j];
}
}
}

/*:18*/
#line 716 "./wiener.w"

/*19:*/
#line 965 "./wiener.w"

double**dmatrix(long nrl,long nrh,long ncl,long nch)
{
long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
double**m;
m= (double**)malloc((size_t)((nrow+1)*sizeof(double*)));
if(!m){
fprintf(stderr,"%s: Allocation failure 1 in dmatrix()\n",progname);
exit(FAILURE);
}
m+= 1;
m-= nrl;
m[nrl]= (double*)malloc((size_t)((nrow*ncol+1)*sizeof(double)));
if(!m[nrl]){
fprintf(stderr,"%s: Allocation failure 2 in dmatrix()\n",progname);
exit(FAILURE);
}
m[nrl]+= 1;
m[nrl]-= ncl;
for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
return m;
}

/*:19*/
#line 717 "./wiener.w"

/*20:*/
#line 992 "./wiener.w"

void free_dmatrix(double**m,long nrl,long nrh,long ncl,long nch){
free((char*)(m[nrl]+ncl-1));
free((char*)(m+nrl-1));
}

/*:20*/
#line 718 "./wiener.w"

/*21:*/
#line 1000 "./wiener.w"

void display_help_message(void){
fprintf(stderr,"Usage: %s M [options]\n",progname);
fprintf(stderr,"Options:\n");
fprintf(stderr," -h, --help\n");
fprintf(stderr,"    Display this help message and exit clean.\n");
fprintf(stderr," -v, --verbose\n");
fprintf(stderr,"    Toggle verbose mode. Default: off.\n");
fprintf(stderr," -M, --num_samples <M>\n");
fprintf(stderr,"    Generate M samples of data. Here M should always be\n");
fprintf(stderr,"    an even number, greater than the long lag KK=%d.\n",KK);
fprintf(stderr,"    If an odd number is specified, the program will\n");
fprintf(stderr,"    automatically adjust this to the next higher\n");
fprintf(stderr,"    integer. Default: M=%d.\n",KK);
fprintf(stderr," -D, --dimension <D>\n");
fprintf(stderr,"    Generate D-dimensional samples of data. Default: D=1.\n");
fprintf(stderr," -s, --seed <seed>\n");
fprintf(stderr,"    Define a custom seed number for the initialization\n");
fprintf(stderr,"    of the random number generator. Default: seed=%d.\n",
DEFAULT_SEED);
fprintf(stderr," -u, --uniform\n");
fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
fprintf(stderr,"    to simply generate a uniform distribution of\n");
fprintf(stderr,"    D-dimensional points, with each element distributed\n");
fprintf(stderr,"    over the interval [0,1].\n");
fprintf(stderr," -n, --normal\n");
fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
fprintf(stderr,"    to simply generate a normal distribution of\n");
fprintf(stderr,"    D-dimensional points, with each element distributed\n");
fprintf(stderr,"    with zero expectation value and unit variance.\n");
}

/*:21*/
#line 719 "./wiener.w"

/*22:*/
#line 1039 "./wiener.w"

short pathcharacter(int ch){
return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
(ch=='-')||(ch=='+'));
}

/*:22*/
#line 720 "./wiener.w"

/*23:*/
#line 1051 "./wiener.w"

char*strip_away_path(char filename[]){
int j,k= 0;
while(pathcharacter(filename[k]))k++;
j= (--k);
while(isalnum((int)(filename[j])))j--;
j++;
return(&filename[j]);
}

/*:23*/
#line 721 "./wiener.w"


/*:14*/
#line 632 "./wiener.w"


int main(int argc,char*argv[])
{
/*24:*/
#line 1092 "./wiener.w"

double**aa;
unsigned long mm= KK;
unsigned dd= 1;
int seed= DEFAULT_SEED;
short mode= MODE_WIENER_PROCESS,verbose= 0;
int no_arg;
register int j,k;

/*:24*/
#line 636 "./wiener.w"

/*25:*/
#line 1106 "./wiener.w"

progname= strip_away_path(argv[0]);
no_arg= argc;
while(--argc){
if(cmd_match("-h","--help",argv[no_arg-argc])){
display_help_message();
exit(SUCCESS);
}else if(cmd_match("-v","--verbose",argv[no_arg-argc])){
verbose= (verbose?0:1);
}else if(cmd_match("-M","--num_samples",argv[no_arg-argc])){
--argc;
if(!sscanf(argv[no_arg-argc],"%lu",&mm)){
fprintf(stderr,"%s: Error detected when parsing the number of "
"samples.\n",progname);
display_help_message();
exit(FAILURE);
}
if(mm<KK){
fprintf(stderr,"%s: The M number of data points must be at least "
"the long lag of the generator, M >= KK = %d.\n",progname,KK);
exit(FAILURE);
}
if(is_odd(mm))mm++;
}else if(cmd_match("-D","--dimension",argv[no_arg-argc])){
--argc;
if(!sscanf(argv[no_arg-argc],"%ud",&dd)){
fprintf(stderr,"%s: Error detected when parsing dimension.\n",
progname);
display_help_message();
exit(FAILURE);
}
if(dd<1){
fprintf(stderr,"%s: Dimension D should be at least 1.\n",progname);
exit(FAILURE);
}
}else if(cmd_match("-s","--seed",argv[no_arg-argc])){
--argc;
if(!sscanf(argv[no_arg-argc],"%d",&seed)){
fprintf(stderr,"%s: Error detected when parsing the seed of the "
"initializer.\n",progname);
display_help_message();
exit(FAILURE);
}
}else if(cmd_match("-u","--uniform",argv[no_arg-argc])){
mode= MODE_LOCKED_UNIFORM_DISTRIBUTION;
}else if(cmd_match("-n","--normal",argv[no_arg-argc])){
mode= MODE_LOCKED_NORMAL_DISTRIBUTION;
}else{
fprintf(stderr,"%s: Sorry, I do not recognize option '%s'.\n",
progname,argv[no_arg-argc]);
display_help_message();
exit(FAILURE);
}
}
if(verbose)
fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);

/*:25*/
#line 637 "./wiener.w"

/*26:*/
#line 1169 "./wiener.w"

aa= dmatrix(0,mm-1,0,dd-1);

/*:26*/
#line 638 "./wiener.w"

/*27:*/
#line 1174 "./wiener.w"

wiener(aa,mm,dd,seed,mode);

/*:27*/
#line 639 "./wiener.w"

/*28:*/
#line 1179 "./wiener.w"

for(k= 0;k<mm;k++){
for(j= 0;j<dd-1;j++)
printf("%.20f ",aa[k][j]);
printf("%.20f\n",aa[k][j]);
}

/*:28*/
#line 640 "./wiener.w"

/*29:*/
#line 1188 "./wiener.w"

free_dmatrix(aa,0,mm-1,0,dd-1);

/*:29*/
#line 641 "./wiener.w"

return(SUCCESS);
}

/*:10*/