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
Generated by ::viewsrc::