Search:

Return to previous page

Contents of file 'wiener/wiener.tex':



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

Return to previous page

Generated by ::viewsrc::

Last modified Wednesday 15 Feb 2023