Search:

Return to previous page

Contents of file 'wiener/wiener.c':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   /*10:*/
    2   #line 628 "./wiener.w"
    3   
    4   /*11:*/
    5   #line 653 "./wiener.w"
    6   
    7   #include <math.h> 
    8   #include <stdio.h> 
    9   #include <stdlib.h> 
   10   #include <string.h> 
   11   #include <ctype.h> 
   12   
   13   /*:11*/
   14   #line 629 "./wiener.w"
   15   
   16   /*12:*/
   17   #line 679 "./wiener.w"
   18   
   19   #define VERSION "1.0"
   20   #define COPYRIGHT "Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>"
   21   #define SUCCESS (0)
   22   #define FAILURE (1)
   23   #define QUALITY (1009)
   24   #define KK (100)
   25   #define LL  (37)
   26   #define DEFAULT_SEED (310952)
   27   #define cmd_match(s,l,c) ((!strcmp((c),(s)))||(!strcmp((c),(l))))
   28   #define MODE_WIENER_PROCESS (0)
   29   #define MODE_LOCKED_UNIFORM_DISTRIBUTION (1)
   30   #define MODE_LOCKED_NORMAL_DISTRIBUTION (2)
   31   
   32   /*:12*/
   33   #line 630 "./wiener.w"
   34   
   35   /*13:*/
   36   #line 701 "./wiener.w"
   37   
   38   extern char*optarg;
   39   char*progname;
   40   double ranf_arr_buf[QUALITY];
   41   double ranf_arr_dummy= -1.0,ranf_arr_started= -1.0;
   42   double*ranf_arr_ptr= &ranf_arr_dummy;
   43   
   44   /*:13*/
   45   #line 631 "./wiener.w"
   46   
   47   /*14:*/
   48   #line 712 "./wiener.w"
   49   
   50   /*15:*/
   51   #line 760 "./wiener.w"
   52   
   53   #define mod_sum(x,y) (((x)+(y))-(int)((x)+(y)))
   54   double ran_u[KK];
   55   
   56   void ranf_array(double aa[],int n)
   57   {
   58   register int i,j;
   59   
   60   for(j= 0;j<KK;j++)aa[j]= ran_u[j];
   61   for(;j<n;j++)aa[j]= mod_sum(aa[j-KK],aa[j-LL]);
   62   for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK],aa[j-LL]);
   63   for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK],ran_u[i-LL]);
   64   }
   65   
   66   void ranf_matrix(double**aa,int mm,int dd)
   67   {
   68   register int i,j,col;
   69   
   70   for(col= 0;col<dd;col++){
   71   for(j= 0;j<KK;j++)aa[j][col]= ran_u[j];
   72   for(;j<mm;j++)aa[j][col]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
   73   for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
   74   for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],ran_u[i-LL]);
   75   }
   76   }
   77   
   78   /*:15*/
   79   #line 713 "./wiener.w"
   80   
   81   /*16:*/
   82   #line 798 "./wiener.w"
   83   
   84   #define TT  70   
   85   #define is_odd(s) ((s)&1)
   86   
   87   void ranf_start(long seed){
   88   register int t,s,j;
   89   double u[KK+KK-1];
   90   double ulp= (1.0/(1L<<30))/(1L<<22);
   91   double ss= 2.0*ulp*((seed&0x3fffffff)+2);
   92   
   93   for(j= 0;j<KK;j++){
   94   u[j]= ss;
   95   ss+= ss;
   96   if(ss>=1.0)ss-= 1.0-2*ulp;
   97   }
   98   u[1]+= ulp;
   99   for(s= seed&0x3fffffff,t= TT-1;t;){
  100   for(j= KK-1;j> 0;j--)
  101   u[j+j]= u[j],u[j+j-1]= 0.0;
  102   for(j= KK+KK-2;j>=KK;j--){
  103   u[j-(KK-LL)]= mod_sum(u[j-(KK-LL)],u[j]);
  104   u[j-KK]= mod_sum(u[j-KK],u[j]);
  105   }
  106   if(is_odd(s)){
  107   for(j= KK;j> 0;j--)u[j]= u[j-1];
  108   u[0]= u[KK];
  109   u[LL]= mod_sum(u[LL],u[KK]);
  110   }
  111   if(s)s>>= 1;else t--;
  112   }
  113   for(j= 0;j<LL;j++)ran_u[j+KK-LL]= u[j];
  114   for(;j<KK;j++)ran_u[j-LL]= u[j];
  115   for(j= 0;j<10;j++)ranf_array(u,KK+KK-1);
  116   ranf_arr_ptr= &ranf_arr_started;
  117   }
  118   
  119   #define ranf_arr_next() (*ranf_arr_ptr>=0? *ranf_arr_ptr++: ranf_arr_cycle())
  120   double ranf_arr_cycle()
  121   {
  122   if(ranf_arr_ptr==&ranf_arr_dummy)
  123   ranf_start(314159L);
  124   ranf_array(ranf_arr_buf,QUALITY);
  125   ranf_arr_buf[KK]= -1;
  126   ranf_arr_ptr= ranf_arr_buf+1;
  127   return ranf_arr_buf[0];
  128   }
  129   
  130   /*:16*/
  131   #line 714 "./wiener.w"
  132   
  133   /*17:*/
  134   #line 874 "./wiener.w"
  135   
  136   #define PI (3.14159265358979323846264338)
  137   void normdist(double**aa,int mm,int dd){
  138   register int j,k;
  139   register double f,z;
  140   
  141   for(j= 0;j<dd;j++){
  142   for(k= 0;k<mm-1;k+= 2){
  143   if(aa[k][j]> 0.0){
  144   f= sqrt(-2*log(aa[k][j]));
  145   z= 2.0*PI*aa[k+1][j];
  146   aa[k][j]= f*cos(z);
  147   aa[k+1][j]= f*sin(z);
  148   }else{
  149   fprintf(stderr,"%s: Error: Zero element detected!\n",progname);
  150   fprintf(stderr,"%s: (row %d, column %d)\n",progname,k,j);
  151   }
  152   }
  153   }
  154   return;
  155   }
  156   #undef PI
  157   
  158   /*:17*/
  159   #line 715 "./wiener.w"
  160   
  161   /*18:*/
  162   #line 942 "./wiener.w"
  163   
  164   void wiener(double**aa,int mm,int dd,int seed,short mode)
  165   {
  166   register int j,k;
  167   
  168   ranf_start(seed);
  169   ranf_matrix(aa,mm,dd);
  170   if(mode==MODE_LOCKED_UNIFORM_DISTRIBUTION)return;
  171   normdist(aa,mm,dd);
  172   if(mode==MODE_LOCKED_NORMAL_DISTRIBUTION)return;
  173   for(j= 0;j<dd;j++){
  174   aa[0][j]= 0.0;
  175   for(k= 1;k<mm;k++){
  176   aa[k][j]+= aa[k-1][j];
  177   }
  178   }
  179   }
  180   
  181   /*:18*/
  182   #line 716 "./wiener.w"
  183   
  184   /*19:*/
  185   #line 965 "./wiener.w"
  186   
  187   double**dmatrix(long nrl,long nrh,long ncl,long nch)
  188   {
  189   long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
  190   double**m;
  191   m= (double**)malloc((size_t)((nrow+1)*sizeof(double*)));
  192   if(!m){
  193   fprintf(stderr,"%s: Allocation failure 1 in dmatrix()\n",progname);
  194   exit(FAILURE);
  195   }
  196   m+= 1;
  197   m-= nrl;
  198   m[nrl]= (double*)malloc((size_t)((nrow*ncol+1)*sizeof(double)));
  199   if(!m[nrl]){
  200   fprintf(stderr,"%s: Allocation failure 2 in dmatrix()\n",progname);
  201   exit(FAILURE);
  202   }
  203   m[nrl]+= 1;
  204   m[nrl]-= ncl;
  205   for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
  206   return m;
  207   }
  208   
  209   /*:19*/
  210   #line 717 "./wiener.w"
  211   
  212   /*20:*/
  213   #line 992 "./wiener.w"
  214   
  215   void free_dmatrix(double**m,long nrl,long nrh,long ncl,long nch){
  216   free((char*)(m[nrl]+ncl-1));
  217   free((char*)(m+nrl-1));
  218   }
  219   
  220   /*:20*/
  221   #line 718 "./wiener.w"
  222   
  223   /*21:*/
  224   #line 1000 "./wiener.w"
  225   
  226   void display_help_message(void){
  227   fprintf(stderr,"Usage: %s M [options]\n",progname);
  228   fprintf(stderr,"Options:\n");
  229   fprintf(stderr," -h, --help\n");
  230   fprintf(stderr,"    Display this help message and exit clean.\n");
  231   fprintf(stderr," -v, --verbose\n");
  232   fprintf(stderr,"    Toggle verbose mode. Default: off.\n");
  233   fprintf(stderr," -M, --num_samples <M>\n");
  234   fprintf(stderr,"    Generate M samples of data. Here M should always be\n");
  235   fprintf(stderr,"    an even number, greater than the long lag KK=%d.\n",KK);
  236   fprintf(stderr,"    If an odd number is specified, the program will\n");
  237   fprintf(stderr,"    automatically adjust this to the next higher\n");
  238   fprintf(stderr,"    integer. Default: M=%d.\n",KK);
  239   fprintf(stderr," -D, --dimension <D>\n");
  240   fprintf(stderr,"    Generate D-dimensional samples of data. Default: D=1.\n");
  241   fprintf(stderr," -s, --seed <seed>\n");
  242   fprintf(stderr,"    Define a custom seed number for the initialization\n");
  243   fprintf(stderr,"    of the random number generator. Default: seed=%d.\n",
  244   DEFAULT_SEED);
  245   fprintf(stderr," -u, --uniform\n");
  246   fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
  247   fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
  248   fprintf(stderr,"    to simply generate a uniform distribution of\n");
  249   fprintf(stderr,"    D-dimensional points, with each element distributed\n");
  250   fprintf(stderr,"    over the interval [0,1].\n");
  251   fprintf(stderr," -n, --normal\n");
  252   fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
  253   fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
  254   fprintf(stderr,"    to simply generate a normal distribution of\n");
  255   fprintf(stderr,"    D-dimensional points, with each element distributed\n");
  256   fprintf(stderr,"    with zero expectation value and unit variance.\n");
  257   }
  258   
  259   /*:21*/
  260   #line 719 "./wiener.w"
  261   
  262   /*22:*/
  263   #line 1039 "./wiener.w"
  264   
  265   short pathcharacter(int ch){
  266   return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
  267   (ch=='-')||(ch=='+'));
  268   }
  269   
  270   /*:22*/
  271   #line 720 "./wiener.w"
  272   
  273   /*23:*/
  274   #line 1051 "./wiener.w"
  275   
  276   char*strip_away_path(char filename[]){
  277   int j,k= 0;
  278   while(pathcharacter(filename[k]))k++;
  279   j= (--k);
  280   while(isalnum((int)(filename[j])))j--;
  281   j++;
  282   return(&filename[j]);
  283   }
  284   
  285   /*:23*/
  286   #line 721 "./wiener.w"
  287   
  288   
  289   /*:14*/
  290   #line 632 "./wiener.w"
  291   
  292   
  293   int main(int argc,char*argv[])
  294   {
  295   /*24:*/
  296   #line 1092 "./wiener.w"
  297   
  298   double**aa;
  299   unsigned long mm= KK;
  300   unsigned dd= 1;
  301   int seed= DEFAULT_SEED;
  302   short mode= MODE_WIENER_PROCESS,verbose= 0;
  303   int no_arg;
  304   register int j,k;
  305   
  306   /*:24*/
  307   #line 636 "./wiener.w"
  308   
  309   /*25:*/
  310   #line 1106 "./wiener.w"
  311   
  312   progname= strip_away_path(argv[0]);
  313   no_arg= argc;
  314   while(--argc){
  315   if(cmd_match("-h","--help",argv[no_arg-argc])){
  316   display_help_message();
  317   exit(SUCCESS);
  318   }else if(cmd_match("-v","--verbose",argv[no_arg-argc])){
  319   verbose= (verbose?0:1);
  320   }else if(cmd_match("-M","--num_samples",argv[no_arg-argc])){
  321   --argc;
  322   if(!sscanf(argv[no_arg-argc],"%lu",&mm)){
  323   fprintf(stderr,"%s: Error detected when parsing the number of "
  324   "samples.\n",progname);
  325   display_help_message();
  326   exit(FAILURE);
  327   }
  328   if(mm<KK){
  329   fprintf(stderr,"%s: The M number of data points must be at least "
  330   "the long lag of the generator, M >= KK = %d.\n",progname,KK);
  331   exit(FAILURE);
  332   }
  333   if(is_odd(mm))mm++;
  334   }else if(cmd_match("-D","--dimension",argv[no_arg-argc])){
  335   --argc;
  336   if(!sscanf(argv[no_arg-argc],"%ud",&dd)){
  337   fprintf(stderr,"%s: Error detected when parsing dimension.\n",
  338   progname);
  339   display_help_message();
  340   exit(FAILURE);
  341   }
  342   if(dd<1){
  343   fprintf(stderr,"%s: Dimension D should be at least 1.\n",progname);
  344   exit(FAILURE);
  345   }
  346   }else if(cmd_match("-s","--seed",argv[no_arg-argc])){
  347   --argc;
  348   if(!sscanf(argv[no_arg-argc],"%d",&seed)){
  349   fprintf(stderr,"%s: Error detected when parsing the seed of the "
  350   "initializer.\n",progname);
  351   display_help_message();
  352   exit(FAILURE);
  353   }
  354   }else if(cmd_match("-u","--uniform",argv[no_arg-argc])){
  355   mode= MODE_LOCKED_UNIFORM_DISTRIBUTION;
  356   }else if(cmd_match("-n","--normal",argv[no_arg-argc])){
  357   mode= MODE_LOCKED_NORMAL_DISTRIBUTION;
  358   }else{
  359   fprintf(stderr,"%s: Sorry, I do not recognize option '%s'.\n",
  360   progname,argv[no_arg-argc]);
  361   display_help_message();
  362   exit(FAILURE);
  363   }
  364   }
  365   if(verbose)
  366   fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
  367   
  368   /*:25*/
  369   #line 637 "./wiener.w"
  370   
  371   /*26:*/
  372   #line 1169 "./wiener.w"
  373   
  374   aa= dmatrix(0,mm-1,0,dd-1);
  375   
  376   /*:26*/
  377   #line 638 "./wiener.w"
  378   
  379   /*27:*/
  380   #line 1174 "./wiener.w"
  381   
  382   wiener(aa,mm,dd,seed,mode);
  383   
  384   /*:27*/
  385   #line 639 "./wiener.w"
  386   
  387   /*28:*/
  388   #line 1179 "./wiener.w"
  389   
  390   for(k= 0;k<mm;k++){
  391   for(j= 0;j<dd-1;j++)
  392   printf("%.20f ",aa[k][j]);
  393   printf("%.20f\n",aa[k][j]);
  394   }
  395   
  396   /*:28*/
  397   #line 640 "./wiener.w"
  398   
  399   /*29:*/
  400   #line 1188 "./wiener.w"
  401   
  402   free_dmatrix(aa,0,mm-1,0,dd-1);
  403   
  404   /*:29*/
  405   #line 641 "./wiener.w"
  406   
  407   return(SUCCESS);
  408   }
  409   
  410   /*:10*/
  411   

Return to previous page

Generated by ::viewsrc::

Last modified Friday 20 Apr 2012