- Gpc - gnu.de

Tree SSA
by Adriaan van Os 27 May '04

27 May '04

Tree SSA has been merged into gcc mainline, see <http://gcc.gnu.org/projects/tree-ssa/>. This seems to be a major change, e.g. the Ada and Fortran frontends are broken. Any comments from Frank or Waldek or from anybody else on the list ? Regards, Adriaan van Os

3 2

BP (in)compatibility
by Waldek Hebisch 27 May '04

27 May '04

I noticed that the following FPC code (claimed to be BP compatible) is not accepted by GPC: function b:string; begin b[1] := ' ' end; Similarely FPC accepts: type a = record i, j : integer end; function b(i, j: integer) : a; begin b.i := i; b.j := j end; I do not know if we want to implement that feature, but (assuming it is really accepted by BP) we should list it as known incompatiblity. -- Waldek Hebisch hebisch(a)math.uni.wroc.pl

5 9

BP OOP problem
by Frank Heckenbach 25 May '04

25 May '04

I just found another problem in BP's OOP. In the following program, `o.f' in oo.p is ambiguous. It could mean field f of the field o, or function f of the parent type o. Both BP and GPC (currently) do the former (i.e., say 99). I noticed this because the OOE draft forbids this case by requiring that object fields don't clash with parent type names. We do we do here? Allow it in BP mode with a warning, error otherwise? program Foo; type o = object function f: Integer; end; oo = object (o) o: record f: Integer end; procedure p; end; function o.f: Integer; begin f := 42 end; procedure oo.p; begin WriteLn (o.f) end; var v: oo; begin v.o.f := 99; v.p end. Frank -- Frank Heckenbach, frank(a)g-n-u.de, http://fjf.gnu.de/, 7977168E GPC To-Do list, latest features, fixed bugs: http://www.gnu-pascal.de/todo.html GPC download signing key: 51FF C1F0 1A77 C6C2 4482 4DDC 117A 9773 7F88 1707

4 17

BP question
by Frank Heckenbach 20 May '04

20 May '04

BP seems to allow overriding a public method with a private one (and vice versa). Does this seem reasonable? What do Delphi and other compilers do? At first sight, I'd rather forbid this in GPC (as usual, with an exception in `--borland-pascal', but perhaps a warning even there?). program Foo; type a = object procedure p; end; b = object (a) private procedure p; { ??? } end; procedure a.p; begin end; procedure b.p; begin end; begin end. Oh, wait a second, it gets really strange ... unit U; interface type a = object procedure p; end; b = object (a) private procedure p; end; implementation procedure a.p; begin WriteLn ('a.p') end; procedure b.p; begin WriteLn ('b.p') end; var v: b; begin v.p end. program Foo; uses U; var v: b; begin v.p end. Under BP, this says: b.p a.p What does it mean actually? It means that `b' has *two* methods called `p'. Within the unit, one of them shadows the other, and outside of the unit, the first one is hidden (or non-existent because not exported), and the second one becomes visible. That's crazy, isn't it? Frank -- Frank Heckenbach, frank(a)g-n-u.de, http://fjf.gnu.de/, 7977168E GPC To-Do list, latest features, fixed bugs: http://www.gnu-pascal.de/todo.html GPC download signing key: 51FF C1F0 1A77 C6C2 4482 4DDC 117A 9773 7F88 1707

5 11

Re: IEEE 754 Standard (was: gpc-20030830 released)
by Tom Verhoeff 20 May '04

20 May '04

