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