Search:

Return to previous page

Contents of file 'wiener/wiener.w':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   % File:        wiener.w [CWEB source code]
    2   % Residence:   http://jonsson.eu/programs/cweb/wiener/
    3   % Created:     November 10, 2011 [v.1.0]
    4   % Last change: November 10, 2011 [v.1.0]
    5   % Author:      Fredrik Jonsson <http://jonsson.eu>
    6   % Description: The CWEB source code for the WIENER simulator, generating a
    7   %              sequence of data points corresponding to the Wiener process.
    8   %              The program employs Donald Knuth's proposed random number
    9   %              generator, as listed in and the Box-Muller transform. For
   10   %              information on the CWEB programming language, see
   11   %              http://www.literateprogramming.com.
   12   % Compilation: Compile this program by using the enclosed Makefile, or use
   13   %              the blocks of the Makefile as listed in the documentation
   14   %              (file wiener.ps or wiener.pdf). The C source code (as
   15   %              generated from this CWEB code) conforms to the ANSI standard
   16   %              for the C programming language (which is equivalent to the
   17   %              ISO C90 standard for C).
   18   %
   19   % Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>
   20   % Non-commercial copying welcome.
   21   %
   22   \input epsf
   23   %% \input amstex
   24   \def\version{1.0}
   25   \def\lastrevdate{November 11, 2011}
   26   \font\eightcmr=cmr8
   27   \font\tensc=cmcsc10
   28   \font\eightcmssq=cmssq8
   29   \font\eightcmssqi=cmssqi8
   30   \font\twentycmcsc=cmcsc10 at 20 truept
   31   \def\wiener{{\eightcmr WIENER\spacefactor1000}}
   32   \def\poincare{{\eightcmr POINCARE\spacefactor1000}}
   33   \def\ANSI{{\eightcmr ANSI\spacefactor1000}} % The language standard we stick to
   34   \def\GNU{{\eightcmr GNU\spacefactor1000}}       % GNU is Not Unix
   35   \def\GCC{{\eightcmr GCC\spacefactor1000}}       % The GNU C-compiler
   36   \def\CEE{{\eightcmr C\spacefactor1000}}         % The C programming language
   37   \def\ANSICEE{{\eightcmr ANSI~C\spacefactor1000}}% The ANSI-C standard language
   38   \def\CWEB{{\eightcmr CWEB\spacefactor1000}}     % The CWEB programming language
   39   \def\UNIX{{\eightcmr UNIX\spacefactor1000}}
   40   \def\CYGWIN{{\eightcmr CYGWIN\spacefactor1000}}
   41   \def\GNUPLOT{{\eightcmr GNUPLOT\spacefactor1000}}
   42   \def\CTANGLE{{\eightcmr CTANGLE\spacefactor1000}}
   43   \def\CWEAVE{{\eightcmr CWEAVE\spacefactor1000}}
   44   \def\mod{\mathop{\rm mod}\nolimits}  % The real part of a complex number
   45   \def\re{\mathop{\rm Re}\nolimits}  % The real part of a complex number
   46   \def\im{\mathop{\rm Im}\nolimits}  % The imaginary part of a complex number
   47   %% \def\bbr{\mathbb{R}}             % The blackboard bold 'R' as in 'R^3'
   48   \def\backslash{\char'134}          % The `\' character
   49   \def\vertbar{\char'174}            % The `|' character
   50   \def\dollar{\char'044}             % The `$' character
   51   \def\tilde{\char'176}              % The `~' character
   52   \def\tothepower{\char'136}         % The `^' character
   53   \def\onehalf{{\textstyle{{1}\over{2}}}}
   54   \def\threefourth{{\textstyle{{3}\over{4}}}}
   55   \def\endalg{\vrule height 1.4ex width .6ex depth .4ex} % Rectangular bullet
   56   \def\ie{i.\thinspace{e.}~\ignorespaces}
   57   \def\eg{e.\thinspace{g.}~\ignorespaces}
   58   \let\,\thinspace
   59   %% \def\subsection[#1]{\goodbreak\noindent{\it #1}\par\nobreak\noindent}
   60   %
   61   % Define a handy macro for listing the steps of an algorithm.
   62   %
   63   \newdimen\aitemindent \aitemindent=26pt
   64   \newdimen\aitemleftskip \aitemleftskip=36pt
   65   \def\aitem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   66     \hbox to\aitemindent{\bf #1\hfill}%
   67     \hangindent\aitemleftskip\ignorespaces}
   68   %
   69   % Define a handy macro for bibliographic references.
   70   %
   71   \newdimen\refitemindent \refitemindent=18pt
   72   \def\refitem[#1]{\smallbreak\noindent%
   73     \hbox to\refitemindent{[#1]\hfill}%
   74     \hangindent\refitemindent\ignorespaces}
   75   \def\references{\medskip\goodbreak\noindent{\bf References}\nobreak\smallskip}
   76   %
   77   % Define a handy macro for nicely typeset variable descriptions.
   78   %
   79   \newdimen\varitemindent \varitemindent=100pt
   80   \def\varitem[#1]#2{\smallbreak\noindent\hbox to 20pt{}%
   81     \hbox to\varitemindent{#1\hfill}%
   82     \hangindent 120pt\ignorespaces#2\par}
   83   %
   84   % Define a handy macro for nicely typeset descriptions of command line options.
   85   %
   86   \newdimen\optitemindent \optitemindent=80pt
   87   \def\optitem[#1]#2{\smallbreak\noindent\hbox to\hsize{\hskip40pt{#1}\hfill}
   88     \hangindent\optitemindent\ignorespaces{\hskip\optitemindent{#2}}\smallskip}
   89   %
   90   % Define a handy macro for the list of program revisions.
   91   %
   92   \newdimen\citemindent \citemindent=80pt
   93   \newdimen\citemleftskip \citemleftskip=90pt
   94   \def\citem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   95     \hbox to\citemindent{\bf #1\hfill}%
   96     \hangindent\citemleftskip\ignorespaces}
   97   \def\citindent{\smallbreak\noindent\hbox to 10pt{}%
   98     \hbox to\citemindent{\hfil}%
   99     \hangindent\citemleftskip\ignorespaces}
  100   %
  101   % Define a handy macro for the typesetting of figure captions. The double
  102   % '{{' and '}}' are here in place to prevent the increased \leftskip and
  103   % \rightskip to apply globally after the macro has been expanded.
  104   %
  105   \newdimen\figcapindent \figcapindent=40pt
  106   \def\figcaption[#1]#2{{\advance\leftskip by\figcapindent
  107     \advance\rightskip by\figcapindent
  108     \smallskip\noindent{\bf Figure #1.} {#2}\smallskip}}
  109   %
  110   % Define the \beginvrulealign and \endvrulealign macros as described in
  111   % Donald Knuth's The TeXbook, Appendix D: Dirty Tricks. These macros are
  112   % used in typesetting nicely looking tables.
  113   %
  114   \def\beginvrulealign{\setbox0=\vbox\bgroup}
  115   \def\endvrulealign{\egroup % now \box0 holds the entire alignment
  116     \setbox0=\vbox{\setbox2=\hbox{\vrule height\ht0 depth\dp0 width0pt}
  117       \unvbox0 \setbox0=\lastbox % now \box0 is the bottom row
  118       \nointerlineskip \copy0 % put it back
  119       \global\setbox1=\hbox{} % initialize box that will contain rules
  120       \setbox4=\hbox{\unhbox0 % now open up the bottom row
  121         \loop \skip0=\lastskip \unskip % remove tabskip glue
  122         \advance\skip0 by-.4pt % rules are .4 pt wide
  123         \divide\skip0 by 2
  124         \global\setbox1=\hbox{\hskip\skip0\vrule\hskip\skip0
  125           \unhbox2\unhbox1}%
  126         \setbox2=\lastbox % remove alignment entry
  127         \ifhbox2 \setbox2=\hbox{\kern\wd2}\repeat}}%
  128     \hbox{\rlap{\box0}\box1}} % superimpose the alignment on the rules
  129   %
  130   % Make sure that the METAPOST logo can be loaded even in plain TeX.
  131   %
  132   \ifx\MP\undefined
  133       \ifx\manfnt\undefined
  134               \font\manfnt=logo10
  135       \fi
  136       \ifx\manfntsl\undefined
  137               \font\manfntsl=logosl10
  138       \fi
  139       \def\MP{{\ifdim\fontdimen1\font>0pt \let\manfnt = \manfntsl \fi
  140         {\manfnt META}\-{\manfnt POST}}\spacefactor1000 }%
  141   \fi
  142   \ifx\METAPOST\undefined \let\METAPOST=\MP \fi
  143   
  144   \datethis
  145   
  146   @*Introduction.
  147   \vskip 60pt
  148   \centerline{\twentycmcsc Wiener}
  149   \vskip 20pt
  150   \centerline{A computer program for the generation of numerical data
  151     simulating a Wiener process.}
  152   \vskip 2pt
  153   \centerline{(Version \version\ of \lastrevdate)}
  154   \vskip 10pt
  155   \centerline{Written by Fredrik Jonsson}
  156   \vskip 10pt
  157   \centerline{\epsfxsize=88mm\epsfbox{fig3.1}}
  158   \vskip 10pt
  159   \noindent
  160   This \CWEB\footnote{${}^\dagger$}{For information on the \CWEB\ programming
  161   language by Donald E.~Knuth, as well as samples of \CWEB\ programs, see
  162   {\tt http://www-cs-faculty.stanford.edu/\~\ \kern -5pt knuth/cweb.html}.
  163   For general information on literate programming, see
  164   {\tt http://www.literateprogramming.com}.}
  165   computer program computes a series of floating-point numbers
  166   corresponding to a Wiener process in $D$ dimensions. It relies on the random
  167   number generator as proposed by Donald Knuth in {\sl The Art of Computer
  168   Programming, Volume 1 -- Fundamental Algorithms}, 3rd edition
  169   (Addison-Wesley, Boston, 1998), generating numbers which are fed into the
  170   Box--Muller transform to generate the normal distribution associated with the
  171   Wiener process.
  172   Besides providing a simulator of the Wiener process, the \wiener\ program
  173   can also be used in a ``lock-in'' mode with zero expectation value for each
  174   data point, providing a pretty good random number generator for large series
  175   of stochastic data, not relying on the (rather poor) generators available in
  176   standard \CEE\ libraries.
  177   
  178   The \wiener\ program does not solve any problem {\sl per se}, but is merely
  179   to be considered as a generator of statistical data to be used by other
  180   applications in modeling of physical, chemical or financial processes.
  181   \bigskip
  182   \noindent Copyright \copyright\ Fredrik Jonsson, 2011.
  183   All rights reserved. Non-commercial copying welcome.
  184   \vfill
  185   
  186   @*The Wiener process.
  187   ``In mathematics, the Wiener process is a continuous-time stochastic process
  188   named in honor of Norbert Wiener. It is often called standard Brownian motion,
  189   after Robert Brown. It is one of the best known L\'evy processes (c\`adl\`ag
  190   stochastic processes with stationary independent increments) and occurs
  191   frequently in pure and applied mathematics, economics and physics.
  192   
  193   The Wiener process plays an important role both in pure and applied
  194   mathematics. In pure mathematics, the Wiener process gave rise to the study
  195   of continuous time martingales. It is a key process in terms of which more
  196   complicated stochastic processes can be described. As such, it plays a vital
  197   role in stochastic calculus, diffusion processes and even potential theory.
  198   It is the driving process of Schramm--Loewner evolution.
  199   In applied mathematics, the Wiener process is used to represent the integral
  200   of a Gaussian white noise process, and so is useful as a model of noise in
  201   electronics engineering, instrument errors in filtering theory and unknown
  202   forces in control theory.
  203   
  204   The Wiener process has applications throughout the mathematical sciences. In
  205   physics it is used to study Brownian motion, the diffusion of minute particles
  206   suspended in fluid, and other types of diffusion via the Fokker--Planck and
  207   Langevin equations. It also forms the basis for the rigorous path integral
  208   formulation of quantum mechanics (by the Feynman--Kac formula, a solution to
  209   the Schr\"odinger equation can be represented in terms of the Wiener process)
  210   and the study of eternal inflation in physical cosmology. It is also prominent
  211   in the mathematical theory of finance, in particular the Black--Scholes option
  212   pricing model.''
  213   \smallskip
  214   {\eightcmssq --\,Wikipedia, ``Wiener process'' (2011)}
  215   
  216   @ What the \wiener\ program does (and doesn't).
  217   The present \CWEB\ program does not solve any problems related to any of the
  218   processes described by models involving the Wiener process, but is merely an
  219   attempt to produce an as-good-as-possible result when simulating the Wiener
  220   process as such.
  221   In the \wiener\ program, special attention has been paid to the generation
  222   of random numbers, as this is a crucial and rather tricky problem when it
  223   comes to generating large non-recurring series of data.
  224   In the present program, the random number generator proposed by Donald
  225   Knuth\footnote{$\dagger$}{Donald E.~Knuth, {\sl The Art of Computer
  226   Programming, Volume 1 -- Fundamental Algorithms}, 3rd edition (Addison-Wesley,
  227   Boston, 1998).} has been employed, generating uniformly distributed numbers
  228   which are fed into the Box--Muller transform to generate the normal
  229   distribution associated with the Wiener process.
  230   
  231   Apart from being a pretty good and reliable generator of statistical data, to
  232   be used by other applications in modeling of physical, chemical or financial
  233   processes, the \wiener\ program does not solve any problems {\sl per se}.
  234   
  235   @ The CWEB programming language.
  236   For the reader who might not be familiar with the concept of the
  237   \CWEB\ programming language, the following citations hopefully will
  238   be useful. For further information, as well as freeware compilers for
  239   compiling \CWEB\ source code, see {\tt http://www.literateprogramming.com}.
  240   \bigskip
  241   
  242   {\narrower\narrower\narrower\narrower\eightcmssqi\noindent
  243   I believe that the time is ripe for significantly better documentation of
  244   programs, and that we can best achieve this by considering programs to be
  245   works of literature. Hence, my title: `Literate Programming.'
  246   
  247   Let us change our traditional attitude to the construction of programs:
  248   Instead of imagining that our main task is to instruct a computer what to
  249   do, let us concentrate rather on explaining to human beings what we want
  250   a computer to do.
  251   
  252   The practitioner of literate programming can be regarded as an essayist,
  253   whose main concern is with exposition and excellence of style. Such an
  254   author, with thesaurus in hand, chooses the names of variables carefully
  255   and explains what each variable means. He or she strives for a program
  256   that is comprehensible because its concepts have been introduced in an
  257   order that is best for human understanding, using a mixture of formal and
  258   informal methods that reinforce each other.
  259   \smallskip
  260   {\eightcmssq --\,Donald Knuth,}
  261   {\eightcmssqi The CWEB System of Structured Documentation}
  262   {\eightcmssq (Addison-Wesley, Massachusetts, 1994)}
  263   }
  264   \bigskip
  265   
  266   {\narrower\narrower\narrower\narrower\eightcmssqi\noindent
  267   The philosophy behind CWEB is that an experienced system programmer, who
  268   wants to provide the best possible documentation of his or her software
  269   products, needs two things simultaneously: a language like \TeX\ for
  270   formatting, and a language like C for programming. Neither type of language
  271   can provide the best documentation by itself; but when both are appropriately
  272   combined, we obtain a system that is much more useful than either language
  273   separately.
  274   
  275   The structure of a software program may be thought of as a `WEB' that is
  276   made up of many interconnected pieces. To document such a program we want to
  277   explain each individual part of the web and how it relates to its neighbors.
  278   The typographic tools provided by \TeX\ give us an opportunity to explain the
  279   local structure of each part by making that structure visible, and the
  280   programming tools provided by languages like C make it possible for us to
  281   specify the algorithms formally and unambiguously. By combining the two,
  282   we can develop a style of programming that maximizes our ability to perceive
  283   the structure of a complex piece of software, and at the same time the
  284   documented programs can be mechanically translated into a working software
  285   system that matches the documentation.
  286   
  287   Besides providing a documentation tool, CWEB enhances the C language by
  288   providing the ability to permute pieces of the program text, so that a
  289   large system can be understood entirely in terms of small sections and
  290   their local interrelationships. The CTANGLE program is so named because
  291   it takes a given web and moves the sections from their web structure into
  292   the order required by C; the advantage of programming in CWEB is that the
  293   algorithms can be expressed in ``untangled'' form, with each section
  294   explained separately. The CWEAVE program is so named because it takes a
  295   given web and intertwines the \TeX\ and C portions contained in each
  296   section, then it knits the whole fabric into a structured document.
  297   \smallskip
  298   {\eightcmssq --\,Donald Knuth, ``Literate Programming'', in}
  299   {\eightcmssqi Literate Programming}
  300   {\eightcmssq (CSLI Lecture Notes, Stanford, 1992)}
  301   }
  302   
  303   @ Revision history of the program.
  304   \medskip
  305   
  306   \citem[2011-11-11]{[v.1.0]} {\tt <http://jonsson.eu/programs/cweb/>}\hfill\break
  307   First properly working version of the \wiener\ program, written in \CWEB\
  308   and (\ANSI-conformant) \CEE.
  309   
  310   @*Compiling the source code. The program is written in \CWEB, generating
  311   \ANSICEE\ (ISO C90) conforming source code and documentation as plain
  312   \TeX-source, and is to be compiled using the sequences as outlined in the
  313   {\tt Makefile} listed below.
  314   \bigskip
  315   {\obeylines\obeyspaces\tt
  316   \#
  317   \# Makefile designed for use with ctangle, cweave, gcc, and plain TeX.
  318   \#
  319   \# Copyright (C) 2002-2011, Fredrik Jonsson <http://jonsson.eu>
  320   \#
  321   \# The CTANGLE program converts a CWEB source document into a C program
  322   \# which may be compiled in the usual way. The output file includes \#line
  323   \# specifications so that debugging can be done in terms of the CWEB source
  324   \# file.
  325   \#
  326   \# The CWEAVE program converts the same CWEB file into a TeX file that may
  327   \# be formatted and printed in the usual way. It takes appropriate care of
  328   \# typographic details like page layout and the use of indentation, italics,
  329   \# boldface, etc., and it supplies extensive cross-index information that it
  330   \# gathers automatically.
  331   \#
  332   \# CWEB allows you to prepare a single document containing all the informa-
  333   \# tion that is needed both to produce a compilable C program and to produce
  334   \# a well-formatted document describing the program in as much detail as the
  335   \# writer may desire.  The user of CWEB ought to be familiar with TeX as well
  336   \# as C.
  337   \#
  338   PROJECT  = wiener
  339   FIGURES  = fig1.eps fig2.eps fig3.eps
  340   ~ ~
  341   CTANGLE  = ctangle
  342   CWEAVE   = cweave
  343   CC       = gcc
  344   CCOPTS   = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
  345   LNOPTS   = -lm
  346   TEX      = tex
  347   DVIPS    = dvips
  348   DVIPSOPT = -ta4 -D1200
  349   PS2PDF   = ps2pdf
  350   METAPOST = mpost
  351   ~ ~
  352   all: \$(PROJECT) \$(FIGURES) \$(PROJECT).pdf
  353   ~ ~
  354   \$(PROJECT): \$(PROJECT).o
  355   ~        \$(CC) \$(CCOPTS) -o \$(PROJECT) \$(PROJECT).o \$(LNOPTS)
  356   ~ ~
  357   \$(PROJECT).o: \$(PROJECT).c
  358   ~        \$(CC) \$(CCOPTS) -c \$(PROJECT).c
  359   ~ ~
  360   \$(PROJECT).c: \$(PROJECT).w
  361   ~        \$(CTANGLE) \$(PROJECT)
  362   ~ ~
  363   \$(PROJECT).pdf: \$(PROJECT).ps
  364   ~        \$(PS2PDF) \$(PROJECT).ps \$(PROJECT).pdf
  365   ~ ~
  366   \$(PROJECT).ps: \$(PROJECT).dvi
  367   ~        \$(DVIPS) \$(DVIPSOPT) \$(PROJECT).dvi -o \$(PROJECT).ps
  368   ~ ~
  369   \$(PROJECT).dvi: \$(PROJECT).tex
  370   ~        \$(TEX) \$(PROJECT).tex
  371   ~ ~
  372   \$(PROJECT).tex: \$(PROJECT).w
  373   ~        \$(CWEAVE) \$(PROJECT)
  374   ~ ~
  375   \#
  376   \# Generate the Encapsulated PostScript image fig1.eps for the documentation.
  377   \# This is a 2D scatter plot of the uniformly distributed pseudo-random numbers
  378   \# prior to having been fed into the Box-Muller transform.
  379   \#
  380   fig1.eps: Makefile \$(PROJECT).w
  381   ~        wiener --uniform -D 2 -M 10000 > fig1.dat;
  382   ~        @echo "input graph;{\backslash}
  383   ~        ~        def mpdot = btex{\backslash}
  384   ~        ~        ~  {{\backslash}vrule height 0.5pt width 1.0pt depth 0.5pt}{\backslash}
  385   ~        ~        etex enddef;{\backslash}
  386   ~        ~        beginfig(1);{\backslash}
  387   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  388   ~        ~        setrange(0,0,1,1);{\backslash}
  389   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  390   ~        ~        gdraw {\backslash}"fig1.dat{\backslash}" plot mpdot;{\backslash}
  391   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  392   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  393   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  394   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  395   ~        ~        endgraph; endfig; end" > fig1.mp
  396   ~        ~        \$(METAPOST) fig1.mp
  397   ~        \$(TEX) -jobname=fig1 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  398   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig1.1}}{\backslash}bye"
  399   ~        \$(DVIPS) -D1200 -E fig1.dvi -o fig1.eps
  400   ~ ~
  401   \#
  402   \# Generate the Encapsulated PostScript image fig2.eps for the documentation.
  403   \# This is a 2D scatter plot of the normally distributed pseudo-random numbers
  404   \# resulting from the Box-Muller transform.
  405   \#
  406   fig2.eps: Makefile \$(PROJECT).w
  407   ~        wiener --normal -D 2 -M 10000 > fig2.dat;
  408   ~        @echo "input graph;{\backslash}
  409   ~        ~        def mpdot = btex{\backslash}
  410   ~        ~        ~  {{\backslash}vrule height 0.5pt width 1.0pt depth 0.5pt}{\backslash}
  411   ~        ~        etex enddef;{\backslash}
  412   ~        ~        beginfig(1);{\backslash}
  413   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  414   ~        ~        setrange(whatever,whatever,whatever,whatever);{\backslash}
  415   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  416   ~        ~        gdraw {\backslash}"fig2.dat{\backslash}" plot mpdot;{\backslash}
  417   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  418   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  419   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  420   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  421   ~        ~        endgraph; endfig; end" > fig2.mp
  422   ~        ~        \$(METAPOST) fig2.mp
  423   ~        \$(TEX) -jobname=fig2 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  424   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig2.1}}{\backslash}bye"
  425   ~        \$(DVIPS) -D1200 -E fig2.dvi -o fig2.eps
  426   ~ ~
  427   \#
  428   \# Generate the Encapsulated PostScript image fig3.eps for the documentation.
  429   \# This is a 2D graph showing the resulting simulated Wiener process.
  430   \#
  431   fig3.eps: Makefile \$(PROJECT).w
  432   ~        wiener -D 2 -M 10000 > fig3.dat;
  433   ~        @echo "input graph;{\backslash}
  434   ~        ~        beginfig(1);{\backslash}
  435   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  436   ~        ~        setrange(whatever,whatever,whatever,whatever);{\backslash}
  437   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  438   ~        ~        gdraw {\backslash}"fig3.dat{\backslash}";{\backslash}
  439   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  440   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  441   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  442   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  443   ~        ~        endgraph; endfig; end" > fig3.mp
  444   ~        ~        \$(METAPOST) fig3.mp
  445   ~        \$(TEX) -jobname=fig3 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  446   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig3.1}}{\backslash}bye"
  447   ~        \$(DVIPS) -D1200 -E fig3.dvi -o fig3.eps
  448   ~ ~
  449   clean:
  450   ~        -rm -Rf \$(PROJECT) *~ *.c *.o *.exe *.dat *.pdf *.mp *.trj *.mpx
  451   ~        -rm -Rf *.tex *.aux *.log *.toc *.idx *.scn *.dvi *.ps *.1 *.eps
  452   ~ ~
  453   archive:
  454   ~        make -ik clean
  455   ~        tar --gzip --directory=../ -cf ../\$(PROJECT).tgz \$(PROJECT)
  456   }
  457   \bigskip
  458   \noindent
  459   This {\tt Makefile} essentially executes two major calls. First, the \CTANGLE\
  460   program parses the \CWEB\ source document {\tt wiener.w} to extract a \CEE\
  461   source file {\tt wiener.c} which may be compiled in the usual way using any
  462   \ANSICEE\ conformant compiler. The output source file includes {\tt \#line}
  463   specifications so that any debugging can be done conveniently in terms of
  464   the original \CWEB\ source file.
  465   Second, the \CWEAVE\ program parses the same \CWEB\ source file to extract a
  466   plain \TeX\ source file {\tt wiener.tex} which may be compiled in the usual
  467   way.
  468   It takes appropriate care of typographic details like page layout and the use
  469   of indentation, italics, boldface, and so on, and it supplies extensive
  470   cross-index information that it gathers automatically.
  471   
  472   After having executed {\tt make} in the same catalogue where the files
  473   {\tt wiener.w} and {\tt Makefile} are located, one is left with an
  474   executable file {\tt wiener}, being the ready-to-use compiled program,
  475   and a PostScript file {\tt wiener.ps}
  476   which contains the full documentation of the program, that is to say the
  477   document you currently are reading.
  478   On platforms running any operating system by Microsoft, the executable file
  479   will instead automatically be called {\tt wiener.exe}.
  480   This convention also applies to programs compiled under the \UNIX-like
  481   environment \CYGWIN.
  482   
  483   @*Running the program. The program is entirely controlled by the command
  484   line options supplied when invoking the program. The syntax for executing the
  485   program is\par
  486   \medskip
  487   \line{\hskip 40pt{\tt wiener [options]}\hfill}\par
  488   \medskip
  489   \noindent
  490   where {\tt options} include the following, given in their long (`{\tt --}') as
  491   well as their short (`{\tt -}') forms:
  492   \medskip
  493   \optitem[{\tt --help}, {\tt -h}]{Display a brief help message and exit clean.}
  494   \optitem[{\tt --verbose}, {\tt -v}]{Toggle verbose mode. Default: off.}
  495   \optitem[{\tt --num\_samples}, {\tt -M} $M$]{Generate $M$ samples of data. Here
  496     $M$ should always be an even number, greater than the long lag |KK|. If an
  497     odd number is specified, the program will automatically adjust this to the
  498     next higher integer. Default: $M=|KK|=100$.}
  499   \optitem[{\tt --dimension}, {\tt -D} $D$]{Specifies the dimension $D$ of the
  500     Wiener process, that is to say generating a set of $D$ numbers for each of
  501     the $M$ data points in the seqence. Default: $D=1$.}
  502   \optitem[{\tt --seed}, {\tt -s} $S$]{Define a custom seed number $S$ for the
  503     initialization of the random number generator.
  504     Default: $S=|DEFAULT_SEED|=310952$.}
  505   \optitem[{\tt --uniform}, {\tt -u}]{Instead of generating a sequence of data
  506     corresponding to a Wiener process, lock the program to simply generate a
  507     uniform distribution of $D$-dimensional points, with each element distributed
  508     over the interval $[0,1]$.}
  509   \optitem[{\tt --normal}, {\tt -n}]{Instead of generating a sequence of data
  510     corresponding to a Wiener process, lock the program to simply generate a
  511     normal distribution of $D$-dimensional points, with each element distributed
  512     with zero expectation value and unit variance.}
  513   \noindent
  514   One may look upon the two last options as verification options, generating data
  515   suitable for spectral tests on the quality of the generator of pseudo-random
  516   numbers.
  517   
  518   @ Plotting the results using \GNUPLOT.
  519   The data sets generated by \wiener\ may be displayed
  520   and saved as Encapsulated PostScript images using, say, \GNUPLOT\ or \METAPOST.
  521   While I personally prefer MetaPost, mainly due to the possibility of
  522   incorporating the same typygraphic elements as in \TeX, many people consider
  523   \GNUPLOT\ to be easier in operation.
  524   
  525   In order to save a scatter graph as Encapsulated PostScript using \GNUPLOT, you
  526   may include the following calls in, say, a {\tt Makefile} or a shell script:
  527   \bigskip
  528   {\obeylines\obeyspaces\tt
  529   ~      wiener -D 2 -M 10000 > figure.dat;
  530   ~      echo "set term postscript eps;\backslash
  531   ~             set output \\"figure.eps\\";\backslash
  532   ~             set size square;\backslash
  533   ~             set nolabel;\backslash
  534   ~             plot \\"figure.dat\\" with dots notitle;\backslash
  535   ~             quit" \vertbar\ gnuplot
  536   }
  537   
  538   @ Plotting the results using \METAPOST.
  539   Another choice is to go for the \METAPOST\ way. This is illustrated with the
  540   following blocks, taken directly from the enclosed Makefile and generating
  541   the figures which can be seen in the section relating to the generation of
  542   normally distributed variables (routine |normdist|):
  543   \bigskip
  544   {\obeylines\obeyspaces\tt
  545   PROJECT  = wiener
  546   TEX      = tex
  547   DVIPS    = dvips
  548   METAPOST = mpost
  549   ~ ~
  550   \#
  551   \# Generate the Encapsulated PostScript image fig1.eps for the documentation.
  552   \# This is a 2D scatter plot of the uniformly distributed pseudo-random numbers
  553   \# prior to having been fed into the Box-Muller transform.
  554   \#
  555   fig1.eps: Makefile \$(PROJECT).w
  556   ~        wiener --uniform -D 2 -M 10000 > fig1.dat;
  557   ~        @echo "input graph;{\backslash}
  558   ~        ~        def mpdot = btex{\backslash}
  559   ~        ~        ~  {{\backslash}vrule height 0.5pt width 1.0pt depth 0.5pt}{\backslash}
  560   ~        ~        etex enddef;{\backslash}
  561   ~        ~        beginfig(1);{\backslash}
  562   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  563   ~        ~        setrange(0,0,1,1);{\backslash}
  564   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  565   ~        ~        gdraw {\backslash}"fig1.dat{\backslash}" plot mpdot;{\backslash}
  566   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  567   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  568   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  569   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  570   ~        ~        endgraph; endfig; end" > fig1.mp
  571   ~        ~        \$(METAPOST) fig1.mp
  572   ~        \$(TEX) -jobname=fig1 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  573   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig1.1}}{\backslash}bye"
  574   ~        \$(DVIPS) -D1200 -E fig1.dvi -o fig1.eps
  575   ~ ~
  576   \#
  577   \# Generate the Encapsulated PostScript image fig2.eps for the documentation.
  578   \# This is a 2D scatter plot of the normally distributed pseudo-random numbers
  579   \# resulting from the Box-Muller transform.
  580   \#
  581   fig2.eps: Makefile \$(PROJECT).w
  582   ~        wiener --normal -D 2 -M 10000 > fig2.dat;
  583   ~        @echo "input graph;{\backslash}
  584   ~        ~        def mpdot = btex{\backslash}
  585   ~        ~        ~  {{\backslash}vrule height 0.5pt width 1.0pt depth 0.5pt}{\backslash}
  586   ~        ~        etex enddef;{\backslash}
  587   ~        ~        beginfig(1);{\backslash}
  588   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  589   ~        ~        setrange(whatever,whatever,whatever,whatever);{\backslash}
  590   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  591   ~        ~        gdraw {\backslash}"fig2.dat{\backslash}" plot mpdot;{\backslash}
  592   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  593   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  594   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  595   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  596   ~        ~        endgraph; endfig; end" > fig2.mp
  597   ~        ~        \$(METAPOST) fig2.mp
  598   ~        \$(TEX) -jobname=fig2 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  599   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig2.1}}{\backslash}bye"
  600   ~        \$(DVIPS) -D1200 -E fig2.dvi -o fig2.eps
  601   ~ ~
  602   \#
  603   \# Generate the Encapsulated PostScript image fig3.eps for the documentation.
  604   \# This is a 2D graph showing the resulting simulated Wiener process.
  605   \#
  606   fig3.eps: Makefile \$(PROJECT).w
  607   ~        wiener -D 2 -M 10000 > fig3.dat;
  608   ~        @echo "input graph;{\backslash}
  609   ~        ~        beginfig(1);{\backslash}
  610   ~        ~        draw begingraph(86mm,86mm);{\backslash}
  611   ~        ~        setrange(whatever,whatever,whatever,whatever);{\backslash}
  612   ~        ~        pickup pencircle scaled .5pt;{\backslash}
  613   ~        ~        gdraw {\backslash}"fig3.dat{\backslash}";{\backslash}
  614   ~        ~        pickup pencircle scaled .25pt;{\backslash}
  615   ~        ~        autogrid(itick bot,itick lft);{\backslash}
  616   ~        ~        glabel.bot(btex \$\$ x\$\$ etex,OUT);{\backslash}
  617   ~        ~        glabel.lft(btex \$\$ y\$\$ etex,OUT);{\backslash}
  618   ~        ~        endgraph; endfig; end" > fig3.mp
  619   ~        ~        \$(METAPOST) fig3.mp
  620   ~        \$(TEX) -jobname=fig3 "{\backslash}input epsf{\backslash}nopagenumbers{\backslash}
  621   ~        ~   {\backslash}centerline{{\backslash}epsfbox{fig3.1}}{\backslash}bye"
  622   ~        \$(DVIPS) -D1200 -E fig3.dvi -o fig3.eps
  623   ~ ~
  624   }
  625   
  626   @*The main program. Here follows the general outline of the main program.
  627   
  628   @c
  629   @<Library inclusions@>@;
  630   @<Global definitions@>@;
  631   @<Global variables@>@;
  632   @<Routines@>@;
  633   
  634   int main(int argc, char *argv[])
  635   {
  636      @<Declaration of local variables@>@;
  637      @<Parse command line for parameters@>@;
  638      @<Allocate memory for a vector containing $M\times D$ elements@>@;
  639      @<Fill vector with $M$ number of $D$:tuples describing the Wiener process@>@;
  640      @<Print the generated vector at standard terminal output@>@;
  641      @<Deallocate the memory occupied by the vector of $M\times D$ elements@>@;
  642      return(SUCCESS);
  643   }
  644   
  645   @ Library dependencies.
  646   The standard \ANSICEE\ libraries included in this program are:\par
  647   \varitem[{\tt math.h}]{For access to common mathematical functions.}
  648   \varitem[{\tt stdio.h}]{For file access and any block involving |fprintf|.}
  649   \varitem[{\tt stdlib.h}]{For access to the |exit| function.}
  650   \varitem[{\tt string.h}]{For string manipulation, |strcpy|, |strcmp| etc.}
  651   \varitem[{\tt ctype.h}]{For access to the |isalnum| function.}
  652   
  653   @<Library inclusions@>=
  654   #include <math.h>
  655   #include <stdio.h>
  656   #include <stdlib.h>
  657   #include <string.h>
  658   #include <ctype.h>
  659   
  660   @ Global definitions.
  661   These are the global definitions present in the \wiener\ program:\par
  662   \varitem[{\tt VERSION}]{The current program revision number.}
  663   \varitem[{\tt COPYRIGHT}]{The copyright banner.}
  664   \varitem[{\tt SUCCESS}]{The return code for successful program termination.}
  665   \varitem[{\tt FAILURE}]{The return code for unsuccessful program termination.}
  666   \varitem[|QUALITY|]{The recommended ``quality level'' for high-resolution
  667            use, according to Knuth. Used by the |ranf_array| routine.}
  668   \varitem[{\tt KK}]{The ``long lag'' used by routines |ranf_array| and
  669                      |ranf_start|.}
  670   \varitem[{\tt LL}]{The ``short lag'' used by routines |ranf_array| and
  671                      |ranf_start|.}
  672   \varitem[|DEFAULT_SEED|]{The default seed to use when initializing the random
  673            number generator, using the routine |ranf_start|. The seed can be
  674            hand-picked using the {\tt --seed} command-line option.}
  675   \varitem[{\tt cmd\_match(s,l,c)}]{Check if the string |s| and/or |l| matches
  676            the string |c|. This short-hand macro is used when parsing the
  677            command line for options.}
  678   
  679   @<Global definitions@>=
  680   #define VERSION "1.0"
  681   #define COPYRIGHT "Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>"
  682   #define SUCCESS (0)
  683   #define FAILURE (1)
  684   #define QUALITY (1009)
  685   #define KK (100)
  686   #define LL  (37)
  687   #define DEFAULT_SEED (310952)
  688   #define cmd_match(s,l,c) ((!strcmp((c),(s)))||(!strcmp((c),(l))))
  689   #define MODE_WIENER_PROCESS (0)
  690   #define MODE_LOCKED_UNIFORM_DISTRIBUTION (1)
  691   #define MODE_LOCKED_NORMAL_DISTRIBUTION (2)
  692   
  693   @ Declaration of global variables. Usually, the only global variables allowed
  694   in my programs are |optarg|, which is a pointer to the the string of characters
  695   that specified the call from the command line, and |progname|, which simply
  696   is a pointer to the string containing the name of the program, as it was
  697   invoked from the command line. However, as Donald Knuth has a {\sl faiblesse}
  698   for global variables, I have for the sake of consistency with the routines
  699   related to the random number generator kept his definitions in a global scope.
  700   
  701   @<Global variables@>=
  702   extern char *optarg;
  703   char *progname;
  704   double ranf_arr_buf[QUALITY];
  705   double ranf_arr_dummy=-1.0, ranf_arr_started=-1.0;
  706   double *ranf_arr_ptr=&ranf_arr_dummy; /* the next random fraction, or -1 */
  707   
  708   @*Declaration of routines. These routines exclusively concern the generation
  709   of random numbers and the generation of normally distributed data points in
  710   the Wiener process.
  711   
  712   @<Routines@>=
  713     @<Routine for generation of uniformly distributed random numbers@>@;
  714     @<Routine for initialization of the random number generator@>@;
  715     @<Routine for generation of normally distributed variables@>@;
  716     @<Routine for generation of numerical data describing a Wiener process@>@;
  717     @<Routine for memory allocation@>@;
  718     @<Routine for memory de-allocation@>@;
  719     @<Routine for displaying a help message at terminal output@>@;
  720     @<Routine for checking valid path characters@>@;
  721     @<Routine for stripping away path string@>@;
  722   
  723   @ Generation of uniformly distributed random numbers. We here
  724   avoid relying on the standard functions available in \CEE,\footnote{$\dagger$}{
  725   Here only included as a reference, the (primitive) standard \CEE\ way
  726   of generating random integer numbers, in this case initializing with a simple
  727   {\tt srand(time(NULL))} and generating random numbers between 0 and |RAND_MAX|,
  728   is
  729   {\obeylines\obeyspaces\tt
  730      int rand\_stdc() {
  731         srand(time(NULL)); /* Initialize random number generator */
  732         return(rand());
  733      }
  734   }
  735   \noindent
  736   I have personally {\sl not} (yet) checked this approach using the spectral
  737   tests recommended by Knuth.} but rather take resort to the algorithm devised by
  738   Donald Knuth in {\sl The Art of Computer Programming, Volume 1 -- Fundamental
  739   Algorithms}, 3rd edition, Section 3.6. (Addison--Wesley, Boston, 1998).%
  740   \footnote{$\ddagger$}{
  741   The credit for this random number generator, which employs double
  742   floating-point precision rather than the original integer version,
  743   goes entirely to Donald Knuth and
  744   the persons who contributed. The original source code %and full list of credits
  745   are available at
  746   {\tt http://www-cs-faculty.stanford.edu/{\tilde}knuth/programs/rng-double.c}.
  747   The current routine takes into account changes as explained in the errata to
  748   the 2nd edition, see
  749   {\tt http://www-cs-faculty.stanford.edu/{\tilde}knuth/taocp.html} in the
  750   changes to Volume 2 on pages 171 and following.}
  751   The variables to the routine are as follows:
  752   \varitem[{\tt double aa[]}]{On return, this array contains |n| new random
  753   numbers, following the recurrence $X_j=(X_{j-100}-X_{j-37})\mod2^{30}.$
  754   Not used on input.}
  755   \varitem[{\tt double n}]{The number |n| of random numbers to be generated.
  756   This array length must be at least |KK| elements.}
  757   \noindent
  758   The |mod_sum(x,y)| macro is here merely a shorthand for ``$(x+y)\mod1.0$.''
  759   
  760   @<Routine for generation of uniformly distributed random numbers@>=
  761   #define mod_sum(x,y) (((x)+(y))-(int)((x)+(y)))
  762   double ran_u[KK];           /* the generator state */
  763   
  764   void ranf_array(double aa[], int n)  /* put n new random fractions in aa */
  765   {
  766     register int i,j;
  767   
  768     for (j=0;j<KK;j++) aa[j]=ran_u[j];
  769     for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
  770     for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
  771     for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK],ran_u[i-LL]);
  772   }
  773   
  774   void ranf_matrix(double **aa, int mm, int dd)
  775   {
  776      register int i,j,col;
  777   
  778      for (col=0;col<dd;col++) {
  779         for (j=0;j<KK;j++) aa[j][col]=ran_u[j];
  780         for (;j<mm;j++) aa[j][col]=mod_sum(aa[j-KK][col],aa[j-LL][col]);
  781         for (i=0;i<LL;i++,j++) ran_u[i]=mod_sum(aa[j-KK][col],aa[j-LL][col]);
  782         for (;i<KK;i++,j++) ran_u[i]=mod_sum(aa[j-KK][col],ran_u[i-LL]);
  783      }
  784   }
  785   
  786   @ Initialization of the random number generator. To quote Knuth,
  787   ``The tricky thing about using a recurrence like [$X_j=(X_{j-100}-X_{j-37})
  788   \mod 2^{30}$] is, of course, to get everuthing started properly in the
  789   first place, by setting up suitable values of $X_0,\ldots,X_{99}$.
  790   The following routine |ran_start| initializes the generator nicely when given
  791   any seed number between $0$ and $2^{30}-3=1,073,741,821$ inclusive.''
  792   Here we rather employ the double precision variant of the initialization to
  793   match the data type of |ran_array|.
  794   
  795   Again, the credits for the |ranf_start| routine goes entirely to Donald Knuth
  796   and the persons who contributed.
  797   
  798   @<Routine for initialization of the random number generator@>=
  799   #define TT  70   /* guaranteed separation between streams */
  800   #define is_odd(s) ((s)&1)
  801   
  802   void ranf_start(long seed) {
  803      register int t,s,j;
  804      double u[KK+KK-1];
  805      double ulp=(1.0/(1L<<30))/(1L<<22);               /* 2 to the -52 */
  806      double ss=2.0*ulp*((seed&0x3fffffff)+2);
  807   
  808      for (j=0;j<KK;j++) {
  809         u[j]=ss;                                /* bootstrap the buffer */
  810         ss+=ss;
  811         if (ss>=1.0) ss-=1.0-2*ulp;  /* cyclic shift of 51 bits */
  812      }
  813      u[1]+=ulp;                     /* make u[1] (and only u[1]) "odd" */
  814      for (s=seed&0x3fffffff,t=TT-1; t; ) {
  815         for (j=KK-1;j>0;j--)
  816            u[j+j]=u[j],u[j+j-1]=0.0;                         /* "square" */
  817         for (j=KK+KK-2;j>=KK;j--) {
  818            u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]);
  819            u[j-KK]=mod_sum(u[j-KK],u[j]);
  820         }
  821         if (is_odd(s)) {                             /* "multiply by z" */
  822            for (j=KK;j>0;j--) u[j]=u[j-1];
  823            u[0]=u[KK];                    /* shift the buffer cyclically */
  824            u[LL]=mod_sum(u[LL],u[KK]);
  825         }
  826         if (s) s>>=1; else t--;
  827      }
  828      for (j=0;j<LL;j++) ran_u[j+KK-LL]=u[j];
  829      for (;j<KK;j++) ran_u[j-LL]=u[j];
  830      for (j=0;j<10;j++) ranf_array(u,KK+KK-1);  /* warm things up */
  831      ranf_arr_ptr=&ranf_arr_started;
  832   }
  833   
  834   #define ranf_arr_next() (*ranf_arr_ptr>=0? *ranf_arr_ptr++: ranf_arr_cycle())
  835   double ranf_arr_cycle()
  836   {
  837     if (ranf_arr_ptr==&ranf_arr_dummy)
  838       ranf_start(314159L); /* the user forgot to initialize */
  839     ranf_array(ranf_arr_buf,QUALITY);
  840     ranf_arr_buf[KK]=-1;
  841     ranf_arr_ptr=ranf_arr_buf+1;
  842     return ranf_arr_buf[0];
  843   }
  844   
  845   @ Generation of normally distributed variables. Accepting uniformly
  846   distributed random numbers as input, the Box--Muller transform creates a set of
  847   normally distributed numbers. This transform originally appeared in
  848   G.~E. P.~Box and Mervin E.~Muller, {\sl A Note on the Generation of Random
  849   Normal Deviates}, Annals Math.~Stat. {\bf 29}, 610--611 (1958).
  850   
  851   In Donald Knuth's {\sl The Art of Computer Programming, Volume 1 -- Fundamental
  852   Algorithms}, 3rd edition (Addison--Wesley, Boston, 1998), the
  853   Box--Muller method is denoted {the polar method} and is described in detail
  854   in Section 3.4.1, Algorithm P, on page 122.
  855   
  856   \centerline{\epsfbox{fig1.1}}
  857   \figcaption[1]{Sample two-dimensional output from the |ranf_matrix| routine,
  858     in this case 10\,000 data points uniformly distributed over the domain
  859     $0\le\{x,y\}\le 1$. The data for this graph was generated by \wiener\
  860     using {\tt wiener --uniform -D 2 -M 10000 > fig1.dat}.
  861     See the {\tt fig1.eps} block in the enclosed {\tt Makefile} for details on
  862     how \METAPOST\ was used in the generation of the encapsulated PostScript
  863     image of the graph.}
  864   
  865   \centerline{\epsfbox{fig2.1}}
  866   \figcaption[2]{The same data points as in Fig.~1, but after having applied
  867     the Box--Muller transform to yield a normal distribution of pseudo-random
  868     numbers. The data for this graph was generated by \wiener\
  869     using {\tt wiener --normal -D 2 -M 10000 > fig2.dat}.
  870     See the {\tt fig2.eps} block in the enclosed {\tt Makefile} for details on
  871     how \METAPOST\ was used in the generation of the encapsulated PostScript
  872     image of the graph.}
  873   
  874   @<Routine for generation of normally distributed variables@>=
  875   #define PI (3.14159265358979323846264338)
  876   void normdist(double **aa, int mm, int dd) {
  877      register int j, k;
  878      register double f, z;
  879   
  880      for (j=0;j<dd;j++) {
  881         for (k=0;k<mm-1;k+=2) {
  882            if (aa[k][j]>0.0) {
  883               f=sqrt(-2*log(aa[k][j]));
  884               z=2.0*PI*aa[k+1][j];
  885               aa[k][j]=f*cos(z);
  886               aa[k+1][j]=f*sin(z);
  887            } else {
  888               fprintf(stderr,"%s: Error: Zero element detected!\n",progname);
  889               fprintf(stderr,"%s: (row %d, column %d)\n",progname,k,j);
  890            }
  891         }
  892      }
  893      return;
  894   }
  895   #undef PI
  896   
  897   @ Routine for generation of numerical data describing a Wiener process. This
  898   is the core routine of the \wiener\ program.
  899   After having initialized the random number generator with the supplied seed
  900   value (calling |ranf_start(seed)|), the generation of a sequence of
  901   numbers describing a Wiener process starts with the generation of $M\times D$
  902   random and uniformly distributed floating-point numbers ($M\times D/2$ pairs,
  903   assuming $M\times D$ to be an even number), from which $M\times D$ normally
  904   distributed floating-point numbers are computed using the Box--Muller
  905   transform\footnote{$\dagger$}{For example, see
  906   {\tt http://en.wikipedia.org/wiki/Box\%E2\%80\%93Muller\_transform}.}
  907   
  908   The actual computation of the random numbers (at first corresponding to a
  909   uniform distribution in the interval $[0,1]$) is done by the routine
  910   |ranf_matrix(aa,mm,dd)|, which fills an $[M\times D]$ array.
  911   
  912   \centerline{\epsfbox{fig3.1}}
  913   \figcaption[3]{The same data points as in Fig.~2, but after having chained
  914     the normally distributed points to form the simulated Wiener process.
  915     The data for this graph was generated by \wiener\
  916     using {\tt wiener -D 2 -M 10000 > fig3.dat}.
  917     The trajectory starts with data point 1 at $(0,0)$ and end up with data
  918     point 10\,000 at approximately $(89.9,12.6)$
  919     See the {\tt fig3.eps} block in the enclosed {\tt Makefile} for details on
  920     how \METAPOST\ was used in the generation of the encapsulated PostScript
  921     image of the graph.}
  922   
  923   The variables interfacing the |wiener| routine are as follows:
  924   \medskip
  925   \varitem[|aa|]{[Output] The $M\times D$ matrix $A$, on return containing the
  926     $M$ number of $D$-dimensional data points for the generated Wiener process.}
  927   \varitem[|mm|]{[Input] The $M$ number of data points to generate. This equals
  928     to the number of rows in the |aa| array, and must be at least |KK| elements.}
  929   \varitem[|dd|]{[Input] The dimension of each generated data point.}
  930   \varitem[|seed|]{[Input] The seed to use when initializing the random number
  931     generator, using the routine |ranf_start|.}
  932   \varitem[|mode|]{[Input] Determines if the sequence should be
  933     locked to simply generate a uniform distribution over the interval $[0,1]$
  934     (|MODE_LOCKED_UNIFORM_DISTRIBUTION|) or a normal (Gaussian) distribution with
  935     expectation value zero and unit variance (|MODE_LOCKED_NORMAL_DISTRIBUTION|).
  936     Otherwise, the series of data will be generated to simulate a Wiener process,
  937     as is the main purpose of the \wiener\ program.
  938     One may look upon the two first modes as verification options, generating
  939     data suitable for spectral tests on the quality of the generator of
  940     pseudo-random numbers.}
  941   
  942   @<Routine for generation of numerical data describing a Wiener process@>=
  943   void wiener(double **aa, int mm, int dd, int seed, short mode)
  944   {
  945      register int j, k;
  946   
  947      ranf_start(seed);
  948      ranf_matrix(aa,mm,dd); /* Uniform distribution over $[0,1]$ */
  949      if (mode==MODE_LOCKED_UNIFORM_DISTRIBUTION) return;
  950      normdist(aa, mm, dd); /* Normal distribution of unit variance around zero */
  951      if (mode==MODE_LOCKED_NORMAL_DISTRIBUTION) return;
  952      for (j=0;j<dd;j++) {
  953         aa[0][j]=0.0;
  954         for (k=1;k<mm;k++) {
  955            aa[k][j]+=aa[k-1][j];
  956         }
  957      }
  958   }
  959   
  960   @ Memory allocation.
  961   The |dmatrix| routine allocates an array of double floating-point precision,
  962   with array row index ranging from |nrl| to |nrh| and column index ranging
  963   from |ncl| to |nch|.
  964   
  965   @<Routine for memory allocation@>=
  966   double **dmatrix(long nrl, long nrh, long ncl, long nch)
  967   {
  968      long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  969      double **m;
  970      m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
  971      if (!m) {
  972         fprintf(stderr,"%s: Allocation failure 1 in dmatrix()\n",progname);
  973         exit(FAILURE);
  974      }
  975      m += 1;
  976      m -= nrl;
  977      m[nrl]=(double*) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
  978      if (!m[nrl]) {
  979         fprintf(stderr,"%s: Allocation failure 2 in dmatrix()\n",progname);
  980         exit(FAILURE);
  981      }
  982      m[nrl] += 1;
  983      m[nrl] -= ncl;
  984      for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  985      return m;
  986   }
  987   
  988   @ Memory de-allocation.
  989   The |free_dmatrix| routine {\sl releases} the memory occupied by the
  990   double floating-point precision matrix |v[nl..nh]|, as allocated by |dmatrix()|.
  991   
  992   @<Routine for memory de-allocation@>=
  993   void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch) {
  994      free((char*) (m[nrl]+ncl-1));
  995      free((char*) (m+nrl-1));
  996   }
  997   
  998   @ Displaying a help message at terminal output.
  999   
 1000   @<Routine for displaying a help message at terminal output@>=
 1001   void display_help_message(void) {
 1002      fprintf(stderr,"Usage: %s M [options]\n",progname);
 1003      fprintf(stderr,"Options:\n");
 1004      fprintf(stderr," -h, --help\n");
 1005      fprintf(stderr,"    Display this help message and exit clean.\n");
 1006      fprintf(stderr," -v, --verbose\n");
 1007      fprintf(stderr,"    Toggle verbose mode. Default: off.\n");
 1008      fprintf(stderr," -M, --num_samples <M>\n");
 1009      fprintf(stderr,"    Generate M samples of data. Here M should always be\n");
 1010      fprintf(stderr,"    an even number, greater than the long lag KK=%d.\n",KK);
 1011      fprintf(stderr,"    If an odd number is specified, the program will\n");
 1012      fprintf(stderr,"    automatically adjust this to the next higher\n");
 1013      fprintf(stderr,"    integer. Default: M=%d.\n",KK);
 1014      fprintf(stderr," -D, --dimension <D>\n");
 1015      fprintf(stderr,"    Generate D-dimensional samples of data. Default: D=1.\n");
 1016      fprintf(stderr," -s, --seed <seed>\n");
 1017      fprintf(stderr,"    Define a custom seed number for the initialization\n");
 1018      fprintf(stderr,"    of the random number generator. Default: seed=%d.\n",
 1019         DEFAULT_SEED);
 1020      fprintf(stderr," -u, --uniform\n");
 1021      fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
 1022      fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
 1023      fprintf(stderr,"    to simply generate a uniform distribution of\n");
 1024      fprintf(stderr,"    D-dimensional points, with each element distributed\n");
 1025      fprintf(stderr,"    over the interval [0,1].\n");
 1026      fprintf(stderr," -n, --normal\n");
 1027      fprintf(stderr,"    Instead of generating a sequence of D-dimensional\n");
 1028      fprintf(stderr,"    corresponding to a Wiener process, lock the program\n");
 1029      fprintf(stderr,"    to simply generate a normal distribution of\n");
 1030      fprintf(stderr,"    D-dimensional points, with each element distributed\n");
 1031      fprintf(stderr,"    with zero expectation value and unit variance.\n");
 1032   }
 1033   
 1034   @ Checking for a valid path character.
 1035   The |pathcharacter| routine takes one character |ch| as argument, and returns
 1036   1 (``true'') if the character is valid character of a path string, otherwise 0
 1037   (``false'') is returned.
 1038   
 1039   @<Routine for checking valid path characters@>=
 1040   short pathcharacter(int ch) {
 1041      return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
 1042          (ch=='-')||(ch=='+'));
 1043   }
 1044   
 1045   @ Stripping path string from a file name.
 1046   The |strip_away_path| routine takes a character string |filename| as
 1047   argument, and returns a pointer to the same string but without any
 1048   preceding path segments. This routine is, for example, useful for
 1049   removing paths from program names as parsed from the command line.
 1050   
 1051   @<Routine for stripping away path string@>=
 1052   char *strip_away_path(char filename[]) {
 1053      int j,k=0;
 1054      while (pathcharacter(filename[k])) k++;
 1055      j=(--k); /* this is the uppermost index of the full path+file string */
 1056      while (isalnum((int)(filename[j]))) j--;
 1057      j++; /* this is the lowermost index of the stripped file name */
 1058      return(&filename[j]);
 1059   }
 1060   
 1061   @ Declaration of local variables of the |main| program.
 1062   In \CWEB\ one has the option of adding variables along the program, for example
 1063   by locally adding temporary variables related to a given sub-block of code.
 1064   However, the philosophy in the \wiener\ program is to keep all variables of
 1065   the |main| section collected together, so as to simplify tasks as, for example,
 1066   tracking down a given variable type definition. Exceptions to this rule are
 1067   dummy variables merely used for the iteration over loops, not participating in
 1068   any other tasks.
 1069   The local variables of the program are as follows:
 1070   \medskip
 1071   \varitem[|aa|]{The $M\times D$ matrix $A$, containing the $M$ number of
 1072     $D$-dimensional data points for the generated Wiener process.}
 1073   \varitem[|mm|]{The $M$ number of data points to generate. This equals to the
 1074     number of rows in the |aa| array (of dimension $[M\times D]$), and must be
 1075     at least |KK| elements. The default initialization is $|mm|=|KK|$; however
 1076     this may change depending on supplied command-line parameters.}
 1077   \varitem[|dd|]{The dimension of each generated data point. Default value: 1.}
 1078   \varitem[|seed|]{The seed to use when initializing the random number generator,
 1079            using the routine |ranf_start|. The seed can be hand-picked using the
 1080            {\tt --seed} command-line option. Default value: 310952.}
 1081   \varitem[|mode|]{Determines if the sequence should be
 1082     locked to simply generate a uniform distribution over the interval $[0,1]$
 1083     (|MODE_LOCKED_UNIFORM_DISTRIBUTION|) or a normal (Gaussian)
 1084     distribution with expectation value zero and unit variance
 1085     (|MODE_LOCKED_NORMAL_DISTRIBUTION|).
 1086     Otherwise, the series of data will be generated to simulate a Wiener process,
 1087     as is the main purpose of the \wiener\ program.
 1088     One may look upon the two first modes as verification options, generating
 1089     data suitable for spectral tests on the quality of the generator of
 1090     pseudo-random numbers.}
 1091   
 1092   @<Declaration of local variables@>=
 1093      double **aa;
 1094      unsigned long mm=KK;
 1095      unsigned dd=1;
 1096      int seed=DEFAULT_SEED;
 1097      short mode=MODE_WIENER_PROCESS, verbose=0;
 1098      int no_arg;
 1099      register int j, k;
 1100   
 1101   @ Parse command line for parameters.
 1102   We here use the possibility open in \CWEB\ to add {\tt getopt.h} to the
 1103   inclusions of libraries, as this block is the only one making use of the
 1104   definitions therein.
 1105   
 1106   @<Parse command line for parameters@>=
 1107      progname=strip_away_path(argv[0]);
 1108      no_arg=argc;
 1109      while(--argc) {
 1110         if(cmd_match("-h", "--help", argv[no_arg-argc])) {
 1111            display_help_message();
 1112            exit(SUCCESS);
 1113         } else if(cmd_match("-v", "--verbose", argv[no_arg-argc])) {
 1114            verbose=(verbose?0:1);
 1115         } else if(cmd_match("-M", "--num_samples", argv[no_arg-argc])) {
 1116            --argc;
 1117            if(!sscanf(argv[no_arg-argc],"%lu",&mm)) {
 1118               fprintf(stderr,"%s: Error detected when parsing the number of "
 1119                  "samples.\n", progname);
 1120               display_help_message();
 1121               exit(FAILURE);
 1122            }
 1123            if (mm<KK) {
 1124               fprintf(stderr,"%s: The M number of data points must be at least "
 1125                  "the long lag of the generator, M >= KK = %d.\n", progname, KK);
 1126               exit(FAILURE);
 1127            }
 1128            if (is_odd(mm)) mm++; /* If odd, then make it even */
 1129         } else if(cmd_match("-D", "--dimension", argv[no_arg-argc])) {
 1130            --argc;
 1131            if(!sscanf(argv[no_arg-argc],"%ud",&dd)) {
 1132               fprintf(stderr,"%s: Error detected when parsing dimension.\n",
 1133                  progname);
 1134               display_help_message();
 1135               exit(FAILURE);
 1136            }
 1137            if (dd<1) {
 1138               fprintf(stderr,"%s: Dimension D should be at least 1.\n",progname);
 1139               exit(FAILURE);
 1140            }
 1141         } else if(cmd_match("-s", "--seed", argv[no_arg-argc])) {
 1142            --argc;
 1143            if(!sscanf(argv[no_arg-argc],"%d",&seed)) {
 1144               fprintf(stderr,"%s: Error detected when parsing the seed of the "
 1145                  "initializer.\n", progname);
 1146               display_help_message();
 1147               exit(FAILURE);
 1148            }
 1149         } else if(cmd_match("-u", "--uniform", argv[no_arg-argc])) {
 1150            mode=MODE_LOCKED_UNIFORM_DISTRIBUTION;
 1151         } else if(cmd_match("-n", "--normal", argv[no_arg-argc])) {
 1152            mode=MODE_LOCKED_NORMAL_DISTRIBUTION;
 1153         } else {
 1154            fprintf(stderr,"%s: Sorry, I do not recognize option '%s'.\n",
 1155               progname, argv[no_arg-argc]);
 1156            display_help_message();
 1157            exit(FAILURE);
 1158         }
 1159      }
 1160      if (verbose)
 1161         fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
 1162   
 1163   @ Allocate memory for a vector containing $M\times D$ elements. We here make a
 1164   call to the operating system for the allocation of a sufficcient amount of
 1165   memory to accommodate $M$ data points, each of dimension $D$.
 1166   We here apply the common convention in \CEE, starting the indexing of the
 1167   allocated array at zero.
 1168   
 1169   @<Allocate memory for a vector containing $M\times D$ elements@>=
 1170      aa=dmatrix(0,mm-1,0,dd-1);
 1171   
 1172   @ Fill vector with $M$ number of $D$:tuples describing the Wiener process.
 1173   
 1174   @<Fill vector with $M$ number of $D$:tuples describing the Wiener process@>=
 1175      wiener(aa,mm,dd,seed,mode);
 1176   
 1177   @ Print the generated vector at standard terminal output.
 1178   
 1179   @<Print the generated vector at standard terminal output@>=
 1180      for (k=0;k<mm;k++) {
 1181         for (j=0;j<dd-1;j++)
 1182            printf("%.20f ", aa[k][j]);
 1183         printf("%.20f\n", aa[k][j]);
 1184      }
 1185   
 1186   @ Deallocate the memory occupied by the vector of $M\times D$ elements.
 1187   
 1188   @<Deallocate the memory occupied by the vector of $M\times D$ elements@>=
 1189      free_dmatrix(aa,0,mm-1,0,dd-1);
 1190   
 1191   @*Index.
 1192   

Return to previous page

Generated by ::viewsrc::

Last modified Friday 20 Apr 2012