Adriaan van Os wrote: > >Frank Heckenbach wrote: > >> ISO Pascal wants an error (though, in principle, I think we could >> use NaN and document the error as "not detected" -- we might do this >> sometime, perhaps optionally). For now, I'm making it runtime >> errors. (emil27*.pas) > >Special cases for the pow function in IEEE 754 (see e.g. Apple >Numerics Manual, second edition, page 64, or the PowerPC Numerics >manual available from ><http://developer.apple.com/documentation/mac/PPCNumerics/PPCNumerics- >2.html>). > >Operation Result Exceptions raised [ List suppressed ] It would be a good thing if gpc provided full support of the IEEE 754 Standard. Unfortunatly (for numerical applications, at least), many compilers either provide no support at all, or provide only partial support. The latter is often even worse than the former. (Apple has indeed been a loyal supporter of the standard from the very beginning.) The need for some features in the IEEE 754 standard is less obvious to the uninitiated. Good articles in this respect are Kahan, W. Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic. May 1996. <http://www.cs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF> and Kahan, W. Why do we need a floating-point arithmetic standard? Unpublished note, February 1981. <http://www.cs.berkeley.edu/~dbindel/class/cs279/why-ieee.pdf> William Kahan is the intellectual father of the standard. Another important references is (the title says it all): Goldberg, D. ``What every computer scientist should know about floating-point arithmetic'', Appendix~D in Numerical Computation Guide. Sun Microsystems, Revision A, May 2000; originally published in ACM Computing Surveys (March 1991). <http://docs.sun.com/source/806-3568/ncg_goldberg.html> The following (short) book is a good introduction to the use of IEEE 754: Overton, M. L. Numerical Computing with IEEE Floating Point Arithmetic. SIAM, 2001. For those with less time, I recommend an article I recently wrote: Gyula Horvath, Tom Verhoeff. ``Numerical Difficulties in Pre-University Informatics Education and Competitions'', Informatics in Education, Vol. 2, Number 1, pp.21-38. <http://www.win.tue.nl/~wstomv/publications/INFE012.pdf> Here is one simple example contained in that article, to explain the need for avoiding unnecessary exception raising. Suppose you want to write a piece of program to calculate the effective resistance of two (and later more than two) parallel resistors, R_1 and R_2, both non-negative real numbers (Ohms). The basic formula is 1 R_Eff = --------- (i.e. the "harmonic" mean of R_1 and R_2) 1 1 --- + --- R_1 R_2 Of course, it fails (under normal math rules), when R_1 or R_2 (or both) are zero. However, when either (or both) are zero, the effective resistance is well-defined, viz. zero. A slightly better formula may seem to be R_1 * R_2 R_Eff = --------- R_1 + R_2 because it also "works" when either of R_1 and R_2 are zero. However, there are two problems with that formula: (1) when both are zero, it gives 0/0, and (2) it does not generalize easily to more than two resistors. Now, under the IEEE 754 Standard, the first formula works perfectly well. Why? Because the standard includes signed infinities and NaN (not-a-number) that are well-behaved: 1 / +0 = +inf (zero is also signed!) a + inf = inf 1 / +inf = +0 0 / 0 = NaN 0 * +inf = NaN NaN <any-operator> a = NaN With these rules, the first formula works correctly under all circumstances, without requiring any case analysis in the program. However, if division by zero would raise an exception, then you need to catch it in the program, and you get a clumsy program. The second formula fails under these rules for R_1 = R_2 = 0, because it involves 0 / 0. In (numerical) practice, there are even "worse" examples (such as rational function with multiple "false" poles; see Kahan's articles), where a case analysis would be horrible compared to the straightforward approach via IEEE 754. All in all, the committee that developed the IEEE 754 Standard thought hard about the features, and they came up with a carefully balanced set. Partial implementation is not recommended, because it destroys some important properties. Tom

4 3

Someone's getting old
by Neil Santos 20 May '04

20 May '04

Someone forgot to change the homepage on gnu-pascal.de about a month ago. Not that I'm complaining; the message is an important one. But maybe the stuff with dates should be taken out? But what do I know? ;b -- Political Dissident Neil Santos Free Software Advocate Jabber: nsantos(a)jabber.org Please avoid sending me attachments encoded in secret formatting. See http://www.fsf.org/philosophy/no-word-attachments.html

2 1

Re: gp-0.54 released
by Frank Heckenbach 18 May '04

18 May '04

Thorsten Glaser wrote: > Dixitur illum frank(a)g-n-u.de scribere... > > >one of the next releases, and automake support will later be dropped > >from GPC. > > Would it be too difficult to keep --automake support, Yes! > If one would add a feature that prints out the names > of the generated (or not, depending on if they have > been rebuilt) .o files, this feature integrated nicely > with BSD-style Makefile templates, and I was already > planning to code that. Such a feature could also be added to GP. In fact, I've thought about it, any it may be faster than doing it with GPC/automake (because GP only takes a quick look at source files, whereas GPC/automake parses them fully). My use for it would be building the RTS and GP itself (which obviously can't use GP in the normal way), though it's not too important for me, since the list of files here changes very rarely, so maintaining an explicit makefile will be possible. But note that even if you use Makefiles, you can still use GP (unless you absolutely want to avoid installing it :-). Just use a rule such as: foo: foo.pas force gp <options> foo.pas Where `foo.pas' is the main source file (GP will find all other source files -- if the necessary paths are given, of course), and `force' makes sure that GP is called each time, so it can check dependencies. (Consider it on the same level as a recursive sub-make, e.g.) If nothing needs to be done, GP will not compile or link anything, so it isn't costly. If your program depends on generated files, you'll have to list them in the dependencies of the rule, but nothing more is required AFAICS. > Of course, you need not retain automake support if > you could make better compiler code otherwise, but > as long as it's no real effort, I'd appreciate the > feature being left in. It is much effort, and as I wrote, automake is still buggy and incomplete (e.g., it doesn't work well with cyclic dependencies). Frank -- Frank Heckenbach, frank(a)g-n-u.de, http://fjf.gnu.de/, 7977168E GPC To-Do list, latest features, fixed bugs: http://www.gnu-pascal.de/todo.html GPC download signing key: 51FF C1F0 1A77 C6C2 4482 4DDC 117A 9773 7F88 1707

2 1

Re: gp-0.54 released
by Frank Heckenbach 18 May '04

18 May '04

Thorsten Glaser wrote: > Dixitur illum ih8mj(a)fjf.gnu.de scribere... > > >Such a feature could also be added to GP. In fact, I've thought > >about it, any it may be faster than doing it with GPC/automake > >(because GP only takes a quick look at source files, whereas > >GPC/automake parses them fully). > > Hmm, but the pascal compiler is still gpc... To extract dependencies you don't need a compiler. > >But note that even if you use Makefiles, you can still use GP > > Our Makefiles look like > > --- snip > # $Id$ > > PROG= foo > SRCS= foo.pas > > .include <bsd.prog.mk> > --- snip > > which triggers ${PC} ${PFLAGS} -c foo.pas, and linking. I don't know exactly what it does. My point was that the Makefile should invoke GP unconditionally. Of course, it could explicitly list all dependencies, but it's easier (and won't need amending when dependencies change) to call GP and let it find out when nothing needs to be done. Frank -- Frank Heckenbach, frank(a)g-n-u.de, http://fjf.gnu.de/, 7977168E GPC To-Do list, latest features, fixed bugs: http://www.gnu-pascal.de/todo.html GPC download signing key: 51FF C1F0 1A77 C6C2 4482 4DDC 117A 9773 7F88 1707

1 0

Results for 3.3.4 20040516 (prerelease) testsuite on i686-pc-linux-gnu
by David Kredba 18 May '04

18 May '04

Native configuration is i686-pc-linux-gnu === gpc tests === Running target any === gpc Summary === # of tests 3910 # of expected passes 3907 # of unsupported tests 3 /home/dave/gcc/gcc-3.3.4/gcc/xgpc version 20040516, based on gcc-3.3.4 20040516 (prerelease) === gcc tests === Running target unix FAIL: gcc.dg/duff-2.c (test for excess errors) === gcc Summary === # of expected passes 21868 # of unexpected failures 1 # of expected failures 69 # of unsupported tests 100 /home/dave/gcc/gcc-3.3.4/gcc/xgcc version 3.3.4 20040516 (prerelease) === g++ tests === Running target unix FAIL: g++.other/decl5.C caused compiler crash XPASS: g++.other/init5.C Execution test === g++ Summary === # of expected passes 8242 # of unexpected failures 1 # of unexpected successes 1 # of expected failures 91 # of untested testcases 9 # of unsupported tests 15 /home/dave/gcc/gcc-3.3.4/gcc/testsuite/../g++ version 3.3.4 20040516 (prerelease) === g77 tests === Running target unix === g77 Summary === # of expected passes 1736 /home/dave/gcc/gcc-3.3.4/gcc/testsuite/../g77 version 3.3.4 20040516 (prerelease) === objc tests === Running target unix === objc Summary === # of expected passes 1166 /home/dave/gcc/gcc-3.3.4/gcc/xgcc version 3.3.4 20040516 (prerelease) === libjava tests === Running target unix WARNING: program timed out. FAIL: SyncTest execution - gij test WARNING: program timed out. FAIL: SyncTest execution - gij test FAIL: Thread_Interrupt output - source compiled test FAIL: Thread_Interrupt output - gij test FAIL: Thread_Interrupt -O output - source compiled test FAIL: Thread_Interrupt output - gij test === libjava Summary === # of expected passes 2983 # of unexpected failures 6 # of expected failures 16 # of untested testcases 14 === libstdc++-v3 check-abi Summary === # of added symbols: 3 # of missing symbols: 0 # of incompatible symbols: 0 === libstdc++-v3 tests === Running target unix FAIL: 22_locale/codecvt_members_char_char.cc execution test FAIL: 22_locale/codecvt_members_wchar_t_char.cc execution test FAIL: 22_locale/ctor_copy_dtor.cc execution test FAIL: 22_locale/ctype_to_char.cc execution test FAIL: 22_locale/ctype_to_wchar_t.cc execution test FAIL: 22_locale/facet.cc execution test FAIL: 27_io/filebuf_virtuals.cc execution test FAIL: 27_io/ostream_inserter_arith.cc execution test FAIL: 27_io/streambuf_members.cc execution test FAIL: 27_io/stringbuf_virtuals.cc execution test === libstdc++-v3 Summary === # of expected passes 443 # of unexpected failures 10 # of expected failures 28 Compiler version: 3.3.4 20040516 (prerelease) Platform: i686-pc-linux-gnu configure flags: --host=i686-pc-linux-gnu --prefix=/opt --enable-shared --enable-threads=posix --enable-checking --enable-long-long --enable-__cxa_atexit --enable-clocale=gnu --enable-languages=c,c++,ada,f77,objc,java,pascal --with-system-zlib -- David Kredba <kredba at ibot.cas.cz> GPG: ID 1024D/5B6B02DE Fingerprint: F0B3312596BEDCF91DFB 0699E06AACD75B6B02DE

1 0

random number unit
by Toby Ewing 14 May '04

14 May '04

Hi, all I've been working with the Mersenne Twister random number generator for a while. I translated it from the original C to Pascal several years ago, and the GPC team incorporated it into the GPC unit. But since then the original authors have found and fixed a bug in their initilization procedures, and other users have suggested various speedups. I've now updated my original translation, validated it against the authors' test output, performed some other minor tests, and done what small optimizations I know how to do. Source for the unit and a test program are appended below. I don't know much about optimization, and there may be some compiler directives that would greatly speed the unit. I'd appreciate any comments on that front: currently I only use: {$no-pack-struct, maximum-field-alignment 0} and more could probably be done. regards, Toby ------------------------------------------ Unit MTrand; { Mersenne Twister random number generator -- a Pascal unit. Based on code by Makoto Matsumoto, Takuji Nishimura, Shawn Cokus, Richard J. Wagner, and others. This code is based on the Jan 26 2002 version, found on Matsumoto's MT website: http://konatu.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html See also the original website: http://www.math.keio.ac.jp/~matumoto/emt.html Translation to a GPC unit by Toby Ewing (ewing(a)iastate.edu) March 2004, with assistance from James Houston (jhouston(a)tsargent.com) The Mersenne Twister is an algorithm for generating random numbers. It was designed with consideration of the flaws in various other generators. The period, 2^19937-1, and the order of equidistribution, 623 dimensions, are far greater. The generator is also fast; it avoids multiplication and division, and it benefits from caches and pipelines. Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. The names of its contributors may not be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. The original code included the following notice: When you use this, send an email to: matumoto(a)math.keio.ac.jp with an appropriate reference to your work. } {$no-pack-struct, maximum-field-alignment 0} Interface const mtN = 624; {length of state vector} type MTcard = cardinal attribute (Size = 32); PMTary = ^Buffary; Buffary(m : medCard) = array[0..m] of medCard; {Generate uniform random 32-bit cardinal} Function randInt : MTcard; {cardinal in [0,2^32-1]} Function randrInt(range : MTcard) : MTcard; {cardinal in [0,n-1]} {Generate uniform random double} Function rand : double; {double in [0,1]} Function randr(range : double) : double; {double in [0,r]} Function randEx1 : double; {double in [0,1)} Function randrEx1(range : double) : double; {double in [0,r)} Function randEx0 : double; {double in (0,1]} Function randrEx0(range : double) : double; {double in (0,r]} Function randDblExc : double; {double in (0,1)} Function randrDblExc(range : double) : double; {double in (0,r)} {Generate 53-bit uniform random double (capacity of IEEE double precision)} Function rand53 : double; {double in [0,1)} {Generate random double from normal distribution, given mean and stdev} Function randNorm(Mean, stdev : double) : double;{double from (mean, stdev)} {Initialization} Procedure SeedRand(RandSeed : MTcard); {initialize with an MTcard} Procedure SeedRandArray( keylen : medCard; {initialize with an array} var initkey : PMTary); Implementation const mtM = 397; {period parameter} vectA = 16#9908b0df; {constant vector a} UMask = 16#80000000; {most significant bits} LMask = 16#7fffffff; {least significant bits} TempB = 16#9d2c5680; {Tempering parameters B and C} TempC = 16#efc60000; type MTary = array[0..mtN] of MTcard; magary = array[0..1] of MTcard; var mt : MTary; {array for MT state vector} mti : MTcard; {array index} Procedure Reload; {This is the basic vector generator, making mtN random MTcards at a time} var y, kk : MTcard; mag01 : magary = (0, vectA); attribute (static); begin for kk := 0 to mtN-mtM-1 do begin y := (mt[kk] AND UMask) OR (mt[kk+1] AND LMask); mt[kk] := mt[kk+mtM] XOR (y shr 1) XOR mag01[y AND 16#00000001]; end; for kk := mtN-mtM to mtN-1 do begin y := (mt[kk] AND UMask) OR (mt[kk+1] AND LMask); mt[kk] := mt[kk+(mtM-mtN)] XOR (y shr 1) XOR mag01[y AND 16#00000001]; end; y := (mt[mtN-1] AND UMask) OR (mt[0] AND LMask); mt[mtN-1] := mt[mtM-1] XOR (y shr 1) XOR mag01[y AND 16#00000001]; mti := 0; end; Function RandInt : MTcard; {card in [0,2^32-1]} {This is the root function for returning a single random number: all of the other functions draw on this one.} var twist : MTcard; begin if (mti >= mtN) then reload; twist := mt[mti]; inc(mti); twist := twist XOR (twist shr 11); twist := twist XOR ((twist shl 7) AND TempB); twist := twist XOR ((twist shl 15) AND TempC); randInt := twist XOR (twist shr 18); end; {Real versions (in C) by Isaku Wada} Function rand : double; {double in [0,1]} begin rand := randInt * (1.0/4294967295.0); end; Function randr(range : double) : double; {double in [0,r]} begin randr := range * randInt * (1.0/4294967295.0); end; Function randEx1 : double; {double in [0,1)} begin randEx1 := randInt * (1.0/4294967296.0); end; Function randrEx1(range : double) : double; {double in [0,r)} begin randrEx1 := range * randInt * (1.0/4294967296.0); end; Function randEx0 : double; {double in (0,1]} begin randEx0 := (randInt + 1.0) * (1.0/4294967295.0); end; Function randrEx0(range : double) : double; {double in (0,r]} begin randrEx0 := range * (randInt + 1.0) * (1.0/4294967295.0); end; Function randDblExc : double; {double in (0,1)} begin randDblExc := (randInt + 0.5) * (1.0/4294967296.0); end; Function randrDblExc(range : double) : double; {double in (0,r)} begin randrDblExc := range * (randInt + 0.5) * (1.0/4294967296.0); end; Function randrInt(range : MTcard) : MTcard; {card in [0,n-1]} { Return a MTcard within [0..range-1]. The code first finds which bits are used within that range. Optimized by Magnus Jonsson (magnus(a)smartelectronix.com) } var used, i : MTcard; begin used := range; used := used OR (used shr 1); used := used OR (used shr 2); used := used OR (used shr 4); used := used OR (used shr 8); used := used OR (used shr 16); repeat i := randInt AND used {Draw numbers until one is found in [0,n]} until (i <= range); {toss unused bits to shorten search} randrInt := i; end; Function rand53 : double; {double in [0,1)} var a, b : MTcard; begin {code by Isaku Wada} a := randInt shr 5; b := randInt shr 6; rand53 := ((a * 67108864.0) + b) * (1.0/9007199254740992.0); end; Function randNorm(mean, stdev : double) : double;{double from (mean, stdev)} { Draw a double from a normal (Gaussian) distribution with specified mean and standard deviation, using the Box-Muller method. This method is adapted from Press et al.'s Numerical Recipes. Wagner used the following instead: r := sqrt(-2.0 * ln(1.0-randDblExc)) * variance; phi := 2.0 * PI * randEx0; randNorm := mean + (variance * r * cos(phi)); It is simpler, but it calls for (almost) twice as many randInts, and uses a trigonometric function. According to my time trials, the Wagner method takes approximately 3 times longer than the Box-Muller I use here. } var RF, r1, r2 : double; rsecond : double; attribute (static); bsecond : Boolean = False; attribute (static); begin if bsecond then begin bsecond := False; randNorm := rsecond; end else begin bsecond := True; repeat r1 := randr(2.0) - 1.0; r2 := randr(2.0) - 1.0; RF := sqr(r1) + sqr(r2); until ((RF < 1.0) and (RF <> 0.0)); RF := stdev * sqrt(-2.0 * ln(RF)/RF); rsecond := mean + (r2 * RF); randNorm := mean + (r1 * RF); end; end; Procedure SeedRand(RandSeed : MTcard); { Initialize generator state with seed. See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. In previous versions, the most significant bits (MSBs) of the seed affected only the MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. } begin mt[0] := RandSeed AND 16#ffffffff; for mti := 1 to mtN do mt[mti] := (1812433253 * (mt[mti-1] XOR (mt[mti-1] shr 30))+mti) AND 16#ffffffff; mti := mtN + 1; end; Procedure SeedRandArray( keylen : medCard; var initkey : PMTary); var i, j, k : MTcard; begin SeedRand(19650218); i := 1; j := 0; k := Max(mtN, keylen); repeat dec(k); mt[i] := (mt[i] XOR ((mt[i-1] XOR (mt[i-1] shr 30)) * 1664525)) + initkey^[j+1] + j; {non-linear} mt[i] := mt[i] AND 16#ffffffff; {if machine wordsize > 32} inc(i); inc(j); if (i >= mtN) then begin mt[0] := mt[mtN-1]; i := 1; end; if (j >= keylen) then j := 0; until (k = 0); for k := mtN-1 downto 0 do begin mt[i] := (mt[i] XOR ((mt[i-1] XOR (mt[i-1] shr 30)) * 1566083941)) - i; mt[i] := mt[i] AND 16#ffffffff; {if machine wordsize > 32} inc(i); if (i >= mtN) then begin mt[0] := mt[mtN-1]; i := 1; end; end; mt[0] := 16#80000000; {force MSB = 1} end; to begin do {Unit Initialization} begin {initialize state vector} SeedRand(5489); {this is Matsumoto's default} end; end. {unit MTrand} ------------------------------------------ program TestMTrand; uses MTrand, GPC; const big = 100000000; notsobig = 10000000; bin = 40; mean = 100.0; std = 20.0; type binary = array[1..bin] of longCard; var time1, time2, mult, i : longCard; m, s, rj : double; keylen, jj : medCard; bins : binary; init : PMTary; {This type is native to MTrand: PMTary = ^Buffary; Buffary(m : medCard) = array[0..m] of medCard;} begin writeln; writeln('Testing the Mersenne Twister Unit ''MTrand.pas'''); writeln; (* {This part passed March 23, 2004 in test against Matsumoto's original test output. Not implemented in this test program because of subsequent changes in the MT initialization procedures.} writeln('Precision: compare with Matsumoto''s 1000 original integer output'); SeedRand(4357); for i := 1 to 1000 do begin write(randInt:10, ' '); if ((i mod 5) = 0) then writeln; end; writeln; *) {This part passed March 24, 2004 in test against Matsumoto's output} writeln('Precision: compare with Matsumoto''s 1000 integer output, posted at'); writeln(' http://www.math.sci.hiroshima-u.ac.jp/~'); writeln(' m-mat/MT/MT2002/CODES/MTARCOK/mt19937ar-cok.out'); New(init, 4); init^[1] := 16#00000123; init^[2] := 16#00000234; init^[3] := 16#00000345; init^[4] := 16#00000456; keylen := 4; SeedRandArray(keylen, init); Dispose(init); for i := 1 to 1000 do begin write(randInt:10, ' '); if ((i mod 5) = 0) then writeln; end; writeln; {This part passed March 23, 2004} writeln('Uniformity of random 32-bit cardinals:'); for i := 1 to bin do bins[i] := 0; mult := bin; SeedRand(1234567890); for i := 1 to big do inc(bins[1+(randInt*mult) div 16#ffffffff]); for i := 1 to bin do begin write(i:2, ': ', bins[i]:7, ' '); if ((i mod 5) = 0) then writeln; end; writeln; writeln('Uniformity of range-specified cardinals:'); for i := 1 to bin do bins[i] := 0; mult := bin; SeedRand(1234567890); for i := 1 to big do inc(bins[1+(randrInt(bin))]); for i := 1 to bin do begin write(i:2, ': ', bins[i]:7, ' '); if ((i mod 5) = 0) then writeln; end; writeln; writeln; writeln('Distribution of doubles in [0..1]:'); for i := 1 to bin do bins[i] := 0; m := 0.0; s := 0.0; SeedRand(987654321); for i := 1 to notsobig do begin rj := rand; m := m + rj; s := s + sqr(rj); end; m := m/notsobig; s := sqrt((s-(notsobig*sqr(m))) / (notsobig-1)); writeln('Theoretical mean: 0.5000, standard deviation: 0.2886'); writeln('Returned mean : ', m:7:4, ', standard deviation: ', s:7:4); writeln; writeln('Distribution of Random normal deviates with known mean and standard deviation:'); for i := 1 to bin do bins[i] := 0; m := 0.0; s := 0.0; SeedRand(2468097531); for i := 1 to notsobig do begin rj := randNorm(mean, std); m := m + rj; s := s + sqr(rj); end; m := m/notsobig; s := sqrt((s-(notsobig*sqr(m))) / (notsobig-1)); writeln('Requested mean: ', mean:7:4, ', standard deviation: ', std:7:4); writeln('Returned mean: ', m:7:4, ', standard deviation: ', s:7:4); writeln; {This part tested and used to optimize the unit, March 23-24, 2004} writeln('Timing:'); SeedRand(1357908642); time1 := GetMicroSecondTime; for i := 1 to big do jj := randInt; time2 := GetMicroSecondTime; writeln('Time for ', big:1, ' randInts: ', ((time2-time1)*0.000001):7:2, ' seconds.'); SeedRand(564739201); time1 := GetMicroSecondTime; for i := 1 to big do rj := Rand; time2 := GetMicroSecondTime; writeln('Time for ', big:1, ' rands: ', ((time2-time1)*0.000001):7:2, ' seconds.'); SeedRand(1029384756); time1 := GetMicroSecondTime; for i := 1 to big do rj := randNorm(mean, std); time2 := GetMicroSecondTime; writeln('Time for ', big:1, ' randNorms: ', ((time2-time1)*0.000001):7:2, ' seconds.'); writeln; writeln('done.'); end.

2 5