IMPLEMENTATION MODULE STATUNIF; (****************************************************************************** Ce module permet a l'usager de soumettre un generateur de nombres pseudo- aleatoires a une serie de tests statistiques. Environnements : SUN et MS-Dos sur PC. -------------- Programmeurs : ------------ Pierre L'Ecuyer Pierre Audet Pierre Poulin Eric Boucher Brigitte Mercier Bruno Pelletier Richard Simard Responsable, support technique et modifications : ----------------------------------------------- Pierre L'Ecuyer Departement d'informatique et de recherche operationnelle Universite de Montreal C.P. 6128, Succ. A, Montreal, Canada H3C 3J7 ******************************************************************************) (***************************************************************************** PRENDRE CES IMPORTs POUR LE VAX FROM MYINOUT IMPORT WriteLongInt; FROM InOut IMPORT WriteLn, WriteString, WriteInt, WriteCard, WriteReal; FROM CHRONO IMPORT Chrono, InitChrono, WriteChrono, TimeFormat; FROM NUM IMPORT IToLI, IToLR, LIToLR, REALToLR, Trunc, WriteLR, RoundLR, Round, Pow, Exp, Sqrt, Log, TwoExp, LONGREAL, LONGINT, LONGCARD, OneLR, ZeroLR, TwoLR; ******************************************************************************) (***************************************************************************** PRENDRE CES IMPORTs POUR LE SUN *) FROM Files IMPORT File (*, StandardOutput*); FROM MYINOUT IMPORT WriteLongInt, WriteLn, WriteString, WriteInt, WriteCard, WriteReal, SetOutput, GetOutput, WriteLongReal ; FROM CHRONO IMPORT Chrono, InitChrono, WriteChrono, TimeFormat; FROM NUM IMPORT IToLI, IToLR, LIToLR, StringToLR, Trunc, WriteLR, RoundLR, Round, Pow, Exp, Sqrt, Ln, TwoExp, OneLR, ZeroLR, TwoLR, FourLR; FROM LongMathLib IMPORT Entier ; (*****************************************************************************) FROM SYSTEM IMPORT TSIZE, WORD, ADDRESS; FROM Storage IMPORT ALLOCATE, DEALLOCATE; FROM UNIFORM IMPORT Gener01, WriteState; (* FROM RESTRICT IMPORT (*SetSize,*) MaxIndex2, MaxIndex8; *) CONST LongTitre = 79; (* Longueur maximale du titre d'un *) (* generateur. *) Separateur = '************************************************************'; (* MaxIndexI est l'indice maximum d'un tableau d'elements de I octets, nu- *) (* merotes de [0..MaxIndexI]. *) (* Valeurs pour (*SUN*); pour (*MS-DOS*) *) MaxIndex2 = 262143 (*2^18 - 1*) (*32766*); MaxIndex4 = 262143 (*16382*); MaxIndex8 = 131071 (*8190*); (*** Sur SUN, on pourra choisir des valeurs plus grandes pour les dimensions des tableaux MaxIndex? , si on veut ****) TYPE STRING = ARRAY [0..31] OF CHAR; IVector = POINTER TO ARRAY [0..MaxIndex2] OF INTEGER; LIVector = POINTER TO ARRAY [0..MaxIndex4] OF LONGINT; RVector = POINTER TO ARRAY [0..MaxIndex8] OF LONGREAL; BSetVector = POINTER TO ARRAY [0..MaxIndex2] OF BITSET; BoolVector = POINTER TO ARRAY [0..MaxIndex2] OF BOOLEAN; (* BoolVector pourrait etre [0..65534] *) (* car les BOOLEANs sont representes *) (* sur 1 octet. (*MS-DOS*) *) VectorType = (I, LI, R, BSet, Bool); FlexTab = POINTER TO InfoFlexTab; InfoFlexTab = RECORD CASE Type : VectorType OF I : Vi : IVector | LI : Vli : LIVector | R : Vr : RVector | BSet : VBSet : BSetVector | Bool : Vbool : BoolVector END; Dimension : CARDINAL; (* Nombre d'elements du tableau *) END; FctnRepart = PROCEDURE (LONGREAL, LONGREAL) : LONGREAL; (* Fonction de repartition theorique.*) VAR Khi, (* Tableau des valeurs de Khi-Deux. *) Count, (* Compte le nombre de fois par *) (* sequence qu'un meme objet est *) (* genere. *) NbExpected, (* Nombre moyen de fois qu'un objet *) (* devrait etre genere. *) Loca : FlexTab; (* Loca[r] indique la classe corres-*) (* pondante en fonction des regroupe-*) (* ments. *) Uniforme01 : Gener01; (* Generateur que l'on teste *) TitreExperience : ARRAY [0..LongTitre] OF CHAR; C : Chrono; (* Duree d'execution du test *) LnDeux, (* ln 2 *) UnSurLnDeux, (* 1 / (ln 2) *) DixExpMoinsVingt, Stat : LONGREAL; (* Valeur de la statistique *) CorrCalculees : INTEGER; (* Taille des correl. deja calculees *) (* pour le test PeriodsInStrings *) pCorrCalculees : ADDRESS; (* pointeur sur la chaine de *) (* correlations de longueur CorrCalc.*) fichierTemp : File; (* fichier temporaire *) (****************************************************************************) PROCEDURE NotInit () : LONGREAL; BEGIN WriteString ("Erreur module STATUNIF: generateur Uniforme01 non initialise."); RETURN ZeroLR END NotInit; (****************************************************************************) PROCEDURE SetGener01 (G : Gener01; Titre : ARRAY OF CHAR); (* Cette procedure permet de definir le generateur a tester. *) (* Ce sera celui implante par la procedure G. Le parametre Titre *) (* est une chaine decrivant le generateur et apparait dans les rapports. *) VAR i, BorneSup : INTEGER ; BEGIN Uniforme01 := G; BorneSup := HIGH ( Titre ); IF (BorneSup > LongTitre) THEN BorneSup := LongTitre END; FOR i:= 0 TO BorneSup DO TitreExperience[i] := Titre[i] END; END SetGener01; (******************************************************************************) PROCEDURE MessageErreur ( Number : INTEGER (* Numero du message d'erreur a *) ); (* afficher. *) (* N.B. Apres un appel a "MessageErreur", le programme se termine et le *) (* controle retourne a l'utilisateur. *) BEGIN WriteLn; WriteLn; WriteString (Separateur); WriteLn; WriteLn; WriteString ("Erreur provenant du module TESTUNIF"); WriteLn; WriteLn; CASE Number OF 0 : WriteString ("Appel d'une fonction de repartition du Khi-Deux avec"); WriteString (" un degre de liberte inferieur a 1."); WriteLn |1 : WriteString ("Appel d'une fonction de combinaison (r dans n) avec"); WriteString (" r < 0."); WriteLn |2 : WriteString ("Demande d'allocation d'un tableau trop grand."); WriteLn |10 : WriteString ("Appel SerialTest avec NbExpected < MinExpected."); WriteLn |11 : WriteString ("Appel SerialTest avec d^k > 2^18."); (*SUN*) WriteLn; WriteLn |20 : WriteString ("Appel GapTest avec"); WriteString (" Alpha hors de l'intervalle [0..1]."); WriteLn |21 : WriteString ("Appel GapTest avec"); WriteString (" Beta hors de l'intervalle [Alpha..1]."); WriteLn |22 : WriteString ("Appel GapTest avec"); WriteString (" NbEsper < MinExpected."); WriteLn |23 : WriteString ("Appel GapTest avec"); WriteString (" t hors de l'intervalle [0..1023]."); WriteLn |30 : WriteString ("Appel SimpPoker avec d >= 128 ou k >= 128."); WriteLn |40 : WriteString ("Appel CouponCollectorTest avec t >= 50."); WriteLn |41 : WriteString ("Appel CouponCollectorTest avec d >= t ou d < 2."); WriteLn |50 : WriteString ("Appel PermutationTest avec NbExpected < MinExpected."); WriteLn |51 : WriteString ("Appel PermutationTest avec t > 8."); WriteLn (*SUN*) |60 : WriteString ("Appel SerialOver avec NbExpected < MinExpected."); WriteLn |61 : WriteString ("Appel SerialOver avec NbCases trop eleve."); WriteLn; WriteLn |62 : WriteString ("Appel SerialOver avec r + s > 31."); (*SUN*) WriteLn; WriteLn |63 : WriteString ("Appel SerialOver avec k*s > 18"); (*SUN*) WriteLn; WriteLn |70 : WriteString ("Appel MaxOft avec NbExpected < MinExpected."); WriteLn |80 : WriteString ("Appel de Collision Test avec d^k > 2^21."); (*SUN*) WriteLn |81: WriteString ("Appel de Collision Test avec n > 2^17"); (*SUN*) WriteLn |90 : WriteString ("Appel de PeriodsInStrings avec (r + s) > 31."); (*SUN*) WriteLn |91 : WriteString ("Appel de PeriodsInStrings avec s > 30 ou s < 2."); WriteLn; |100: WriteString ("Appel de NearestPair avec k > 12"); WriteLn |101: WriteString ("Appel de NearestPair avec n > 100000"); (*SUN*) WriteLn |110: WriteString ("Appel de Maurer avec L non multiple de s "); WriteLn |111: WriteString ("Appel de Maurer avec r + s > 31 "); WriteLn (*SUN*) |113: WriteString ("Appel de Maurer avec L > 18"); WriteLn (*SUN*) |120: WriteString ("Appel Opso avec r < 0 ou r > 20."); WriteLn |121: WriteString ("Appel Opso avec p non dans {1, 2, 3}"); WriteLn (*SUN*) |130: WriteString ("Appel LongestHeadRun avec r + s > 31."); (*SUN*) WriteLn; |131: WriteString ("Erreur dans LongestHeadRun: implantation du test"); WriteString ("incompatible avec l'adressage des elements d'un"); WriteString ("ensemble dans cet ordinateur -> transfert de type"); WriteLn |132: WriteString ("Appel LongestHeadRun avec n trop petit."); WriteLn |133: WriteString ("Appel LongestHeadRun avec l < 1000"); WriteLn |140: WriteString ("Appel gtransform avec g < 1.0 ou g > 10.0"); WriteLn |152: WriteString ("Appel BirthdaySpacings avec k < 20."); WriteLn |150: WriteString ("Appel BirthdaySpacings avec m trop grand."); WriteLn |151: WriteString ("Appel BirthdaySpacings avec m^3/n trop grand.");WriteLn END; WriteLn; WriteString ("Interruption du test."); WriteLn; WriteLn;WriteString (Separateur); WriteLn; HALT (* Termine le programme *) END MessageErreur; (****************************************************************************) PROCEDURE NewFlexTab ( VAR F : FlexTab; T : VectorType); (* Cette procedure cree une variable de type InfoFlexTab et affecte *) (* l'adresse de cette nouvelle variable au parametre F. Le parametre T *) (* indique le type d'elements que contient le tableau de la structure *) (* InfoFlexTab nouvellement cree. Apres un appel a NewFlexTab, ce tableau *) (* n'existe pas. Il sera cree, en temps opportun, par la procedure *) (* CreateTab. *) BEGIN NEW (F); WITH F^ DO Type := T; Dimension := 0 END END NewFlexTab; (****************************************************************************) PROCEDURE DeleteTab ( VAR Tab : FlexTab); (* Detruit le tableau de la structure InfoFlexTab pointee par Tab. *) BEGIN WITH Tab^ DO CASE Type OF I : DEALLOCATE (Vi, Dimension * TSIZE (INTEGER)) | LI : DEALLOCATE (Vli, Dimension * TSIZE (LONGINT)) | R : DEALLOCATE (Vr, Dimension * TSIZE (LONGREAL)) | BSet : DEALLOCATE (VBSet, Dimension * TSIZE (BITSET)) | Bool : DEALLOCATE (Vbool, Dimension * TSIZE(BOOLEAN)) END; Dimension := 0 END END DeleteTab; (****************************************************************************) PROCEDURE CreateTab ( NbElement : INTEGER; VAR Tableau : FlexTab ); (* Procedure appelee pour creer un tableau necessaire a un test. Le *) (* parametre NbElement indique le nombre d'elements du tableau. Ce nombre *) (* ne peut pas depasser MaxIndex2, MaxIndex4 ou MaxIndex8 selon le cas. *) BEGIN WITH Tableau^ DO (* Verifie la taille du tableau *) IF ( (Type = R) AND (NbElement > (MaxIndex8 + 1)) ) OR ( (Type = LI) AND (NbElement > (MaxIndex4 + 1)) ) OR ( (Type # R) AND (Type # LI) AND (NbElement > (MaxIndex2 + 1)) ) THEN MessageErreur (2); END; IF (NbElement = 0) THEN (* Cas possible pour un tableau de *) Dimension := 0 (* BITSET *) ELSE Dimension := NbElement END; CASE Type OF I : ALLOCATE (Vi, Dimension * TSIZE (INTEGER)) | LI : ALLOCATE (Vli, Dimension * TSIZE (LONGINT)) | R : ALLOCATE (Vr, Dimension * TSIZE (LONGREAL)) | BSet : ALLOCATE (VBSet, Dimension * TSIZE (BITSET)) | Bool : ALLOCATE (Vbool, Dimension * TSIZE (BOOLEAN)) END END END CreateTab; (****************************************************************************) PROCEDURE Entete (Test : STRING); BEGIN InitChrono (C); WriteLn; WriteString (Separateur); WriteLn; WriteLn; WriteString (Test); WriteString (' :'); WriteLn; WriteString ('-----------------------'); WriteLn; WriteString (TitreExperience); WriteLn; WriteLn END Entete; (****************************************************************************) PROCEDURE WriteResults (Stat, KS : LONGREAL); BEGIN WriteLn; WriteLn; WriteString ('Valeur de la statistique de K-S : '); WriteLR (Stat, 6, 4, 2); WriteLn; WriteString ('Niveau descriptif du test : '); WriteLR (KS, 6, 4, 2); WriteLn; WriteString ('Temps CPU consomme : '); WriteChrono (C, hms); WriteLn; WriteState; END WriteResults; (****************************************************************************) PROCEDURE Factorielle (N : INTEGER) : LONGREAL; (* Cette fonction retourne la factorielle de N. *) VAR T : LONGREAL; C : INTEGER; BEGIN T := OneLR; FOR C := 1 TO N DO T := T * IToLR (C) END; RETURN T; END Factorielle; (****************************************************************************) PROCEDURE Combinaison (* / n \ n! *) ( r, (* \ r / = ------------------ *) n : INTEGER (* r! (n-r)! *) ) : LONGREAL; (* *) (* *) (* Si r < (n - r) : *) (* *) (* n(n-1) ... (n-r+1) *) (* = ------------------ *) (* r(r-1) ... 1 *) (* *) (* *) (* Si r >= (n - r) : *) (* *) (* n(n-1) ... (r+1) *) (* = ------------------ *) (* (n-r) (n-r-1)... 1 *) (* Retourne la valeur d'une combinaison de r dans n valeurs; i.e. le *) (* nombre de facons de choisir r valeurs parmi n. *) VAR Res : LONGREAL; Diff, i : INTEGER; BEGIN IF r = 0 THEN RETURN OneLR END; IF (r < 0) THEN MessageErreur (1) ELSE IF (r > n) THEN RETURN ZeroLR; ELSE Res := OneLR; Diff := n - r; IF (r < Diff) THEN FOR i := 1 TO r DO Res := Res * IToLR (Diff + i) / IToLR (i) END; ELSE FOR i := 1 TO Diff DO Res := Res * IToLR (r + i) / IToLR (i) END; END; RETURN Res; END END END Combinaison; (****************************************************************************) (* les constantes dans le calcul de FctnRepartNormale *) TYPE RecordFRNormale = RECORD p10, p11, p12, p13, p20, p21, p22, p23, p24, p25, p26, p27, p30, p31, p32, p33, p34, q10, q11, q12, q13, q20, q21, q22, q23, q24, q25, q26, q27, q30, q31, q32, q33, q34, RacineUnSurPi, RacineDeDeux, POINT477, DEUX, QUATRE : LONGREAL END; ptrToFRNormale = POINTER TO RecordFRNormale; VAR DataFRNormale : ptrToFRNormale; (****************************************************************************) (* pour initialiser les constantes de la FctnRepartNormale *) PROCEDURE InitFRNormale ( VAR DataFRNormale:ptrToFRNormale); BEGIN WITH DataFRNormale^ DO p10 := StringToLR ("2.4266795523053175E2"); p11 := StringToLR ("2.1979261618294152E1"); p12 := StringToLR ("6.9963834886191355E0"); p13 := StringToLR ("-3.5609843701815385E-2"); p20 := StringToLR ("3.004592610201616005E2"); p21 := StringToLR ("4.519189537118729422E2"); p22 := StringToLR ("3.393208167343436870E2"); p23 := StringToLR ("1.529892850469404039E2"); p24 := StringToLR ("4.316222722205673530E1"); p25 := StringToLR ("7.211758250883093659E0"); p26 := StringToLR ("5.641955174789739711E-1"); p27 := StringToLR ("-1.368648573827167067E-7"); p30 := StringToLR ("-2.99610707703542174E-3"); p31 := StringToLR ("-4.94730910623250734E-2"); p32 := StringToLR ("-2.26956593539686930E-1"); p33 := StringToLR ("-2.78661308609647788E-1"); p34 := StringToLR ("-2.23192459734184686E-2"); q10 := StringToLR ("2.1505887586986120E2"); q11 := StringToLR ("9.1164905404514901E1"); q12 := StringToLR ("1.5082797630407787E1"); q13 := StringToLR ("1.0"); q20 := StringToLR ("3.004592609569832933E2"); q21 := StringToLR ("7.909509253278980272E2"); q22 := StringToLR ("9.313540948506096211E2"); q23 := StringToLR ("6.389802644656311665E2"); q24 := StringToLR ("2.775854447439876434E2"); q25 := StringToLR ("7.700015293522947295E1"); q26 := StringToLR ("1.278272731962942351E1"); q27 := StringToLR ("1.0"); q30 := StringToLR ("1.06209230528467918E-2"); q31 := StringToLR ("1.91308926107829841E-1"); q32 := StringToLR ("1.05167510706793207E0"); q33 := StringToLR ("1.98733201817135256E0"); q34 := StringToLR ("1.0"); RacineUnSurPi := StringToLR ("0.56418958354775628694E0"); RacineDeDeux := StringToLR ("1.4142135623730950488E0"); POINT477 := StringToLR ("0.477"); DEUX := StringToLR ("2.0"); QUATRE := StringToLR ("4.0"); END; END InitFRNormale; (****************************************************************************) PROCEDURE FctnRepartNormale ( Bidon, (* Pour respecter le type FctnRepart *) X : LONGREAL ) : LONGREAL; (* Fonction de repartition d'une v.a. suivant une loi normale N(0,1). *) (* Ref. : Kennedy and Gentle, p.90 *) (* La fonction de repartition de la loi normale est utilisee pour appro- *) (* ximer la fonction de repartition d'une Khi-deux. *) VAR YCarre, UnSurY2, Y, R, Erf : LONGREAL; BEGIN IF (X < ZeroLR) THEN RETURN ( OneLR - FctnRepartNormale (ZeroLR, -X)) END; (* Cas ou X >= 0 *) WITH DataFRNormale^ DO Y := X / RacineDeDeux; YCarre := X * X / DEUX; IF (Y <= POINT477) THEN (* Cas ou 0 <= X <= 0.477 *) R := (p10 + YCarre * (p11 + YCarre * (p12 + YCarre * p13))) / (q10 + YCarre * (q11 + YCarre * (q12 + YCarre * q13))); Erf := Y * R ELSIF (Y <= QUATRE) THEN (* Cas ou 0.477 <= X <= 4.0 *) R := (p20 + Y * (p21 + Y * (p22 + Y * (p23 + Y * (p24 + Y * (p25 + Y * (p26 + Y * p27))))))) / (q20 + Y * (q21 + Y * (q22 + Y * (q23 + Y * (q24 + Y * (q25 + Y * (q26 + Y * q27))))))); Erf := OneLR - Exp (- YCarre) * R ELSE (* Cas ou X > 4.0 *) UnSurY2 := OneLR / YCarre; R := (p30 + UnSurY2 * (p31 + UnSurY2 * (p32 + UnSurY2 * (p33 + UnSurY2 * p34)))) / (q30 + UnSurY2 * (q31 + UnSurY2 * (q32 + UnSurY2 * (q33 + UnSurY2 * q34)))); Erf := OneLR - (Exp (- YCarre) / Y) * (RacineUnSurPi + R / YCarre) END; END; RETURN ((OneLR + Erf) / TwoLR) END FctnRepartNormale; (****************************************************************************) (* les constantes dans le calcul de FctnRepartKhiDeux *) TYPE RecordFRKhiDeux = RECORD Pt2, TIERS, MOINSDIXHUIT, CINQ, Gam :LONGREAL; END; ptrToFRKhiDeux = POINTER TO RecordFRKhiDeux; VAR DataFRKhiDeux : ptrToFRKhiDeux; (****************************************************************************) (* pour initialiser les constantes de la FctnRepartKhiDeux *) PROCEDURE InitFRKhiDeux ( VAR DataFRKhiDeux :ptrToFRKhiDeux); BEGIN WITH DataFRKhiDeux^ DO Pt2 := StringToLR ("0.22222222"); TIERS := StringToLR ("0.33333333"); CINQ := StringToLR ("5.0"); MOINSDIXHUIT := StringToLR ("-18.8055"); Gam := StringToLR ("0.8862269254527579825931E0"); END; (* Racine de PI, / 2 = Gamma(3/2) *) END InitFRKhiDeux; (******************************************************************************) PROCEDURE FctnRepartKhiDeux ( DegLib, (* Deg_lib est entier *) X : LONGREAL ) : LONGREAL; (* Valeur de la fonction de repartition d'une Khi-Deux a "DegLib" degres *) (* de liberte. Ref. : Kennedy and Gentle, p.116 *) (* La fonction de repartition de la Khi-deux est approximee par l'utilisa- *) (* tion de la fonction de repartition de la loi normale. *) VAR H, H2, E, DemiX, Terme, Sommation : LONGREAL; DL, n : LONGINT; BEGIN DL := Trunc (DegLib); IF (DL > 66) THEN (* Approx. par la loi normale lorsque*) (* le degre de liberte est grand *) IF (X < TwoLR) THEN RETURN ZeroLR END; X := (Pow ((X / DegLib), DataFRKhiDeux^.TIERS) - (OneLR - DataFRKhiDeux^.Pt2 / DegLib)) / Sqrt (DataFRKhiDeux^.Pt2 / DegLib); IF (X > DataFRKhiDeux^.CINQ) THEN RETURN OneLR END; IF (X < DataFRKhiDeux^.MOINSDIXHUIT) THEN RETURN ZeroLR END; RETURN FctnRepartNormale (ZeroLR, X) ELSE DemiX := X / TwoLR ; IF NOT (ODD (DL)) THEN (* Cas ou degre de liberte est pair *) Terme := Exp (-DemiX); Sommation := Terme; FOR n := 1 TO (DL DIV 2) - 1 DO Terme := Terme * DemiX / LIToLR (n); Sommation := Sommation + Terme; END; RETURN OneLR - Sommation ELSE (* Cas ou degre de liberte est impair*) H2 := -OneLR + TwoLR * FctnRepartNormale (ZeroLR, Sqrt(X)); IF DL = 1 THEN RETURN (H2) END; E := Exp (-DemiX); Terme := Sqrt (DemiX) * E / DataFRKhiDeux^.Gam; n := 3; H := H2; WHILE (n < DL) DO H := H - Terme; Terme := Terme * DemiX * TwoLR / LIToLR (n); n := n + 2; END; RETURN (H - Terme); END END END FctnRepartKhiDeux; (***************************************************************************) (* les constantes dans le calcul de FctnRepartKS *) TYPE RecordFRKS = RECORD (* constantes de FctnRepartKS *) P1, P2, P3, P4, P5, POINT571, POINT520, POINT481, POINT441, POINT10, POINT05, POINT025, POINT01, UnPointTrois, DeuxPointSix, MILLE, DEUXCINQUANTE, MOINS4, POINTQUINZE, POINT275, POINT04 : LONGREAL; END; ptrToFRKS = POINTER TO RecordFRKS; VAR DataFRKS : ptrToFRKS; (****************************************************************************) (* pour initialiser les constantes de la FctnRepartKS *) PROCEDURE InitFRKS ( VAR DataFRKS:ptrToFRKS); BEGIN WITH DataFRKS^ DO P1 := StringToLR ("1.282051282"); P2 := StringToLR ("0.980392156"); P3 := StringToLR ("0.641025641"); P4 := StringToLR ("0.375000000"); P5 := StringToLR ("0.022675736"); POINT571 := StringToLR ("0.571"); POINT520 := StringToLR ("0.520"); POINT481 := StringToLR ("0.481"); POINT441 := StringToLR ("0.441"); UnPointTrois := StringToLR ("1.33333333"); DeuxPointSix := StringToLR ("2.66666667"); POINT10 := StringToLR ("0.1"); POINT05 := StringToLR ("0.05"); POINT025 := StringToLR ("0.025"); POINT01 := StringToLR ("0.01"); MILLE := StringToLR ("1000.0"); DEUXCINQUANTE := StringToLR ("250.0"); MOINS4 := StringToLR ("-4.0"); POINTQUINZE := StringToLR ("0.15"); POINT275 := StringToLR ("0.275"); POINT04 := StringToLR ("0.04"); END; END InitFRKS; (****************************************************************************) PROCEDURE FctnRepartKS ( n, (* Nombre de valeurs *) X : LONGREAL ) : LONGREAL; (* Approximation de la fonction de repartition de la statistique *) (* de Kolmogorov-Smirnov : Dn = Max (Fn - F) ou Fn et F sont respecti- *) (* vement la fonction de repartition empirique d'une v.a. et la fonc- *) (* tion de repartion theorique de la v.a. . On utilise le corollaire *) (* Z de D.A. Darling, dans "Theory of Probability and its Applications" *) (* 1960, vol. V, 4, p.356. L'approximation devient meilleure lorsque *) (* n augmente ou encore lorsque X s'eloigne de zero. Lorsque la valeur *) (* FctnRepartKS est < 0.15, on utilise le tableau "Lower Tail" de *) (* Stephens, Use of the Kolmogorov-Smirnov, Cramer - von Mises and *) (* Related Statistics without Extensive Tables, J. of the Royal Statis- *) (* tical Society, Series B, 33, 1, 115-122, 1970, avec interpolation *) (* lineaire. L'erreur sur la valeur de la fonction ne depasse pas 1% *) (* (environ) lorsque n > 100. *) VAR T, E2, E8, Resultat, RacineCarreen : LONGREAL; PROCEDURE LowerTail (Y : LONGREAL) : LONGREAL; BEGIN WITH DataFRKS^ DO IF (Y > POINT571) THEN RETURN ( POINT10 + ((Y - POINT571) * P1 )) ELSIF (Y > POINT520) THEN RETURN ( POINT05 + ((Y - POINT520) * P2 )) ELSIF (Y > POINT481) THEN RETURN ( POINT025 + ((Y - POINT481) * P3 )) ELSIF ( Y > POINT441) THEN RETURN ( POINT01 + ((Y - POINT441) * P4 )) ELSE RETURN (Y * P5) END; END END LowerTail; BEGIN T := TwoLR * n * X * X; (* Commentaire concernant la version sur PC: *) (* La fonction MATHLIB.Exp a un bug: pour des valeurs de -1500 a -11400 *) (* elle retourne un NaN (Not a Number). D'ou les conditions suivantes: *) IF T > DataFRKS^.MILLE THEN E2 := ZeroLR ELSE E2 := Exp (-T) END; IF T > DataFRKS^.DEUXCINQUANTE THEN E8 := ZeroLR ELSE E8 := Exp (DataFRKS^.MOINS4 * T) END; Resultat := OneLR + (DataFRKS^.UnPointTrois * X - TwoLR) * E2 + (TwoLR - DataFRKS^.DeuxPointSix * X) * E8; IF (Resultat < DataFRKS^.POINTQUINZE) THEN RacineCarreen := Sqrt(n); Resultat := LowerTail (X * (RacineCarreen + DataFRKS^.POINT275 - DataFRKS^.POINT04 / RacineCarreen)) END; RETURN Resultat END FctnRepartKS; (****************************************************************************) PROCEDURE StatKS ( Fctn : FctnRepart; (* Fonction de repartition theorique *) (* d'une v.a. quelconque. (Dans la *) (* pluspart des cas, Khi-Deux). *) Param1 : LONGREAL; (* Si Fctn est Khi-deux alors "Param1" *) (* represente le degre de liberte. *) VAR V : FlexTab; (* Differentes observations d'une v.a. *) (* a partir desquelles on va construire*) (* la fonction de repartition empirique*) IndiceFin : INTEGER (* Nombre d'observations considerees *) (* dans le vecteur V, de l'indice 1 a *) (* l'indice IndiceFin. *) ) : LONGREAL; (* Retourne la valeur de la statistique qui correspond a un test de K-S *) (* pour comparer la fonction de repartition theorique FctnRepart *) (* avec la fonction de repartition empirique des valeurs contenues *) (* dans le vecteur V. *) VAR FOdeX, D, (* Valeur de la statistique de K-S *) D1, D2, (* Statistiques de K-S intermediaires *) UnSurIndiceFin : LONGREAL; i : INTEGER; PROCEDURE QuickSort ( VAR Tab : FlexTab; n : INTEGER (* Dernier indice considere dans le *) (* du tri du tableau Tab. *) ); (* On trie le tableau des observations de la v.a. pour faciliter *) (* l'obtention de la fonction de distribution empirique. *) PROCEDURE Sort (l, r : INTEGER); (* QuickSort recursif tel que presente par Niklaus Wirth, *) (* dans Algorithms + Data Structures = Programs, p.79. *) VAR i, j : INTEGER; (* Indices dans le tableau Tab. *) x, w : LONGREAL; BEGIN (* Sort *) i := l; j := r; x := Tab^.Vr^[(l+r) DIV 2]; REPEAT WHILE (Tab^.Vr^[i] < x) DO i := i + 1 END; WHILE (x < Tab^.Vr^[j]) DO j := j - 1 END; IF (i <= j) THEN w := Tab^.Vr^[i]; Tab^.Vr^[i] := Tab^.Vr^[j]; Tab^.Vr^[j] := w; i := i + 1; j := j - 1 END; UNTIL (i > j); IF (l < j) THEN Sort (l, j) END; IF (i < r) THEN Sort (i, r) END END Sort; BEGIN (* QuickSort *) Sort (1, n) END QuickSort; BEGIN (* StatKS *) (* Logitech ne permet pas de verifier l'egalite Fcnt = FctnRepartKhiDeux *) IF ( (TRUNC (Param1) < 1) (* AND (Fctn = FctnRepartKhiDeux) *) ) THEN MessageErreur (0) END; QuickSort (V, IndiceFin); D := ZeroLR; UnSurIndiceFin := OneLR / IToLR (IndiceFin); (* Recherche de D, la mesure de la plus grande divergence entre les *) (* fonctions de repartition theorique et empirique. *) FOR i := 1 TO IndiceFin DO FOdeX := Fctn ( Param1, (V^.Vr^[i]) ); D1 := ABS (FOdeX - IToLR (i - 1) * UnSurIndiceFin); D2 := ABS (FOdeX - IToLR (i) * UnSurIndiceFin); IF (D1 > D) THEN D := D1 END; IF (D2 > D) THEN D := D2 END END; RETURN D END StatKS; (****************************************************************************) PROCEDURE Serial ( N : INTEGER; n : LONGINT; d, k : INTEGER ); VAR KhiSeq, (* Variable tampon utilisee pour sau-*) (* ver du temps d'acces dans Khi[Seq]*) dReal, (* Valeur de d convertie en reel *) kReal, (* Valeur de k convertie en reel *) NbExpected, (* Nombre moyen de fois qu'un meme *) (* vecteur devrait etre genere. *) (* (Nombre espere d'observations par *) (* case). *) Diff : LONGREAL; j, l, (* Nombre courant de vecteurs generes*) (* dans une sequence *) Base, (* Compteur permettant la conversion *) (* d'un indice tableau en indice *) (* vecteur *) Seq, (* Numero de la sequence courante *) Indice, (* Indice du tableau converti en *) (* indice du vecteur *) NbCases : INTEGER; (* Nombre de cases contenues dans le *) (* tableau defini par d ** Dim. *) L : LONGINT; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput ( StdOut ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Serial (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); WriteInt(d,1); WriteString(', '); WriteInt(k,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); BEGIN dReal := IToLR (d); kReal := IToLR (k); NbCases := Entier (RoundLR (Pow (dReal, kReal))); NbExpected := LIToLR (n) / IToLR (NbCases); END ProcInit; BEGIN ProcInit; (* Imprime l'entete du rapport *) Entete ('Serial Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', n = '); WriteLongInt (n, 4); WriteString (', d = '); WriteInt (d, 4); WriteString (', k = '); WriteInt (k, 2); WriteLn; WriteString (' Nb. de cases = d ** k = '); WriteInt (NbCases, 4); WriteLn; WriteString (' Nb. espere par case = '); WriteLR (NbExpected, 7, 1, 2); WriteLn; WriteLn; IF (NbCases > (MaxIndex2 + 1)) THEN MessageErreur (11) END; IF (NbExpected < MinExpected) THEN MessageErreur (10) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (NbCases, Count); (* Dimension de Count [0..NbCase-1] *) IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR l := 0 TO (NbCases - 1) DO Count^.Vi^[l] := 0 END; FOR L := 1 TO n DO (* Generation d'un vecteur *) Indice := TRUNC (Uniforme01() * dReal); Base := 1; FOR j := 1 TO (k - 1) DO Base := Base * d; Indice := Indice + Entier (Uniforme01() * dReal) * Base END; INC (Count^.Vi^[Indice]) END; (* Calcul de la statistique du *) (* Khi-Deux pour une sequence. *) KhiSeq := ZeroLR; FOR l := 0 TO (NbCases - 1) DO Diff := IToLR (Count^.Vi^[l]) - NbExpected; KhiSeq := KhiSeq + Diff * Diff END; KhiSeq := KhiSeq / NbExpected; Khi^.Vr^[Seq] := KhiSeq; IF Detailed THEN WriteLR (KhiSeq, 7, 2, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultat *) END; END END; (* Test de K-S sur la repartition des*) (* valeurs de Khi[Seq] *) Stat := StatKS (FctnRepartKhiDeux, IToLR(NbCases - 1), Khi, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (Count); END Serial; (******************************************************************************) PROCEDURE SerialOver ( N : INTEGER; n : LONGINT; r, s, k : CARDINAL ); (* 1) On genere les coordonnees de points dans un espace a k dimensions par concatenation de groupes de s bits tires de nombres produits par le gene- rateur de v.a. a tester. Pour un point donne, les s msb (k-ieme coordonnee ) proviennent du (k-1) ieme nombre precedent produit par le generateur de v.a. et les s lsb proviennent du nombre courant produit par le generateur. 2) Pour une implantation donnee (SUN ou MS-DOS), s et k sont etroitement relies et limites par le nombre de bits representant l'indice maximum d'un tableau d'elements d'entiers. Par exemple sur la version PC cet in- dice maximum est 32766 (tableau indice a partir de 0); la relation entre s et k est la suivante: s * k <= plancher( log2(32766) ). - s represente le nombre de bits qu'on observe sur un nombre produit par un generateur de v.a. - s * k est le nombre de bits representant un point dans l'espace a k dimensions. - plancher( log2(32766) ) = 14 est le nombre de bits du maximum de points d'un espace a une ou plusieurs dimensions avec MS-DOS. 3) Sur une implantation ou la dimension des tableaux n'est pas un probleme (i.e. une implant. non MS-DOS), la nouvelle restriction est le nombre de bits de MAX( INTEGER ), puisque NbCasesK numerote les vecteurs, ou bien la dimension maximale d'un tableau: MaxIndex2 + 1. 4) La facon particuliere dont est calcule le khi-deux (etant donne la cor- relation due au recouvrement) est expliquee dans: - Altman, N.S. (1988). Bit-Wise Behavior Of Random Number Generators. Siam J. Sci. Stat. Comput., Vol.9, No 5, 941-949. - Marsaglia, G. (1985). A Current View of Random Number Generation. Computer Science and Statistics, Proceedings of the Sixteenth Symposium on the Interface, Atlanta. march 1484. Elsevier Publ. (North-Holland), 3-10. *) VAR KhiSeqK, KhiSeqKm1, NbExpectedK, (* Nombre moyen de fois qu'un meme *) (* k tuple devrait etre genere. *) NbExpectedKm1, (* Nombre moyen de fois qu'un meme *) (* (k - 1) tuple devrait etre genere. *) Base, (* Base := 2^(r + s) *) DiffKm1, (* DiffKm1 := Nbre obs. de (k - 1) tu- ples - nbre espere. *) DiffK : LONGREAL; (* DiffK := Nbre obs. de k tuples - nbre espere. *) ii, (* Nombre courant de vecteurs generes dans une sequence *) Seq, (* Numero de la sequence courante *) NbCasesKm1, (* Nombre possible de (k-1)tuples diff.*) NbCasesK, (* Nombre possible de k tuples diff. *) DegLib : INTEGER; (* Degres de liberte. *) kLInt, (* k en LONGINT. *) L : LONGINT; Maskm1, (* Masque qui selectionne les s*(k - 1)*) (* premiers bits. *) MaskS, (* Masque qui selectionne les s pre- *) (* miers bits. *) j, element, (* coordonnee d'un point *) point : CARDINAL; (* point a k dimensions (elements) *) longMaskS, (* MaskS en LONGCARD *) rsBits : LONGCARD; (* element en LONGCARD *) (* Tableau des (k - 1) premiers elements du premier tuple. Ces tuples *) (* sont necessaires pour que la sequence de nombres soit circulaire. *) Premier : ARRAY[1..16] OF CARDINAL; PROCEDURE ProcInit (); BEGIN kLInt := k; MaskS := TRUNC (TwoExp [s]); Maskm1 := TRUNC (TwoExp [s * (k - 1)]); longMaskS := VAL( LONGCARD, MaskS); Base := TwoExp [r + s]; NbCasesK := Entier (TwoExp [s * k]); NbCasesKm1 := Entier (TwoExp [s * (k - 1)]); NbExpectedK := LIToLR (n) / IToLR (NbCasesK); NbExpectedKm1 := LIToLR (n) / IToLR (NbCasesKm1); DegLib := NbCasesK - NbCasesKm1; END ProcInit; PROCEDURE Moniteur(); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' SerialOver (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); WriteCard(r,1); WriteString(', '); WriteCard(s,1); WriteString(', '); WriteCard(k,1); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; BEGIN Entete ('SerialOver Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString ('; n = '); WriteLongInt (n, 4); WriteString ('; r = '); WriteCard (r, 2); WriteString ('; s = '); WriteCard (s, 2); WriteString ('; k = '); WriteCard (k, 2); WriteLn; IF (r + s) > 31 THEN MessageErreur (62) END; (*SUN*) IF s*k > 18 THEN MessageErreur (63) END; (*SUN*) ProcInit (); WriteString ('Nb. de k tuples = 2^(s*k) = '); WriteInt (NbCasesK, 4); WriteLn; WriteString ('Nb. espere par tuple = '); WriteLR (NbExpectedK, 7, 1, 2); WriteLn; WriteLn; IF (NbCasesK > (MaxIndex2 + 1)) THEN MessageErreur (61) END; IF (NbExpectedK < MinExpected) THEN MessageErreur (60) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (NbCasesK, Count); (* Dimension de Count [0..NbCase-1] *) CreateTab (NbCasesKm1, Loca); IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; (* Initialisation des tableaux dans lesquels on compte l'occurence des k-tuples (tableau Count) et des (k-1)-tuples (tableau Loca). *) WITH Count^ DO FOR ii := 0 TO NbCasesK - 1 DO Vi^[ii] := 0 END END; WITH Loca^ DO FOR ii := 0 TO NbCasesKm1 - 1 DO Vi^[ii] := 0 END END; (* Generation des (k - 1) premiers elements du premier tuple. On les garde dans le tableau Premier[] pour que la sequence des nombres generes soit circulaire. Ils seront utilises pour obtenir les k-1 derniers tuples. Pour un nombre aleatoire donne, l'element obtenu des bits que le test observe devient l'element le moins significatif du tuple. A chaque nombre aleatoire genere les elements du tuple sont decales (l'element le plus significatif est perdu). *) (* Gener. des (k-1) premiers elements du premier tuple *) point := 0; FOR j := 1 TO (k - 1) DO (* Generation d'un nombre aleatoire; on garde seulement les (r + s) premiers bits *) rsBits := TRUNC (Base * Uniforme01 ()); (* elimination des r premiers bits de gauche *) element := VAL (CARDINAL, (rsBits MOD longMaskS)); Premier[j] := element; (* Decalage du tuple de s bits et insertion du nouvel element *) point := point * MaskS + element; END; (* Generation des n - (k - 1) premiers tuples. *) FOR L := 1 TO (n - (kLInt - 1)) DO (* le (k - 1)-tuple est defini par (k - 1)*s bits *) point := point MOD Maskm1; INC (Loca^.Vi^[point]); (* Generation d'un nombre aleatoire et decalage de (r + s) bits *) rsBits := TRUNC (Base * Uniforme01 ()); (* elimination des r premiers bits de gauche *) element := VAL (CARDINAL, (rsBits MOD longMaskS)); (* Decalage du tuple de s bits et insertion du nouvel element *) point := point * MaskS + element; INC (Count^.Vi^[point]); END; (* Generation des (k - 1) derniers tuples. On utilise les nombres du tableau Premier[] pour que la sequence soit circulaire. *) FOR j := 1 TO (k - 1) DO point := point MOD Maskm1; INC (Loca^.Vi^[point]); element := Premier[j]; point := point * MaskS + element; INC (Count^.Vi^[point]); END; (* Calcul de la statistique du Khi-Deux pour une sequence. *) KhiSeqK := ZeroLR; KhiSeqKm1 := ZeroLR; FOR ii := 0 TO NbCasesK - 1 DO DiffK := IToLR (Count^.Vi^[ii]) - NbExpectedK; KhiSeqK := KhiSeqK + (DiffK * DiffK) END; FOR ii := 0 TO NbCasesKm1 - 1 DO DiffKm1 := IToLR (Loca^.Vi^[ii]) - NbExpectedKm1; KhiSeqKm1 := KhiSeqKm1 + (DiffKm1 * DiffKm1) END; Khi^.Vr^[Seq] := (KhiSeqK / NbExpectedK) - (KhiSeqKm1 / NbExpectedKm1); IF Detailed THEN WriteLR (Khi^.Vr^[Seq], 7, 3, 3); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultat *) END END END; (* Test de K-S sur la repartition des valeurs de Khi[Seq]. *) Stat := StatKS (FctnRepartKhiDeux, IToLR (DegLib), Khi, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (Count); DeleteTab (Loca) END SerialOver; (******************************************************************************) PROCEDURE Gap ( N, n, t : INTEGER; Alpha, Beta : LONGREAL ); VAR r, (* Nombre courant de valeurs qui *) (* sont generees. *) s, (* Nombre courant de Gaps observes *) Seq : INTEGER; (* Numero de la sequence courante. *) Mult, p, (* Probabilite qu'une valeur se re- *) (* trouve entre Alpha et Beta. *) U, KhiSeq, (* Variable tampon utilisee pour sau- *) (* ver du temps d'acces dans Khi[Seq] *) Restriction, (* Min ( np (1-p)**(t-1), n(1-p)**t *) Restriction1, (* n * p * (1 - p) ** (t - 1) *) Restriction2, (* n * (1 - p) ** (t - 1) *) Diff : LONGREAL; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Gap (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(t,1); WriteString(', '); WriteLR (Alpha, 6, 2, 4); WriteString(', '); WriteLR (Beta, 6, 2, 4); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit; BEGIN p := Beta - Alpha; Restriction1 := IToLR (n) * p * Pow (OneLR - p, IToLR (t - 1)); Restriction2 := IToLR (n) * Pow (OneLR - p, IToLR (t)); IF (Restriction1 < Restriction2) THEN Restriction := Restriction1 ELSE Restriction := Restriction2 END; Mult := p * IToLR (n); END ProcInit; BEGIN ProcInit; (* Imprime l'entete du rapport *) Entete ('Gap Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', Alpha = '); WriteLR (Alpha, 6, 4, 2); WriteString (', Beta = ') ; WriteLR (Beta, 6, 4, 2); WriteLn; WriteString (' n = ') ; WriteInt (n, 4); WriteString (' et t = ') ; WriteInt (t, 4); WriteLn; WriteString ('Nb. espere par case = '); WriteLR (Restriction, 7, 1, 2); WriteLn; WriteLn; IF ((Alpha < ZeroLR) OR (Alpha > OneLR)) THEN MessageErreur (20) END; IF ((Beta > OneLR) OR (Beta <= Alpha)) THEN MessageErreur (21) END; IF (t <= 0) OR (t >= 1023) THEN MessageErreur (23) END; IF (Restriction < MinExpected) THEN MessageErreur (22) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (t+1, NbExpected); (* Dimension de NbExpected [0..t] *) CreateTab (t+1, Count); (* Dimension de Count [0..t] *) IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; NbExpected^.Vr^[0] := Mult; (* Initialisation; calcul des proba- *) FOR r := 1 TO (t - 1) DO (* bilites pour chaque longueur de *) Mult := Mult * (OneLR - p); (* 'gap'. *) NbExpected^.Vr^[r] := Mult END; NbExpected^.Vr^[t] := Mult * (OneLR - p) / p; (* Debut du test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR r := 0 TO t DO Count^.Vi^[r] := 0 END; FOR s := 1 TO n DO (* Traitement d'un 'gap' *) r := 0; U := Uniforme01 (); WHILE (U < Alpha) OR (U > Beta) DO INC (r); U := Uniforme01 () END; IF (r >= t) THEN INC ( Count^.Vi^[t]) ELSE INC( Count^.Vi^[r]) END END; (* Calcul du Khi-Deux pour une seq. *) KhiSeq := ZeroLR; FOR r := 0 TO t DO Diff := IToLR(Count^.Vi^[r]) - NbExpected^.Vr^[r]; KhiSeq := KhiSeq + (Diff * Diff) / NbExpected^.Vr^[r] END; IF Detailed THEN WriteLR (KhiSeq, 7, 4, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultats*) END; END; Khi^.Vr^[Seq] := KhiSeq; END; (* Test de K-S sur la repart. des *) (* valeurs de Khi[Seq] *) Stat := StatKS (FctnRepartKhiDeux, IToLR(t), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); DeleteTab (Khi); DeleteTab (NbExpected); DeleteTab (Count) END Gap; (****************************************************************************) PROCEDURE SimpPoker (N, n, d, k : INTEGER); CONST Word = 32; VAR j, r, Seq, (* Numero de la sequence courante. *) Groupe, U, Minkd, NbRegroupements : INTEGER; (* Nombre de "groupes" (classes) de *) (* valeurs de r. *) IndLigne, IndCol : INTEGER; (* Indice de la ligne et la colonne *) (* courants de MatStirling *) dReal, (* valeur de d convertie en reel *) Mult, KhiSeq, (* Variable tampon utilisee pour sau-*) (* ver du temps d'acces dans Khi[Seq]*) dk, (* Nombre de cases. *) Diff : LONGREAL; MatStirling, OccursBool : FlexTab; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' SimpPoker (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(d,1); WriteString(', '); WriteInt(k,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); BEGIN NewFlexTab (MatStirling, R); NewFlexTab (OccursBool, Bool); dReal := IToLR (d); IF (k < d) THEN Minkd := k ELSE Minkd := d END; dk := Pow (dReal, IToLR (k)) END ProcInit; BEGIN Entete ('Simplified Poker Test'); WriteString ('Avec NSeq = '); WriteInt (N, 4); WriteString (', d = '); WriteInt (d, 4); WriteLn; WriteString (', n = '); WriteInt (n, 4); WriteString (' et k = '); WriteInt (k, 4); WriteLn; WriteLn; IF (d >= 128) OR (k >= 128) THEN MessageErreur(30) END; ProcInit; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab ((Minkd + 1) * (k + 1), MatStirling); (* Dim. de MatStirling [0..Minkd, 0..k] *) CreateTab (Minkd + 1, NbExpected);(* Dim. de NbExpected [0...Minkd] *) CreateTab (Minkd + 1, Count); (* Domaine de Count [0..Minkd] *) CreateTab (Minkd + 1, Loca); (* Domaine de Loca [0..Minkd] *) CreateTab (d, OccursBool); NbRegroupements := Minkd; (* Calcul des probabilites et re- *) Loca^.Vi^[0] := 0; (* groupement des classes. *) Mult := IToLR (n) / dk; (* Necessaire pour le WHILE ci-bas. *) (* Calcul des elements de la matrice MatStirling. La valeur de chaque *) (* element [m,n] est calculee en utilisant les elements [m,n-1] et *) (* [m-1,n-1] de la colonne precedente. *) MatStirling^.Vr^[0 * (k+1) + 0] := OneLR; FOR IndLigne := 1 TO Minkd DO MatStirling^.Vr^[IndLigne * (k+1) + 0] := ZeroLR END; FOR IndCol := 1 TO k DO FOR IndLigne := 0 TO Minkd DO IF (IndLigne = 0) THEN MatStirling^.Vr^[IndLigne * (k+1) + IndCol] := ZeroLR ELSE MatStirling^.Vr^[IndLigne * (k+1) + IndCol] := (IToLR (IndLigne) * MatStirling^.Vr^[IndLigne * (k+1) + IndCol - 1]) + MatStirling^.Vr^[(IndLigne - 1) * (k+1) + IndCol - 1]; END END END; (* NbExpected[r] = n * d*(d-1)* ... *(d-r+1) * MatStirling[r,k] / d**k *) FOR r := 1 TO Minkd DO Mult := Mult * IToLR (d - r + 1); NbExpected^.Vr^[r] := (Mult * MatStirling^.Vr^[r * (k+1) + k]); Loca^.Vi^[r] := r END; IF Detailed THEN WriteString ('Les nombres esperes :'); WriteLn; WriteString ('Classe r Loca[r] NbExpected[r]'); WriteLn; FOR r := 1 TO Minkd DO WriteInt (r, 4); WriteInt (Loca^.Vi^[r], 13); WriteLR (NbExpected^.Vr^[r], 20, 4, 2); WriteLn END; WriteLn; END; FOR r := 1 TO Minkd-1 DO (* On regroupe differentes valeurs *) (* de r afin que le nb. espere dans *) (* chaque case soit superieur ou *) (* egal a MinExpected. *) IF NbExpected^.Vr^[r] < MinExpected THEN NbExpected^.Vr^[r+1] := NbExpected^.Vr^[r+1] + NbExpected^.Vr^[r]; NbExpected^.Vr^[r] := ZeroLR; (* On regroupe r (et ceux qui y *) Loca^.Vi^[r] := r + 1; (* etaient regroupes s'il y a lieu) *) j := r - 1; (* avec r + 1. *) DEC (NbRegroupements); WHILE (j > 0) AND (Loca^.Vi^[j] = r) DO Loca^.Vi^[j] := r + 1; j := j - 1 END END END; IF NbExpected^.Vr^[Minkd] < MinExpected THEN (* Cas particulier pour la derniere *) (* classe,Si NbEspere < MinExpexted *) j := Minkd - 1; WHILE (j > 0) AND (NbExpected^.Vr^[j] = ZeroLR) DO j := j-1 (* Recherche la classe precedente *) END; (* valide *) FOR r := j+1 TO Minkd DO (* Change les regroupements incorrects*) Loca^.Vi^[r] := j END; NbExpected^.Vr^[j] := NbExpected^.Vr^[j] + NbExpected^.Vr^[Minkd]; NbExpected^.Vr^[Minkd] := ZeroLR; NbRegroupements := NbRegroupements - 1 END; IF Detailed THEN WriteString ('Nombre de classes de valeurs de r : '); WriteInt (NbRegroupements, 4); WriteLn; WriteLn; WriteString ('Les regroupements :'); WriteLn; WriteLn; WriteString ('Classe r Loca[r] NbExpected[r]'); WriteLn; FOR r := 1 TO Minkd DO WriteInt (r, 4); WriteInt (Loca^.Vi^[r], 13); WriteLR (NbExpected^.Vr^[r], 17, 1, 2); WriteLn END; WriteLn; WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR r := 1 TO Minkd DO Count^.Vi^[r] := 0 END; FOR Groupe := 1 TO n DO (* Tirage d'une main de poker *) FOR j := 0 TO ( d-1 ) DO OccursBool^.Vbool^[j] := FALSE END; r := 0; FOR j := 1 TO k DO U := TRUNC (Uniforme01 () * dReal); IF NOT (OccursBool^.Vbool^[U]) THEN OccursBool^.Vbool^[U] := TRUE; INC (r) (* r = nb. de valeurs differentes *) END; (* dans cette main. *) END; INC (Count^.Vi^[r]) END; (* Regroupement des valeurs obtenues *) FOR r := 1 TO Minkd DO IF Loca^.Vi^[r] <> r THEN INC (Count^.Vi^[Loca^.Vi^[r]], Count^.Vi^[r]) END; END; (* Calcul de la statistique du *) (* Khi-Deux pour un sequence *) KhiSeq := ZeroLR; FOR r := 1 TO Minkd DO IF (Loca^.Vi^[r] = r) THEN Diff := IToLR (Count^.Vi^[r]) - NbExpected^.Vr^[r]; KhiSeq := KhiSeq + (Diff * Diff) / NbExpected^.Vr^[r]; END END; IF Detailed THEN WriteLR (KhiSeq, 8, 3, 2); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultats*) END; END; Khi^.Vr^[Seq] := KhiSeq; END; (* Test de K-S sur la repartition des*) (* valeurs de Khi[Seq] *) Stat := StatKS (FctnRepartKhiDeux, IToLR (NbRegroupements - 1), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (MatStirling); DeleteTab (NbExpected); DeleteTab (Count); DeleteTab (Loca); DeleteTab (OccursBool); DISPOSE (MatStirling); DISPOSE (OccursBool) END SimpPoker; (****************************************************************************) PROCEDURE SimpPokerOver (N, n, d, k : INTEGER); BEGIN END SimpPokerOver; (****************************************************************************) PROCEDURE CouponCollector (N, n, d, t : INTEGER); CONST Word = 32; VAR Exit : BOOLEAN; NbRegroupements, j, k, r, u, Seq, (* Numero de la sequence courante. *) Segm, (* Nombre courant de fois que les d *) (* valeurs ont ete generees dans une *) (* sequence. *) IndLigne, IndCol, (* Indice de la ligne et la colonne *) (* courants de MatStirling *) IndColMax : INTEGER; dReal, (* valeur de d convertie en reel *) Mult, KhiSeq, (* Variable tampon utilisee pour sau-*) (* ver du temps d'acces dans Khi[Seq]*) Moyder, Diff : LONGREAL; OccursBool : FlexTab; (* Durant une collection, les valeurs*) (* qui ont deja ete generees. *) MatStirling : FlexTab; (* Sert a calculer la fonction *) (* recursive Stirling(m,n). *) (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' CouponCollector (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(d,1); WriteString(', '); WriteInt(t,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE NRepet () : INTEGER; (* Compte le nombre de valeurs qui ont du etre generees pour que *) (* chacune des d valeurs possibles apparaissent au moins une fois. *) VAR r, k, u : INTEGER; BEGIN r := 0 ; FOR k:=1 TO d DO REPEAT INC (r); IF (r >= t) THEN RETURN t END; u := TRUNC (Uniforme01 () * dReal) UNTIL NOT (OccursBool^.Vbool^[u]); OccursBool^.Vbool^[u] := TRUE (* Ici, k= nb. de valeurs diffe *) (* rentes observees a date *) END; RETURN r END NRepet; BEGIN NewFlexTab (MatStirling, R); NewFlexTab (OccursBool, Bool); dReal := IToLR (d); Mult := IToLR (n); Entete ('Coupon Collector Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', d = '); WriteInt (d, 4); WriteLn; WriteString (', n = '); WriteInt (n, 4); WriteString (' et t = '); WriteInt (t, 4); WriteLn; WriteLn; IF (t >= 50) THEN MessageErreur (40) END; IF (d >= t) OR (d <= 1) THEN MessageErreur (41) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab ((d+1) * t, MatStirling); (* MatStirling [0..d, 0..(t-1)] *) CreateTab (t+1, NbExpected); (* Dimension de NdExpected [0..t] *) CreateTab (t+1, Count); (* Dimension de NbExpected [0..t] *) CreateTab (t+1, Loca); (* Dimension de NbExpected [0..t] *) CreateTab (d, OccursBool); (* Calcul du NbEspere de segments de chaque longueur. *) (* *) (* NbExpected[r] = n * d! * Stirling (d-1, r-1) / d**r *) (* pour r = d a t - 1 ; *) (* *) (* NbExpected[t] = n * (1 - d! * Stirling (d, t-1) / d**t-1 *) (* *) (* NbExpected[t] doit etre >= 10. *) (* *) FOR r := 1 TO (d-1) DO Mult := Mult * IToLR(r) / dReal END; NbExpected^.Vr^[d] := Mult; (* Calcul des elements de la matrice MatStirling. La valeur de chaque *) (* element [m,n] est calculee en utilisant les elements [m,n-1] et *) (* [m-1,n-1] de la colonne precedente. *) MatStirling^.Vr^[0] := OneLR; (* Stirling (d-1, d-1) = 1 *) FOR IndLigne := 1 TO d DO MatStirling^.Vr^[IndLigne * t + 0] := ZeroLR END; FOR IndLigne := 0 TO d DO FOR IndCol := 1 TO (t - 1 - d) DO IF (IndLigne = 0) THEN MatStirling^.Vr^[IndLigne * t + IndCol] := ZeroLR ELSE MatStirling^.Vr^[IndLigne * t + IndCol] := IToLR (IndLigne) * MatStirling^.Vr^[IndLigne * t + IndCol - 1] + MatStirling^.Vr^[(IndLigne - 1) * t + IndCol - 1] END END END; IndColMax := t - d; FOR IndLigne := 1 TO d DO FOR IndCol := (t - d) TO IndColMax DO MatStirling^.Vr^[IndLigne * t + IndCol] := IToLR (IndLigne) * MatStirling^.Vr^[IndLigne * t + IndCol - 1] + MatStirling^.Vr^[(IndLigne - 1) * t + IndCol - 1] END; IndColMax := IndColMax + 1 END; FOR r := (d+1) TO (t-1) DO Mult := Mult / dReal; NbExpected^.Vr^[r] := (Mult * MatStirling^.Vr^[(d-1) * t + r-1]) END; NbExpected^.Vr^[t] := IToLR(n) - Mult * MatStirling^.Vr^[d * t + t-1]; Moyder := ZeroLR; IF Detailed THEN WriteString ('Tableau des NbExpected[r], pour r = d a t :'); WriteLn; FOR r := d TO t DO WriteInt (r, 4); WriteLR (NbExpected^.Vr^[r], 10, 4, 2); WriteLn; Moyder := Moyder + (IToLR (r) * NbExpected^.Vr^[r]) END; Moyder := Moyder/IToLR(n); WriteString ('Esperance mathematique de r = '); WriteLR (Moyder, 5, 2, 2); WriteLn; END; NbRegroupements := t - d + 1; FOR r := (d - 1) TO t DO Loca^.Vi^[r] := r END; FOR r := d TO (t - 1) DO (* On regroupe differentes valeurs de*) (* r afin que le nb. espere de *) (* valeurs soit >= MinExpected *) IF (NbExpected^.Vr^[r] < MinExpected) THEN NbExpected^.Vr^[r+1] := NbExpected^.Vr^[r+1] + NbExpected^.Vr^[r]; NbExpected^.Vr^[r] := ZeroLR; Loca^.Vi^[r] := r + 1; j := r - 1; DEC (NbRegroupements); WHILE ((j > d - 1) AND (Loca^.Vi^[j] = r)) DO Loca^.Vi^[j] := r + 1; j := j - 1; END; END; END; (* Cas particulier pour la derniere *) (* classe, si NbEspere < MinExpected *) IF (NbExpected^.Vr^[t] < MinExpected) THEN j := t - 1; WHILE ((j > d-1) AND (NbExpected^.Vr^[j] = ZeroLR)) DO j := j-1 (* Recherche de la classe precedente *) END; (* valide. *) FOR r:= j+1 TO t DO (* Change les regroupements incorrects*) Loca^.Vi^[r] := j END; NbExpected^.Vr^[j] := NbExpected^.Vr^[j] + NbExpected^.Vr^[t]; NbExpected^.Vr^[t] := ZeroLR; NbRegroupements := NbRegroupements - 1 END; IF Detailed THEN WriteString ('Nombre de classes de valeurs de r : '); WriteInt (NbRegroupements, 4); WriteLn; WriteLn; WriteString ('Les regroupements :'); WriteLn; WriteString ('Classe r Loca[r] NbExpected[r]'); WriteLn; FOR r := d TO t DO WriteInt (r, 4); WriteInt (Loca^.Vi^[r], 13); WriteLR (NbExpected^.Vr^[r], 17, 1, 2); WriteLn; END; WriteLn; WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn END; (* Debut du test. *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR r := d TO t DO Count^.Vi^[r] := 0 END; FOR Segm := 1 TO n DO (* Debut d'une collection de valeurs.*) FOR k := 0 TO d-1 DO OccursBool^.Vbool^[k] := FALSE END; r := NRepet (); INC (Count^.Vi^[r]) END; (* Regroupements des valeurs obtenues*) FOR r := d TO t DO IF (Loca^.Vi^[r] <> r) THEN INC(Count^.Vi^[Loca^.Vi^[r]], Count^.Vi^[r]) END; END; (* Calcul de la statistique du *) (* Khi-Deux pour une sequence. *) KhiSeq := ZeroLR; FOR r := d TO t DO IF (Loca^.Vi^[r] = r) THEN Diff := IToLR (Count^.Vi^[r]) - NbExpected^.Vr^[r]; KhiSeq := KhiSeq + (Diff * Diff) / NbExpected^.Vr^[r] END END; IF Detailed THEN WriteLR (KhiSeq, 7, 3, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultats*) END; END; Khi^.Vr^[Seq] := KhiSeq; END; (* Test de K-S sur la repartition des*) (* valeurs de Khi[Seq]. *) Stat := StatKS (FctnRepartKhiDeux, IToLR (NbRegroupements - 1), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (MatStirling); DeleteTab (NbExpected); DeleteTab (Count); DeleteTab (Loca); DeleteTab (OccursBool); DISPOSE (MatStirling); DISPOSE (OccursBool) END CouponCollector; (****************************************************************************) PROCEDURE Permutation ( N : INTEGER; n : LONGINT; t : INTEGER (* t < 9. *) ); VAR f, (* 0 <= f <= t! Valeur unique *) (* attribuee selon l'ordre des t *) (* valeurs d'uniformes generees *) k, r, s, (* Indice du vecteur U pris entre 1 *) (* et t ou se retrouve la valeur la *) (* plus forte du vecteur *) Seq, (* Numero de la sequence courante *) Factt : INTEGER; Groupe : LONGINT; U : ARRAY [1..12] OF LONGREAL; (* Vecteur de t valeurs d'uniformes *) (* reparties entre 0 et 1 *) KhiSeq, (* Variable tampon utilisee pour sau-*) (* ver du temps d'acces dans Khi[Seq]*) NbExpected, Diff : LONGREAL; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Permutation (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); WriteInt(t,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) BEGIN Entete ('Permutation Test'); (* Imprime l'entete du rapport *) WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', n = '); WriteLongInt (n, 4); WriteString (' et t = '); WriteInt (t, 4); WriteLn; Factt := TRUNC (Factorielle (t)); NbExpected := LIToLR (n) / IToLR (Factt); WriteString ('Nb. de cases = t! = '); WriteInt (Factt, 7); WriteLn; WriteString ('Nb. espere par case = '); WriteLR (NbExpected, 7, 1, 2); WriteLn; WriteLn; IF (t > 8) THEN MessageErreur (51) END; (*SUN*) IF (NbExpected < MinExpected) THEN MessageErreur (50) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (Factt -1, Count); IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du Test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR k := 0 TO (Factt - 1) DO Count^.Vi^[k] := 0 END; FOR Groupe := 1 TO n DO FOR k := 1 TO t DO (* Generation d'un vecteur et analyse*) U[k] := Uniforme01 () (* de la permutation. *) END; f := 0; FOR r := t TO 2 BY -1 DO (* On trouve U[s], le max entre *) s := 1; (* U[1], U[2], ..., U[r]. *) FOR k := 2 TO r DO IF U[k] > U[s] THEN s := k END END; f := r * f + s - 1; U[s] := U[r] END; INC (Count^.Vi^[f]) END; (* Calcul du Khi-Deux. *) KhiSeq := ZeroLR; FOR r := 0 TO (Factt - 1) DO Diff := IToLR (Count^.Vi^[r]) - NbExpected; KhiSeq := KhiSeq + (Diff * Diff) END; KhiSeq := KhiSeq / NbExpected; IF Detailed THEN WriteLR (KhiSeq, 7, 3, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultats*) END; END; Khi^.Vr^[Seq] := KhiSeq; END; (* Test de K-S sur la repartition *) (* des valeurs de Khi[Seq]. *) Stat := StatKS (FctnRepartKhiDeux, IToLR (Factt - 1), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (Count); END Permutation; (****************************************************************************) PROCEDURE Run ( N : INTEGER; n : LONGINT; Up : BOOLEAN ); VAR nReal, (* convertion du parametre n en reel *) KhiSeq, Uprecedent, (* Valeur de U precedente *) U : LONGREAL; A : ARRAY[1..6], [1..6] OF LONGREAL; B : ARRAY[1..6] OF LONGREAL; (* Les constantes necessaires pour *) (* le calcul de R *) Longueur, (* Longueur courante de la sequence *) Seq, (* Numero de la sequence courante *) i, J : INTEGER; I : LONGINT; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Run (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); IF Up THEN WriteString( 'TRUE' ) ELSE WriteString( 'FALSE' ) END; WriteString (' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); VAR I : INTEGER; BEGIN nReal := LIToLR (n); A[1,1] := StringToLR ("4529.4"); A[1,2] := StringToLR ("9044.9"); A[1,3] := StringToLR ("13568.0"); A[1,4] := StringToLR ("18091.0"); A[1,5] := StringToLR ("22615.0"); A[1,6] := StringToLR ("27892.0"); A[2,2] := StringToLR ("18097.0"); A[2,3] := StringToLR ("27139.0"); A[2,4] := StringToLR ("36187.0"); A[2,5] := StringToLR ("45234.0"); A[2,6] := StringToLR ("55789.0"); A[3,3] := StringToLR ("40721.0"); A[3,4] := StringToLR ("54281.0"); A[3,5] := StringToLR ("67852.0"); A[3,6] := StringToLR ("83685.0"); A[4,4] := StringToLR ("72414.0"); A[4,5] := StringToLR ("90470.0"); A[4,6] := StringToLR ("111580.0"); A[5,5] := StringToLR ("113262.0"); A[5,6] := StringToLR ("139476.0"); A[6,6] := StringToLR ("172860.0"); FOR I := 2 TO 6 DO FOR J := 1 TO (I-1) DO A[I,J] := A[J,I]; END END; B[1] := OneLR/StringToLR ("6.0"); B[2] := StringToLR ("5.0") / StringToLR ("24.0"); B[3] := StringToLR ("11.0") / StringToLR ("120.0"); B[4] := StringToLR ("19.0") / StringToLR ("720.0"); B[5] := StringToLR ("29.0") / StringToLR ("5040.0"); B[6] := OneLR/StringToLR ("840.0"); END ProcInit; BEGIN ProcInit; Entete ('Run Test'); WriteString ('Avec NSeq = '); WriteInt (N, 4); WriteString (', n = '); WriteLongInt (n, 4); WriteString (' et Up = '); IF Up THEN WriteString ('TRUE'); WriteLn ELSE WriteString ('FALSE'); WriteLn END; WriteLn; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (6, Count); IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test. *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR i := 1 TO 6 DO Count^.Vi^[i] := 0 END; Longueur := 1 ; Uprecedent := Uniforme01 (); FOR I := 1 TO n DO U := Uniforme01 (); IF (Up AND (U < Uprecedent)) OR (NOT Up AND (U > Uprecedent)) THEN (* On termine un "Run". *) INC (Count^.Vi^[Longueur]); Longueur := 1 ELSE IF (Longueur < 6) THEN INC (Longueur) END END; Uprecedent := U END; INC (Count^.Vi^[Longueur]); (* Calcul de la statistique du *) (* KhiDeux pour un sequence *) KhiSeq := ZeroLR; FOR i := 1 TO 6 DO FOR J := 1 TO 6 DO KhiSeq := KhiSeq + A[i ,J] * ( IToLR(Count^.Vi^[i]) - nReal * B[i]) * ( IToLR(Count^.Vi^[J]) - nReal * B[J]) END END; Khi^.Vr^[Seq] := KhiSeq / nReal; IF Detailed THEN WriteLR (Khi^.Vr^[Seq], 7, 3, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn (* ligne dans le fichier de resultats*) END END; END; (* Test de K-S sur la repartition des*) (* valeurs de Khi[Seq]. *) Stat := StatKS (FctnRepartKhiDeux, IToLR (6), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (Count); END Run; (****************************************************************************) PROCEDURE MaxOft (N, n, d, t : INTEGER); VAR U, (* Variable contenant la valeur de *) (* l'uniforme courante *) Diff, (* Variable de travail *) MaxU, NbExpected, (* Nombre espere dans chaque classe *) KhiSeq, dReal, (* convertion du parametre d en reel *) tReal : LONGREAL; (* convertion du parametre t en reel *) j, Seq, (* Numero de la sequence courante *) Groupe, Indice : INTEGER; (* Indice du vecteur Count *) (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' MaxOft (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(d,1); WriteString(', '); WriteInt(t,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); BEGIN dReal := IToLR (d); tReal := IToLR (t); NbExpected := IToLR (n) / dReal END ProcInit; BEGIN ProcInit; (* Imprime l'entete du rapport *) Entete ('Maximum Of t Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', d = '); WriteInt (d, 3); WriteString (', n = '); WriteInt (n, 4); WriteString (' et t = '); WriteInt (t, 2); WriteLn; WriteLn; WriteString ('Nb. espere par classe = '); WriteLR (NbExpected, 7, 1, 2); WriteLn; WriteLn; IF (NbExpected < MinExpected) THEN MessageErreur (70) END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (d-1, Count); IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test. *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR Indice := 0 TO (d - 1) DO Count^.Vi^[Indice] := 0 END; FOR Groupe := 1 TO n DO (* Generation d'un vecteur et re- *) MaxU := Uniforme01(); (* cherche de la valeur max. *) FOR j := 1 TO (t - 1) DO U := Uniforme01 (); IF (U > MaxU) THEN MaxU := U END END; Indice := TRUNC (Pow(MaxU, tReal) * dReal); INC (Count^.Vi^[Indice]); END; (* Calcul du Khi-Deux pour une seq. *) KhiSeq := ZeroLR; FOR Indice := 0 TO (d -1) DO Diff := IToLR (Count^.Vi^[Indice]) - NbExpected; KhiSeq := KhiSeq + ( Diff * Diff ) END; KhiSeq := KhiSeq / NbExpected; Khi^.Vr^[Seq] := KhiSeq; IF Detailed THEN WriteLR (KhiSeq, 7, 2, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn END (* ligne dans le fichier de resultats*) END; END; (* Test de K-S sur la repartition *) (* des valeurs de F[Seq]. *) Stat := StatKS (FctnRepartKhiDeux, (dReal - OneLR), Khi, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (Count) END MaxOft; (***************************************************************************) (* Test de Collision. *) VAR NGen : INTEGER; (* Variable utilisee dans la fonction *) (* FctnRepartNbColl qui correspond *) (* au nombre de valeurs generees par *) (* sequence *) A : FlexTab; PROCEDURE FctnRepartNbColl ( Inutilise, (* Pour respecter le type FctnRepart. *) k : LONGREAL ) : LONGREAL; (* Calcule la probabilite, sous H0, d'avoir au plus k collisions. *) VAR kInt, j : INTEGER; R : LONGREAL; BEGIN R := ZeroLR; kInt := Entier (k); FOR j := (NGen - kInt) TO NGen DO R := R + A^.Vr^[j] END; RETURN (R) END FctnRepartNbColl; PROCEDURE Collision (N, n, dint, k : INTEGER); CONST Word = 16; (* Type BITSET contient 16 elements *) long = 15; (* ==> k <= 15 *) VAR m, (* Nombre total de cases dans les- *) (* quelles n vecteurs sont generes. *) d, (* Parametre d en longint *) Indice, (* Indice du tableau converti en *) (* indice du vecteur. *) Base, (* Compteur permettant la conversion *) (* d'un indice tableau en indice *) (* vecteur. *) Basek : LONGINT; (* Base a la k *) Suite : ARRAY[1..long] OF LONGINT; (* Contient la suite de nombre pour *) (* un test avec overlap = TRUE *) Seq, (* Numero de la sequence courante. *) pos, (* Position dans le tableau Suite *) h, i, J0, J1, l, Dim, (* Dimension maximale du tableau *) (* OccursBit. *) Ensemble : INTEGER; (* Ensemble de bits reference *) dReal, (* Convertion du parametre d en reel *) kReal, (* Convertion du parametre k en reel *) mReal, (* Convertion du parametre m en reel *) (* et double precision. *) UnSurm, X : LONGREAL; Position : CARDINAL; (* Pos. d'un bit dans un ensemble *) OccursBSet, NbCollisions : FlexTab; (* Contient pour chaque sequence le *) (* nombre de collisions survenues. *) PROCEDURE ProcInit (); BEGIN d := VAL (LONGINT, dint); dReal := LIToLR (d); kReal := IToLR (k); mReal := Pow (dReal, kReal); m := Trunc(mReal); UnSurm := OneLR / mReal; Dim := VAL (INTEGER, m DIV Word); NewFlexTab (NbCollisions, R); NewFlexTab (A, R); NewFlexTab (OccursBSet, BSet) END ProcInit; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Collision (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(dint,1); WriteString(', '); WriteInt(k,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) BEGIN ProcInit; Entete ('Collision Test'); WriteString ('Avec N = '); WriteInt (N, 4); WriteString (', n = '); WriteInt (n, 4); WriteString (', d : '); WriteInt (dint, 4); WriteString (' et k : '); WriteInt (k, 2); WriteLn; IF n > 131072 (*Si > 2^17*) THEN MessageErreur (81) END; (*SUN*) IF (mReal > StringToLR("2097152.0")) THEN (* Si > (2^17)*Word *) (*SUN*) MessageErreur (80); END; (* Creation des tableaux necessaires pendant le test *) CreateTab (N+1, Khi); CreateTab (N, NbCollisions); CreateTab (n, A); CreateTab (Dim, OccursBSet); IF Detailed THEN WriteString ('Nombre de collisions, pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test. *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR l := 0 TO Dim DO OccursBSet^.VBSet^[l] := {} END; NbCollisions^.Vr^[Seq] := ZeroLR; FOR l := 1 TO n DO Indice := Trunc (Uniforme01 () * dReal); Base := 1; FOR h := 2 TO k DO Base := Base * d; INC (Indice, Trunc (Uniforme01 () * dReal) * Base) END; Ensemble := VAL (INTEGER , Indice DIV Word); Position := VAL (CARDINAL, Indice MOD Word); IF (Position IN OccursBSet^.VBSet^[Ensemble]) THEN NbCollisions^.Vr^[Seq] := NbCollisions^.Vr^[Seq] + OneLR ELSE INCL (OccursBSet^.VBSet^[Ensemble], Position) END; END; IF Detailed THEN WriteLR (NbCollisions^.Vr^[Seq], 7, 0, 7); IF (Seq MOD 5 = 0) THEN WriteLn END END END; (* Traitement des resultats. *) FOR l := 0 TO n DO A^.Vr^[l] := ZeroLR END; A^.Vr^[1] := OneLR; J0 := 1; J1 := 1; FOR l := 1 TO n - 1 DO INC (J1); i := J1; WHILE (i >= J0) DO X := IToLR (i) * UnSurm; A^.Vr^[i] := (X * A^.Vr^[i]) + ((OneLR + UnSurm - X) * A^.Vr^[i-1]); IF (A^.Vr^[i] < DixExpMoinsVingt) THEN IF (i = J1) THEN DEC (J1) ELSE IF (i = J0) THEN INC (J0) END END END; DEC (i); END; END; NGen := n; (* Test de K-S sur le nombre de *) (* collisions. *) Stat := StatKS (FctnRepartNbColl, StringToLR ("99.9"), NbCollisions, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Khi); DeleteTab (NbCollisions); DeleteTab (A); DeleteTab (OccursBSet); DISPOSE (NbCollisions); DISPOSE (A); DISPOSE (OccursBSet) END Collision; (******************************************************************************) PROCEDURE CollisionOver (N, n, dint, k : INTEGER); BEGIN END CollisionOver; (******************************************************************************) PROCEDURE Opso ( N, (* nombre de repetitions du test; *) r, (* r premiers bits de U(0,1) rejetes *) p : INTEGER ); (* ensemble de parametres choisis *) (* pour le test: p = 1 pour le PC. *) (* p = 1, 2, 3 sur SUN *) CONST SetSize = 8 * SIZE(BITSET); (* Nombre de bits dans un BITSET *) (*SUN*) nCasesMax = 4194304; (* Nombre maximum de cases = 2^22 *) (*SUN*) MaxIndex = 131072; (* nCasesMax DIV SetSize *) TYPE BsetVector = POINTER TO ARRAY [0..MaxIndex] OF BITSET; VAR rtBits : LONGREAL; (* r+t bits apres le point de U(0,1) *) t, i, Seq : INTEGER; (* Numero de la sequence courante. *) PosDansMot, (* Le bit choisi dans l'ensemble *) Mot, (* L'ensemble choisi *) nBalles, (* Nombre de balles lancees *) tmask, (* Masque pour les t bits choisis *) Balle, (* Numero de la balle courante *) X, Y, (* Paire de nombres aleatoires *) CasesVides, (* Nombre de cases vides *) Case : CARDINAL; (* La case choisie par la paire X, Y *) CasesNormee : FlexTab; OccursBSet : BsetVector; NBalles, NCases : ARRAY [1..3] OF CARDINAL; Esperance, EcartType : ARRAY [1..3] OF LONGREAL; (*****************) PROCEDURE ProcInit (); BEGIN NCases[1] := 1048576; (* Nombre de cases = 2^20 *) NCases[2] := nCasesMax; (* " = 2^22 *) NCases[3] := nCasesMax; NBalles[1] := 2097152; (* Nombre de balles lancees = 2^21 *) NBalles[2] := 4194304; (* " = 2^22 *) NBalles[3] := 8388608; (* " = 2^23 *) Esperance[1] := 141909.0; (* Esperance theorique du nombre de *) Esperance[2] := 1542998.0; (* cellules vides. *) Esperance[3] := 567637.0; EcartType[1] := 290.26; (* Ecart-type theorique du nombre de *) EcartType[2] := 638.75; (* cellules vides. *) EcartType[3] := 580.80; IF p = 1 THEN t := 10 ELSE t := 11 END; rtBits := TwoExp [r + t]; tmask := TRUNC ( TwoExp [t] ); NewFlexTab (CasesNormee, R); nBalles := NBalles[p]; END ProcInit; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Opso (' ); WriteInt(N,1); WriteString(', '); WriteInt(r,1); WriteString(', '); WriteInt(p,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) BEGIN Entete ('Opso test'); WriteString ('N (Sequences) = '); WriteInt (N, 1); WriteString (' ;r = '); WriteInt (r, 1); WriteString (' ;p = '); WriteInt (p, 1); WriteLn (); IF NOT ((p = 1) OR (p=2) OR (p=3)) THEN MessageErreur (121) END; (*SUN*) IF ( r < 0 ) OR ( r > 20 ) THEN MessageErreur (120) END; ProcInit (); WriteString ('nCases = '); WriteCard (NCases[p], 1); WriteLn (); WriteString ('nBalles = '); WriteCard (NBalles[p], 1); WriteLn (); WriteString ('Moyenne theorique = '); WriteLR (Esperance[p], 10, 3, 8); WriteLn (); WriteString ('Ecart-type theorique = '); WriteLR (EcartType[p], 10, 3, 8); WriteLn (); (* Creation des tableaux necessaires pendant le test. *) CreateTab (N, CasesNormee); ALLOCATE ( OccursBSet, TSIZE (BITSET)*(MaxIndex + 1) ); IF Detailed THEN WriteString('Nombre de cases vides, pour chaque sequence :'); WriteLn; WriteLn END; (* Debut du test. *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR i := 0 TO MaxIndex DO OccursBSet^[i] := {} END; CasesVides := NCases[p]; (* au debut, toutes les cases vides *) X := TRUNC (rtBits * Uniforme01 ()); (* on garde les (r + t) premiers bits apres le point *) X := X MOD tmask; (* ne garde que les t bits choisis *) FOR Balle := 1 TO nBalles DO Y := TRUNC (rtBits * Uniforme01 ()); Y := Y MOD tmask; Case := (X * tmask) + Y; PosDansMot := Case MOD SetSize; Mot := Case DIV SetSize; IF NOT (PosDansMot IN OccursBSet^[Mot]) THEN INCL (OccursBSet^[Mot], PosDansMot); DEC ( CasesVides ); END; X := Y; END; CasesNormee^.Vr^[Seq] := ( IToLR(CasesVides) - Esperance[p] ) / EcartType[p]; IF Detailed THEN WriteCard (CasesVides, 1); WriteString (" "); IF (Seq MOD 5 = 0) THEN WriteLn END; END; END; (* Debut du traitement des resultats *) Stat := StatKS (FctnRepartNormale, 99.9, CasesNormee, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); (* Destruction des tableaux utilises pendant le test. *) DeleteTab (CasesNormee); DISPOSE (CasesNormee); DEALLOCATE ( OccursBSet, TSIZE (BITSET)*(MaxIndex + 1) ); END Opso; (******************************************************************************) PROCEDURE PeriodsInStrings ( N : INTEGER; n : LONGINT; r, s : CARDINAL ); (* Dans cette implantation, les bits apres le point sont numerotes 31, 30, *) (* 29, etc.. Ainsi, le bit 0 est le 32-ieme bit apres le point pour U(0,1) *) (* , et est aussi l'element 0 de l'ensemble BinSet *) (* Etant donne que la procedure Psi retourne un LONGINT, cela impose la limitation suivante: s <= 30. Si la restriction n'est pas respectee, il y a une erreur a la ligne IF (k <= 0) THEN RETURN VAL (LONGINT, TwoExp [-k]) END; de la procedure Psi quand k <= -31. *) CONST max = 31; (* Nbre elements dans type BinSet - 1 *) TYPE BinSet = SET OF [0..max]; ListC = POINTER TO InfoListC; (* Une liste de correlations. *) InfoListC = RECORD Nb : LONGINT; (* Durant les precalculs: nb de *) (* bits de la correlation. *) (* Ensuite: nb d'occurences. *) C : BinSet; (* Une correlation. *) Pop : LONGINT; (* Population associee a C (et(c)). *) Ext : ListC; (* Pointeur sur la plus petite *) (* extension de C (NIL si aucune). *) Ext0 : ListC; (* Si C est une extension de D, *) (* pointe sur la plus petite exten- *) (* sion de D plus longue que C. *) (* = NIL si aucune a date. *) Next : ListC (* Correlation suiv. dans la liste *) (* des correl. de cette taille. *) END; VAR Corr : ARRAY [0..32] OF ListC; (* Listes des correlations pour *) (* chaque taille de chaine. *) Ci, OldCi : ListC; Seq, j : INTEGER; k : CARDINAL; NbReste, i, NbCases : LONGINT; KhiSeq, Exp_, (* Nombre d'occurences attendues d'une *) (* correlation donnee. *) ExpReste, rsBits, (* rsBits := 2^(s+r) *) (* Les chaines de s bits sont vues com-*) (* me des chiffres de 0 a (2^s) - 1. *) Fraction : LONGREAL; MaskS, (* Masque pour extraire les s bits *) maskA : LONGCARD; D : BinSet; (*****************) PROCEDURE WriteCorr (); BEGIN (* Imprimer entete, et imprimer les corr. et leur population. *) WriteString ('Liste des correlations de taille s,'); WriteLn; WriteString ('et leur population : '); WriteLn; WriteLn; Ci := Corr[s]; WHILE (Ci # NIL) DO (* Imprimer Ci^.C et Ci^.Pop *) FOR k := 0 TO s-1 DO IF (k IN Ci^.C) THEN WriteCard (1, 1) ELSE WriteCard (0, 1) END END; WriteLongInt (Ci^.Pop, 10); WriteLn; Ci := Ci^.Next END END WriteCorr; (*****************) PROCEDURE Psi ( C : BinSet; (* Correlation qui est la partie ter- *) (* minale de la correlation k. *) j, (* Longueur de la correlation C *) k : INTEGER (* Correlation consistant d'un 1 suivie*) (* de 0's jusqu'au suffix C, qui est *) ) : LONGINT; (* correlation. k = 100...000C *) BEGIN (* j<=>c & k<=>k de l'article *) IF (k > j) THEN RETURN 0 END; IF (k <= 0) THEN RETURN Trunc (TwoExp [-k]) END; IF (VAL (CARDINAL, (j - k)) IN C) THEN RETURN 1 ELSE RETURN 0 END END Psi; (*****************) PROCEDURE ComputeCorr (); VAR OldCi, Ci, Cj, CjE : ListC; i, j : CARDINAL; (* ne peut pas etre CARDINAL: Peut etre < 0 *) Tmax, p : LONGINT; k : CARDINAL; BEGIN NEW (Corr [0]); WITH Corr [0]^ DO Nb := 0; Pop := 1; Ext := NIL; Ext0 := NIL; Next := NIL END; NEW (Corr [1]); WITH Corr [1]^ DO Nb := 1; INCL (C, 0); Pop := 2; Ext := NIL; Ext0 := NIL; Next := NIL END; FOR i := 2 TO s DO (* i est la taille des chaines. *) (* Compter et former la liste des correlations de taille i.*) NEW (Ci); Corr [i] := Ci; FOR j := 0 TO i-2 DO (* j est la taille des cor. C dans *) (* 100...00C *) Cj := Corr [j]; WHILE (Cj # NIL) DO WITH Cj^ DO (* Calculer le nb. de chaines de long. i et de correl. *) (* "10...0C", et si > 0, ajouter cette corr. a Corr[i].*) p := Pop * Psi (C, j, INTEGER (2 * j - i)); (* Verifier si 1C peut etre une correlation. Vrai seule-*) (* ment si C = 111....1, i.e. derniere de la liste: *) IF (Next = NIL) THEN p := p - 2 * Psi (C, j, INTEGER (2 * j + 2 - i)) END; CjE := Ext; (* Note: j <= j-2, i.e. j+1 <= (i+j) DIV 2 *) Tmax := IToLI ((i+j) DIV 2 ); WHILE ((CjE # NIL) AND (CjE^.Nb <= Tmax)) DO p := p - CjE^.Pop * Psi (C, j, INTEGER (2 * CjE^.Nb - IToLI(i))); CjE := CjE^.Ext0 END; (* p = nb. de chaines voulu. *) IF (p > 0) THEN (* On met cette correlation dans la Ci. *) WITH Ci^ DO Nb := IToLI (i); Pop := p; Ext := NIL; Ext0 := NIL; (* Ci^.C devient Cj^.C decale a droite de i-j *) (* positions, avec un 1 en premiere position. *) C := BinSet{0}; IF (j > 0) THEN FOR k := 0 TO j-1 DO IF (k IN Cj^.C) THEN INCL (C, k + VAL (CARDINAL, i - j)) END END END; END; IF (Ext = NIL) THEN Ext := Ci ELSE CjE := Ext; WHILE (CjE^.Ext0 # NIL) DO CjE := CjE^.Ext0 END; CjE^.Ext0 := Ci END; OldCi := Ci; NEW (Ci); OldCi^.Next := Ci END END; Cj := Cj^.Next END; END; (* Pour j = i-1, on a la correlation "11...1". *) WITH Ci^ DO C := BinSet{}; FOR k := 0 TO i - 1 DO INCL (C, k) END; Nb := IToLI (i); Pop := 2; Ext := NIL; Ext0 := NIL; Next := NIL END END; (* for i *) (* On recupere l'espace qui ne sera plus utilise. *) FOR i := 0 TO s - 1 DO Ci := Corr[i]; WHILE (Ci # NIL) DO OldCi := Ci; Ci := Ci^.Next; DISPOSE (OldCi); END END; IF Detailed THEN WriteCorr END; END ComputeCorr; (*****************) PROCEDURE InitNb (); (* Initialiser les compteurs pour les correlations de longueur s. *) BEGIN Ci := Corr[s]; WHILE (Ci # NIL) DO Ci^.Nb := 0; Ci := Ci^.Next END END InitNb; (*****************) PROCEDURE GenerateC () : BinSet; (* Genere un chaine de bits et en retourne la correlation. *) VAR d ,g : LONGCARD; c : BinSet; (* La determination de la correlation d'une chaine de bits consiste en la comparaison des chaines g et d (initialement egales). A chaque iteration k: 1) elimination du bit le plus a gauche de la chaine g 2) elimination du bit le plus a droite de la chaine d (decalage droit) 3) si g et d sont egales, le bit k de la correlation c est 1, sinon 0. *) BEGIN (* Generation d'un nombre aleatoire; selection des (r+s) bits de gauche *) d := TRUNC (rsBits * Uniforme01 ()); (* Extraction des s bits choisis *) d := d MOD MaskS; g := d; (* Initialisation de la correlation, cas trivial *) c := BinSet{0}; maskA := MaskS DIV 2; FOR k := 1 TO s-1 DO (* elimination du bit le plus a gauche de la chaine g *) g := g MOD maskA; (* elimination du bit le plus a droite de la chaine d *) d := d DIV 2; (* si g et d sont egales, le bit k de la correlation est 1, sinon 0 *) IF g = d THEN INCL(c, k) END; maskA := maskA DIV 2 END; RETURN c END GenerateC; (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' PeriodsInStrings (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); WriteCard(r,1); WriteString(', '); WriteCard(s,1); WriteString(')'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); VAR i : CARDINAL; BEGIN MaskS := 1; FOR i := 1 TO s DO MaskS := MaskS * 2 END; rsBits := TwoExp[s + r]; END ProcInit; (*****************) BEGIN (* Imprimer entete *) Entete ('Test PeriodsInStrings'); WriteString (' N = '); WriteInt (N, 4); WriteString (', n = '); WriteLongInt (n, 2); WriteString (', r = '); WriteCard (r, 4); WriteString (', s = '); WriteCard (s, 4); WriteLn; WriteLn; IF (r + s > 31) THEN MessageErreur (90) END; (*SUN*) IF (s > 30) OR (s < 2) THEN MessageErreur (91) END; ProcInit; (* Ne pas recalculer les correlations si c'est deja fait: *) IF (INTEGER (s) # CorrCalculees) THEN IF (CorrCalculees # 0) THEN (* Liberer l'espace utilise par la liste precedente: *) Ci := pCorrCalculees; WHILE (Ci # NIL) DO OldCi := Ci; Ci := Ci^.Next; DISPOSE (OldCi) END END; ComputeCorr; CorrCalculees := s; pCorrCalculees := Corr[s] ELSE Corr[s] := pCorrCalculees END; IF Detailed THEN WriteLn; WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; Fraction := LIToLR (n) / TwoExp[s]; (* Creation du tableau necessaire pour le test *) CreateTab (N+1, Khi); FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur () END; InitNb; FOR i := 1 TO n DO D := GenerateC (); Ci := Corr[s]; WHILE (Ci^.C # D) DO Ci := Ci^.Next END; INC (Ci^.Nb) END; NbCases := 0; NbReste := 0; ExpReste := ZeroLR; KhiSeq := ZeroLR; Ci := Corr[s]; WHILE (Ci # NIL) DO WITH Ci^ DO Exp_ := LIToLR (Pop) * Fraction; IF (Exp_ < MinExpected) THEN NbReste := NbReste + Nb; ExpReste := ExpReste + Exp_; IF (ExpReste >= MinExpected) THEN KhiSeq := KhiSeq + (LIToLR (NbReste) - ExpReste) * (LIToLR (NbReste) - ExpReste) / ExpReste; INC (NbCases); ExpReste := ZeroLR; NbReste := 0 END; ELSE KhiSeq := KhiSeq + (LIToLR (Nb) - Exp_) * (LIToLR (Nb) - Exp_) / Exp_; INC (NbCases) END; END; Ci := Ci^.Next END; IF Detailed THEN WriteLR (KhiSeq, 7, 3, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Khi par *) WriteLn END; (* ligne dans le fichier resultats *) END; Khi^.Vr^[Seq] := KhiSeq END; Stat := StatKS (FctnRepartKhiDeux, LIToLR (NbCases - 1), Khi, N); (* Test de K-S sur la repartition des*) (* valeurs de Khi[Seq] *) (* Imprimer resultats. *) WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); DeleteTab (Khi); END PeriodsInStrings; (****************************************************************************) (* Test NearestPair. *) VAR mu : LONGREAL; PROCEDURE FctnRepartXExpT (bidon, Valeur : LONGREAL) : LONGREAL; (* Retourne la valeur de la fonction de repartition a x = Valeur *) BEGIN RETURN OneLR - Exp (-Valeur/mu) END FctnRepartXExpT; PROCEDURE NearestPair (N, n, k : INTEGER); (* Test base sur un article de B. D. Ripley, Proc. R. Soc. Lond. A 389, *) (* 197-204(1983), surtout la proposition 2, page 199. *) (* Implantation: Preparata et Shamos: Computationnal Geometry *) (* Springer-Verlag, 1985, pp. 189-198, QA447P745 *) CONST pdim = 12; (* nombre maximum de recursions sur les dimensions: 12 pour (*SUN*) *) dimensionMax = 12; (* nombre maximum de dimensions *) nombreMaxPoints = 100000; (*SUN*) (* pour DOS; le nombre maximum de *) (* points a ete fixe empiriquement a *) (* 12000 DIV 2^{plancher(logBase2 k)} *) (* voir instruction CASE ci-apres *) TYPE (* vecteur^ contient un point dans l'hypercube a k dimensions. *) vecteur = POINTER TO ARRAY [1..dimensionMax] OF LONGREAL; (* pointeur^ contient les pointeurs vers tous les points *) pointeur = POINTER TO ARRAY [1..nombreMaxPoints] OF vecteur; VAR i, j, Seq, PCmax, p : INTEGER; (* la dimension courante *) dmin, (* la distance carree minimale *) PI, (* la constante pi *) Vol : LONGREAL; (* Volume de sphere unite en k dim. *) Points : ARRAY [1..dimensionMax] OF pointeur; (* p tableaux de pointeurs vers les points, un tableau pour chaque *) (* recursion sur les dimensions *) (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' NearestPair (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteInt(k,1); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE volume (); (* Calcul du volume d'une sphere unitaire en k dimensions: (base sur un *) (* developpement en serie d'une integrale; voir notes manuscrites) *) VAR i, j : INTEGER; BEGIN Vol := TwoLR; IF (k > 1) THEN FOR i := 2 TO k DO IF ODD (i) THEN FOR j := 2 TO i-1 BY 2 DO Vol := Vol * IToLR (j) / IToLR (j+1) END; Vol := Vol * TwoLR ELSE (* k est paire > 0 *) FOR j := 1 TO i-1 BY 2 DO Vol := Vol * IToLR (j) / IToLR (j+1) END; Vol := Vol * PI END END END; END volume; (*****************) PROCEDURE QuickSort ( A : pointeur; m, n, p : INTEGER ); (* On trie les points m a n en utilisant la p-ieme coordonnee comme *) (* cle; on permute les pointeurs au lieu de permuter les points. *) PROCEDURE Sort (l, r : INTEGER); VAR i, j : INTEGER; pivot : LONGREAL; vec : vecteur; BEGIN (* Sort *) i := l; j := r; pivot := A^[(l+r) DIV 2]^[p]; REPEAT WHILE (A^[i]^[p] < pivot) DO i := i + 1 END; WHILE (pivot < A^[j]^[p]) DO j := j - 1 END; IF (i <= j) THEN vec := A^[i]; A^[i] := A^[j]; A^[j] := vec; i := i + 1; j := j - 1 END; UNTIL (i > j); IF (l < j) THEN Sort (l, j) END; IF (i < r) THEN Sort (i, r) END END Sort; BEGIN (* QuickSort *) Sort (m, n) END QuickSort; (*****************) PROCEDURE distance( vecA, vecB: vecteur; VAR plusCourte: BOOLEAN ); (* Calcule le carre de la distance entre les points vecA et vecB. *) (* S'il y a lieu, met a jour la valeur de dmin, dans lequel cas, *) (* plusCourte devient TRUE *) VAR i : INTEGER; dist, temp : LONGREAL; BEGIN dist := ZeroLR; i := 1; WHILE (i <= k) AND (dist < dmin) DO temp := vecA^[i] - vecB^[i]; dist := dist + (temp * temp); INC(i); END; IF (i > k) AND (dist < dmin) THEN dmin := dist; plusCourte := TRUE END; END distance; (*****************) PROCEDURE binSearch ( x: LONGREAL; u, v, p: INTEGER) :INTEGER; (* trouve l'indice j dans [u..v] tel que A^[j]^[p] <= x < A^[j+1]^[p] *) (* les points sont ordonnes selon la p-ieme coordonnee. *) VAR j : INTEGER; A : pointeur; BEGIN A := Points[p]; IF A^[u]^[p] >= x THEN RETURN u - 1 END; IF A^[v]^[p] <= x THEN RETURN v END; REPEAT j := (u + v) DIV 2; IF (A^[j]^[p] > x) THEN v := j ELSE u := j END; UNTIL (A^[j]^[p] <= x) AND (A^[j+1]^[p] > x); RETURN j; END binSearch; (*****************) PROCEDURE dminSlice( VAR m: INTEGER; imed, jmed: INTEGER; VAR n: INTEGER; p : INTEGER ); (* E1 = {points indices [m..imed]} , E2 = {...[jmed..n]} *) (* 2 ensembles de points ordonnes selon la p-ieme coordonnee. Reduit *) (* E1 a une tranche de largeur racdmin (dans la direction p) du point *) (* le plus a gauche de E2; et de meme pour E2 par rapport au point le *) (* plus a droite de E1 (toujours par rapport a la p-ieme coordonnee. *) VAR i : INTEGER; temp, racdmin : LONGREAL; A : pointeur; BEGIN IF (m > imed) OR (jmed > n) THEN RETURN END; A := Points[p]; racdmin := Sqrt (dmin); temp := A^[jmed]^[p]; i := imed; WHILE (i >= m) AND (temp - A^[i]^[p] < racdmin) DO DEC(i) END; INC(i); m := i; temp := A^[imed]^[p]; i := jmed; WHILE (i <= n) AND (A^[i]^[p] - temp < racdmin) DO INC(i) END; DEC(i); n := i; END dminSlice; (*****************) PROCEDURE miniproc ( r, s, u, v, p : INTEGER ); (* calcule la distance entre les points de l'ensemble E1 = {indices *) (* [r..s]}, et ceux de E2 = {indices [u..v]}. On considere seulement *) (* les points dont les p-ieme coordonnees <= racdmin. *) VAR plusCourte : BOOLEAN; i, j, inf : INTEGER; racdmin, temp : LONGREAL; ivec : vecteur; A : pointeur; BEGIN plusCourte := TRUE; A := Points[p]; inf := u; FOR i := r TO s DO IF plusCourte THEN racdmin := Sqrt (dmin); plusCourte := FALSE END; ivec := A^[i]; temp := ivec^[p] - racdmin; (* considere seulement les points a *) (* distance racdmin selon la p-ieme *) (* coordonnee. *) j := inf; WHILE (j <= v) AND (A^[j]^[p] <= temp) DO INC(j) END; inf := j; (* indice limite inferieure courante *) temp := temp + TwoLR * racdmin; (* limite superieure courante *) WHILE (j <= v) AND (A^[j]^[p] < temp) DO distance ( ivec, A^[j] , plusCourte); INC(j) END; END; END miniproc; (*****************) PROCEDURE TransClosest ( r, s, u, v, p :INTEGER ); (* calcule la distance minimale entre les points de l'ensemble E1 et *) (* ceux de E2; E1 = {points indices [r..s]}, E2 = {...indices [u..v]} *) VAR i, j, imed : INTEGER; mediane : LONGREAL; A : pointeur; BEGIN IF (r > s) OR (u > v) THEN RETURN END; A := Points[p]; IF p = 1 THEN (* derniere recursion p/r dimensions *) miniproc( r, s, u, v, p ); RETURN END; IF (s - r < 24) OR (v - u < 24) THEN (* petits ensembles de points *) miniproc( r, s, u, v, p ); RETURN END; (* recursion p/r a la mediane de E1 selon la p-ieme coordonnee *) imed := ( r + s ) DIV 2; mediane := A^[imed]^[p]; (* tous les elements de E2 dont la p-ieme coordonnee <= mediane de E1 *) (* ont indice <= j *) j := binSearch( mediane, u, v, p ); TransClosest ( r, imed, u, j, p); (* recursion p/r au nombre de points*) TransClosest ( imed + 1, s, j + 1, v, p ); dminSlice ( r, imed, j + 1, v, p); (* candidats dans tranche racdmin *) dminSlice ( u, j, imed + 1, s, p); (* recursion p/r aux dimensions; introduire le tableau de pointeurs *) (* suivant et ordonner selon la coordonnee suivante *) FOR i := r TO s DO Points[p-1]^[i] := A^[i] END; FOR i := u TO v DO Points[p-1]^[i] := A^[i] END; IF (r <= imed) AND (j + 1 <= v) THEN QuickSort( Points[p-1], r, imed, p-1); QuickSort( Points[p-1], j+1, v, p-1); TransClosest( r, imed, j + 1, v, p-1); END; IF (u <= j) AND (imed + 1 <= s) THEN QuickSort( Points[p-1], imed+1, s, p-1); QuickSort( Points[p-1], u, j, p-1); TransClosest( u, j, imed + 1, s, p-1); END; END TransClosest; (*****************) PROCEDURE ClosestPair( m, n: INTEGER ); (* calcule la distance minimale entre toutes les paires de points *) (* indices [m..n]. *) VAR imediane, i : INTEGER; bidon : BOOLEAN; (* variable inutilisee ici *) A : pointeur; BEGIN A := Points[p]; IF n = m THEN RETURN END; (* ensemble de 1 point *) IF n = m+1 THEN (* ensemble de 2 points *) distance( A^[m], A^[n], bidon ); RETURN END; imediane := (m+n) DIV 2; ClosestPair( m, imediane ); (* recursion avec ensemble E1 *) ClosestPair( imediane+1, n ); (* recursion avec E2 *) IF k = 1 THEN (* 1 dimension: relation ordre total *) distance( A^[imediane], A^[imediane+1], bidon ); RETURN END; (* les seuls candidats restants -> a distance^2 dmin de la mediane *) dminSlice ( m, imediane, imediane + 1, n, p); IF (m > imediane) OR (imediane + 1 > n) THEN RETURN END; (* utiliser le tableau suivant; reordonner selon la coordonnee (p-1) *) FOR i := m TO n DO Points[p - 1]^[i] := A^[i] END; A := Points[p - 1]; QuickSort( A, m, imediane, p - 1); QuickSort( A, imediane+1, n, p - 1); (* calculer la distance minimale entre les points de E1 et ceux de E2 *) TransClosest (m, imediane, imediane + 1, n, p - 1); END ClosestPair; (*****************) BEGIN Entete ('NearestPair'); WriteString (' N = '); WriteInt (N, 4); WriteString ('; n = '); WriteInt (n, 4); WriteString ('; k = '); WriteInt (k, 4); WriteLn; WriteLn; IF (k > dimensionMax) THEN MessageErreur (100) END; (* limite empirique approximative du nombre de points sur un PC *) (* CASE k OF 1 : PCmax := nombreMaxPoints | 2, 3 : PCmax := nombreMaxPoints DIV 2 | 4..7 : PCmax := nombreMaxPoints DIV 4 | 8..12: PCmax := nombreMaxPoints DIV 8 END;*) IF (n > nombreMaxPoints) THEN MessageErreur (101) END; (*SUN*) IF k > pdim THEN p := pdim ELSE p := k END; (* Allocation des p tableaux de pointeurs vers les points, un pour chaque dimension jusqu'a p. *) FOR j := 1 TO p DO ALLOCATE (Points[j], ORD(n) * SIZE (vecteur)); END; (* Allocation de memoire pour les points; initialement, seul le p-ieme tableau de pointeurs pointe vers les points *) FOR i := 1 TO n DO ALLOCATE (Points[p]^[i], ORD(k) * SIZE (LONGREAL)) END; (* Creation du tableau necessaire pendant le test *) CreateTab (N+1, Loca); IF Detailed THEN WriteString ('Distance minimum, exposant '); WriteInt (k, 2); WriteString (', pour chaque sequence :'); WriteLn; WriteLn END; PI := StringToLR ("3.141592653589793"); volume; FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; IF k > pdim THEN p := pdim ELSE p := k END; FOR i := 1 TO n DO (* Generer n points de dimension k *) FOR j := 1 TO k DO Points[p]^[i]^[j] := Uniforme01 (); END END; (* En k dimensions, la plus grande valeur possible de dmin est k: *) (* viz. Uniforme01 *) dmin := IToLR (k); QuickSort( Points[p], 1, n, p ); ClosestPair( 1, n ); (* Elever a k/2 equivaut a extraire la racine et elever a la puiss. k *) dmin := Pow (dmin, IToLR (k) / TwoLR); IF Detailed THEN WriteLR (dmin, 7, 3, 3); WriteString (' '); IF (Seq MOD 10 = 0) THEN WriteLn END END; Loca^.Vr^[Seq] := dmin; END; (* FOR N *) IF k > pdim THEN p := pdim ELSE p := k END; FOR i := 1 TO n DO DEALLOCATE (Points[p]^[i], ORD(k) * SIZE (LONGREAL)) END; FOR j := 1 TO p DO DEALLOCATE (Points[j], ORD(n) * SIZE (vecteur)); END; (* mu est la moyenne de la loi exponentielle que suit dmin^k/2 *) mu := TwoLR / (IToLR (n) * IToLR (n) * Vol); Stat := StatKS (FctnRepartXExpT, IToLR (k), Loca, N); WriteResults (Stat, OneLR - FctnRepartKS (IToLR (N), Stat)); DeleteTab (Loca); END NearestPair; (*****************************************************************************) PROCEDURE Maurer ( N : INTEGER; Q, K : LONGINT; r, s, L : CARDINAL ); (* 1) Le test est base sur l'article publie par Ueli Maurer: - Maurer, U. (1990): A universal statistical test for random bit generator, Crypto'90, Santa Barbara, California, pp.401-413. 2) On construit des suites de n bits en extrayant s bits d'un nombre produit par un generateur. Pour chaque suite, on calcule la moyenne du logBase2 de la plus recente apparition d'un bloc de L bits. On applique un test de KS aux N moyennes observees. 3) Pour une implantation sous MS-DOS, L est limite par le nombre de bits representant l'indice maximum (16382) d'un tableau d'elements de LONGINT (tableau indice a partir de 0), donc L <= 13 . Sur les SUN, la limite est fixee arbitrairement par MaxIndex4. *) VAR Esp : ARRAY [1..32] OF LONGREAL; (* l'esperance theorique du logBase2 de*) (* la plus recente apparition d'un bloc*) KVar : ARRAY [1..32] OF LONGREAL;(* K fois la variance theorique de la *) (* meme quantite *) MoyenneNormee: FlexTab; (* tableau de la moyenne des appari- *) (* tions observees pour chaque suite *) Seq : INTEGER; (* Numero de la suite courante *) LdivS, (* L DIV s *) deuxExpS, (* 2^s selectionne s bits d'un nbre ale*) sBits, (* valeur numerique des s bits choisis *) ii, Nblocs, (* 2^L = Nbre total de blocs distincts *) bloc : CARDINAL; n, (* nbre total de bits dans une suite *) rangBloc : LONGINT; LCdeuxExpS, (* deuxExpS en LONGCARD *) rsBits : LONGCARD; (* nbre aleatoire de 32 bits genere *) Base, (* 2^(r + s) *) ARang, (* plus recente apparition en LONGREAL *) moyenne, (* moyenne observee pour une suite *) somme, sigma : LONGREAL; (* ecart-type theorique = rac(Variance)*) PROCEDURE ProcInit (); BEGIN (* L'esperance theorique de la moyenne et de la variance pour differentes valeurs de L, telles que calculees dans Maurer. *) Esp[1] := StringToLR ("0.73264948"); KVar[1] := StringToLR("0.68977"); Esp[2] := StringToLR ("1.53743829"); KVar[2] := StringToLR("1.33774"); Esp[3] := StringToLR ("2.40160681"); KVar[3] := StringToLR("1.90133"); Esp[4] := StringToLR ("3.31122472"); KVar[4] := StringToLR("2.35774"); Esp[5] := StringToLR ("4.25342659"); KVar[5] := StringToLR("2.70455"); Esp[6] := StringToLR ("5.21770525"); KVar[6] := StringToLR("2.95403"); Esp[7] := StringToLR ("6.19625065"); KVar[7] := StringToLR("3.12539"); Esp[8] := StringToLR ("7.18366555"); KVar[8] := StringToLR("3.23866"); Esp[9] := StringToLR ("8.17642476"); KVar[9] := StringToLR("3.31120"); Esp[10] := StringToLR ("9.17232431"); KVar[10] := StringToLR("3.35646"); Esp[11] := StringToLR ("10.1700323"); KVar[11] := StringToLR("3.38409"); Esp[12] := StringToLR ("11.1687649"); KVar[12] := StringToLR("3.40065"); Esp[13] := StringToLR ("12.1680703"); KVar[13] := StringToLR("3.41043"); Esp[14] := StringToLR ("13.1676926"); KVar[14] := StringToLR("3.41614"); Esp[15] := StringToLR ("14.1674884"); KVar[15] := StringToLR("3.41943"); Esp[16] := StringToLR ("15.1673788"); KVar[16] := StringToLR("3.42130"); IF L > 16 THEN Esp[L] := IToLR(INTEGER(L)) - 0.832746; (* pour L > 16 et proche de 16, la 6-ieme decimale de Esp[L] pourrait etre inexacte *) KVar[L] := 3.423715; (* ...pour KVar[L], la 4-ieme decimale *) END; sigma := Sqrt( KVar[L] /LIToLR(K)); deuxExpS := TRUNC (TwoExp [s]); LCdeuxExpS := VAL( LONGCARD, deuxExpS); Base := TwoExp [r + s]; Nblocs := TRUNC (TwoExp [L]); LdivS := L DIV s; NewFlexTab( MoyenneNormee, R ); END ProcInit; PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput ); *) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' Maurer (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(Q,1); WriteString(', '); WriteLongInt(K,1); WriteString(', '); WriteCard(r,1); WriteString(', '); WriteCard(s,1); WriteString(', '); WriteCard(L,1); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; BEGIN Entete ('Maurer Test'); n := (K + Q) * IToLI(L); WriteString ('N = '); WriteInt (N, 4); WriteString (' suites de '); WriteLongInt( n, 6); WriteString (' bits;'); WriteLn; WriteString ('Q = '); WriteLongInt (Q, 4); WriteString (" blocs d'initialisation; "); WriteLn; WriteString ('K = '); WriteLongInt (K, 4); WriteString (' blocs pour le test;'); WriteLn; WriteString ('r = '); WriteCard (r, 2); WriteString (' ; s = '); WriteCard (s, 2); WriteLn; WriteString ('Blocs de L = '); WriteCard(L,2); WriteString (' bits; '); WriteLn; IF ((r + s) > 31) THEN MessageErreur (111) END; (*SUN*) (* IF ( L > 13 ) THEN MessageErreur (112) END;*)(* sur PC *) IF ( L > 18) THEN MessageErreur (113) END; (*SUN*) IF ( L MOD s # 0 ) THEN MessageErreur (110) END; ProcInit (); CreateTab (N, MoyenneNormee); CreateTab (Nblocs, Count); IF Detailed THEN WriteString ("Valeur de l'apparition moyenne normalisee a N(0,1) pour"); WriteString (' chaque suite:'); WriteLn; WriteLn END; (******** Debut du test ********) FOR Seq := 1 TO N DO (* ecrire le numero de la sequence a l'ecran *) IF surMoniteur THEN Moniteur () END; (* Initialisation des tableaux *) WITH Count^ DO FOR ii := 0 TO Nblocs - 1 DO Vli^[ii] := 0 END END; (* initialisation avec Q blocs *) FOR rangBloc := 0 TO (Q - 1) DO (* Generation d'un bloc *) bloc := 0; FOR ii := 1 TO LdivS DO (* Genere un nombre alea.; garde seulement les r+s premiers bits *) rsBits := TRUNC (Base * Uniforme01 ()); (* elimination des r bits de gauche *) sBits := VAL (CARDINAL, (rsBits MOD LCdeuxExpS)); (* Decalage de s bits et insertion des s nouveaux bits *) bloc := bloc * deuxExpS + sBits; END; Count^.Vli^[bloc] := rangBloc; END; (****** test proprement dit *******) somme := ZeroLR; FOR rangBloc := Q TO (Q + K - 1) DO (* Generation d'un bloc *) bloc := 0; FOR ii := 1 TO LdivS DO (* Genere un nbre alea.; garde seulement les r+s premiers bits *) rsBits := TRUNC (Base * Uniforme01 ()); (* elimination des r bits de gauche *) sBits := VAL (CARDINAL, (rsBits MOD LCdeuxExpS)); (* Decalage de s bits et insertion des s nouveaux bits *) bloc := bloc * deuxExpS + sBits; END; ARang := LIToLR (rangBloc - Count^.Vli^[bloc]); somme := somme + Ln (ARang); Count^.Vli^[bloc] := rangBloc; END; moyenne := UnSurLnDeux * somme / LIToLR(K); MoyenneNormee^.Vr^[Seq] := (moyenne - Esp[L]) / sigma; IF Detailed THEN WriteLR (MoyenneNormee^.Vr^[Seq] , 7, 2, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN (* On ecrit 10 valeurs de Moyenne par*) WriteLn (* ligne dans le fichier de resultat *) END; END; END; (* Test de K-S sur la repartition des valeurs de MoyenneNormee. *) Stat := StatKS (FctnRepartNormale, 99.9 , MoyenneNormee, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); (* Destruction des tableaux utilises pendant le test *) DeleteTab (Count); DeleteTab( MoyenneNormee); END Maurer; (****************************************************************************) PROCEDURE LongestHeadRun ( N, n :INTEGER; l: LONGINT; r, s : CARDINAL ); CONST jmax = 1000; (* Nombre maximum de categories *) lmin = 1000; (* Longueur minimale d'une suite = *) (* nombre de bits *) TYPE BinSet = SET OF [0..31]; VAR Rep, (* Numero de la repetition courante *) Seq : INTEGER; (* Numero de la suite courante *) lIter, i : LONGINT; (* lIter: nombre d'iterations *) nLR, (* n en LONGREAL *) TroisDemies, (* 3/2 *) Lnl, (* ln l *) esperance, (* Esperance theorique de la plus *) (* longue serie de 1 consecutifs *) KhiRep, Diff, (* Variables temporaires *) rsBits : LONGREAL; (* 2^(r+s) = decalage de r+s bits *) ensemble : BinSet; (* Ensemble de bits choisis de U(0,1)*) longueur, (* Longueur temporaire *) longest, (* Plus longue serie de 1 a date *) sm1, (* s - 1 *) jlow, jhigh, (* Indices limites des categories *) j, jmoy : CARDINAL; Prob : ARRAY [0..jmax] OF LONGREAL; (* Prob[j] = probabilite *) (* que la plus longue serie de 1 soit de longueur j *) (******************) PROCEDURE Moniteur(); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput );*) WriteInt(Rep, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' LongestHeadRun (' ); WriteInt(N,1); WriteString(', '); WriteInt(n,1); WriteString(', '); WriteLongInt(l,1); WriteString(', '); WriteCard(r,1); WriteString(', '); WriteCard(s,1); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (******************) PROCEDURE ProbabiliteLHR (j : CARDINAL): LONGREAL; (* Retourne la probabilite que la plus longue serie de 1 consecutifs soit*) (* de longueur j. -> Voir references et notes manuscrites *) VAR temp, x : LONGREAL; BEGIN temp := (IToLR (INTEGER(j + 1))) * LnDeux - Lnl; x := Exp (- Exp (- temp)); temp := temp + LnDeux; x := Exp (- Exp (- temp)) - x; RETURN x END ProbabiliteLHR; (******************) PROCEDURE ProcInit (); BEGIN TroisDemies := StringToLR ("1.5"); Lnl := Ln (LIToLR (l)); nLR := IToLR(n); rsBits := TwoExp [r + s]; sm1 := s - 1; lIter := 1 + l DIV IToLI (s); esperance := (StringToLR ("0.5772156") + Lnl) / LnDeux - TroisDemies; END ProcInit; (******************) PROCEDURE CalculProb (); (* Calcule les probabilites telles que le nombre espere de suites dans *) (* chaque categorie soit plus grand que MinExpected *) VAR jmoy, j : CARDINAL; somme, temp : LONGREAL; BEGIN jmoy := TRUNC (esperance); (* Indice de categorie moyenne *) (* tant que nombre espere >= MinExpected, une longueur = une categorie*) (* trouver ainsi indice jlow *) j := jmoy + 1; somme := ZeroLR; REPEAT DEC(j); Prob[j] := ProbabiliteLHR (j); somme := somme + Prob[j]; UNTIL (Prob[j] * nLR < MinExpected) OR (j = 0); IF Prob[jmoy] * nLR < MinExpected THEN MessageErreur (132) END; (* condenser toutes les categories de longueur <= jlow dans une seule:*) (* j = jlow; voir fichiers long3?.res pour exemples de probabilites *) jlow := j + 1; somme := somme - Prob[j]; temp := Prob[j]; REPEAT DEC(j); Prob[j] := ProbabiliteLHR (j); temp := temp + Prob[j]; UNTIL Prob[j] < StringToLR("0.0000001"); Prob[jlow] := Prob[jlow] + temp; somme := somme + temp; (* tant que nombre espere >= MinExpected, une longueur = une categorie*) (* trouver indice jhigh *) j := jmoy; REPEAT INC(j); Prob[j] := ProbabiliteLHR (j); somme := somme + Prob[j]; UNTIL Prob[j] * nLR < MinExpected; jhigh := j; (* condenser toutes les categories de longueur >= jhigh dans une *) (* seule categorie: j = jhigh; voir fichiers long3?.res *) Prob[jhigh] := Prob[jhigh] + OneLR - somme; FOR j := jlow TO jhigh DO NbExpected^.Vr^[j] := nLR * Prob[j] (* Nombre espere par categorie *) END; END CalculProb; (******************) BEGIN Entete ('Test LongestHeadRun'); WriteString (' N = '); WriteInt (N, 1); WriteString ('; n = '); WriteInt (n, 1); WriteString ('; l = '); WriteLongInt(l, 1); WriteString ('; r = '); WriteCard (r, 1); WriteString ('; s = '); WriteCard (s, 1); WriteLn; WriteLn; IF r + s > 31 THEN MessageErreur (130) END; (*SUN*) IF (LONGCARD (BinSet{0}) # 1) OR (LONGCARD (BinSet{31}) # 80000000H) THEN MessageErreur (131) END; IF l < lmin THEN MessageErreur (133) END; CreateTab (jmax, NbExpected); CreateTab (N+1, Khi); CreateTab (jmax, Count); ProcInit; CalculProb; IF Detailed THEN WriteString (' Esperance theorique de la plus longue suite: '); WriteLR (esperance, 7, 2, 2); WriteLn; WriteLn; WriteString (" Le nombre des plus longues suites de j '1' pour "); WriteString ('chaque sequence:'); WriteLn; WriteLn; WriteString(" j nombre nombre"); WriteLn; WriteString(" espere mesure"); WriteLn; WriteLn; END; FOR Rep := 1 TO N DO (* Debut du test *) IF surMoniteur THEN Moniteur() END; FOR j := jlow TO jhigh DO Count^.Vi^[j] := 0 (* Initialiser chaque categorie a 0 *) END; FOR Seq := 1 TO n DO (* Il y a n suites de l bits pour *) longest := 0; (* chaque Rep *) longueur := 0; (* En fait, les suites Seq sont constituees en tout de l + l MOD s *) (* bits au lieu de l bits -> implantation plus simple et influence *) (* sur les resultats tout a fait negligeable. *) FOR i := 1 TO lIter DO (* On construit une suite de l bits *) (* raison du MessageErreur (131) ici *) ensemble := BinSet (TRUNC (rsBits * Uniforme01 ())); FOR j := 0 TO sm1 DO (* Examiner chaque bit *) IF j IN ensemble THEN INC (longueur) ELSE IF longueur > longest THEN longest := longueur END; longueur := 0 END; END; END; IF longueur > longest THEN longest := longueur END; IF longest <= jlow THEN INC (Count^.Vi^[jlow]) ELSIF longest >= jhigh THEN INC (Count^.Vi^[jhigh]) ELSE INC (Count^.Vi^[longest]) END; END; (* calculer une valeur de KhiDeux pour les n suites de l bits en com- *) (* parant les nombres esperes avec ceux mesures pour chaque categorie *) KhiRep := ZeroLR; FOR j := jlow TO jhigh DO Diff := IToLR (Count^.Vi^[j]) - NbExpected^.Vr^[j]; KhiRep := KhiRep + Diff * Diff / NbExpected^.Vr^[j] END; Khi^.Vr^[Rep] := KhiRep; (* une valeur de KhiDeux obtenue *) IF Detailed THEN FOR j := jlow TO jhigh DO IF j = jlow THEN WriteString (" <= ") ELSIF j = jhigh THEN WriteString (" >= ") ELSE WriteString (" = ") END; WriteCard (j, 3); WriteLR (NbExpected^.Vr^[j], 12, 2, 2); WriteInt (Count^.Vi^[j], 12); WriteLn; END; WriteLn; WriteLn; END; END; (* Test de K-S sur la repartition des N valeurs de KhiDeux obtenues *) Stat := StatKS (FctnRepartKhiDeux, IToLR (jhigh - jlow) , Khi, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); DeleteTab (NbExpected); DeleteTab (Khi); DeleteTab (Count); END LongestHeadRun; (****************************************************************************) PROCEDURE gTransform ( N:INTEGER; n:LONGINT; g:LONGREAL ); CONST jmax = 1000; (* Nombre maximal de categories *) gmax = 10.0; (* Valeur maximale de g *) VAR Seq, (* Numero de la sequence courante *) gTrunc, (* Plancher de g = TRUNC (g) *) temp, (* Variable temporaire *) j, (* Indice d'une categorie *) jlow, (* Indice de la plus basse categorie *) jhigh : INTEGER; (* Indice de la plus haute categorie *) i : LONGINT; nLR, (* n en LONGREAL *) Diff, KhiSeq, (* Variables temporaires *) v : LONGREAL; (* La nouvelle variable aleatoire *) (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput );*) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' gTransform (' ); WriteInt(N,1); WriteString(', '); WriteLongInt(n,1); WriteString(', '); WriteLongReal( g, 4, 6); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProbabiliteG ( j:INTEGER ): LONGREAL; (* Retourne la probabilite que le nombre minimum de nombres aleatoires *) (* dont la somme est plus grande que g egale j+1. -> voir reference *) VAR r :INTEGER; somme, signe, (* +1 ou -1 *) jLR, (* j en LONGREAL *) temp :LONGREAL; BEGIN somme := ZeroLR; signe := OneLR; jLR := IToLR (j); FOR r := 0 TO gTrunc DO temp := signe * Combinaison (r, j+1) ; temp := temp * Pow (g - IToLR(r), jLR); somme := somme + temp; signe := -signe; END; somme := (jLR + OneLR - g) * somme / Factorielle (j + 1); RETURN somme; END ProbabiliteG; (*****************) PROCEDURE CalculNbExpected ( nLR :LONGREAL); (* Calcule le nombre espere dans chaque categorie. *) VAR j : INTEGER; temp, sommeProb : LONGREAL; BEGIN sommeProb := ZeroLR; j := gTrunc - 1; (* reunir toutes les categories inferieures dans une categorie jlow *) (* telle que le nombre espere soit plus grand que MinExpected. *) REPEAT INC(j); temp := ProbabiliteG (j); sommeProb := sommeProb + temp; UNTIL sommeProb >= MinExpected / nLR; jlow := j; NbExpected^.Vr^[jlow] := sommeProb * nLR; (* pour les categories moyennes, 1 valeur de j = 1 categorie *) REPEAT INC(j); temp := ProbabiliteG (j); sommeProb := sommeProb + temp; NbExpected^.Vr^[j] := temp * nLR; UNTIL NbExpected^.Vr^[j] < MinExpected; (* reunir toutes les categories superieures dans une categorie jhigh *) (* telle que le nombre espere soit plus grand que MinExpected. *) jhigh := j - 1; NbExpected^.Vr^[jhigh] := NbExpected^.Vr^[jhigh] + temp * nLR; temp := OneLR - sommeProb; NbExpected^.Vr^[jhigh] := NbExpected^.Vr^[jhigh] + temp * nLR; END CalculNbExpected; (*****************) BEGIN Entete ('Test gTransform'); WriteString (' N = '); WriteInt (N, 4); WriteString ('; n = '); WriteLongInt (n, 4); WriteString ('; g = '); WriteLongReal (g, 4, 6); WriteLn; WriteLn; IF (g < OneLR) OR (g > gmax) THEN MessageErreur (140) END; CreateTab (N+1, Khi); CreateTab (jmax, NbExpected); CreateTab (jmax, Count); gTrunc := TRUNC (g); nLR := LIToLR (n); CalculNbExpected (nLR); IF Detailed THEN WriteString (" j nombre nombre"); WriteString (" Statistique KhiDeux"); WriteLn; WriteString (" espere mesure"); WriteLn; WriteLn; END; (* Debut du test *) FOR Seq := 1 TO N DO IF surMoniteur THEN Moniteur() END; FOR j := 1 TO jmax DO (* Initialise le nombre de chaque *) Count^.Vi^[j] := 0 (* categorie a 0 *) END; FOR i := 1 TO n DO v := ZeroLR; j := 0; REPEAT v := v + Uniforme01() ; INC (j); UNTIL v > g; INC (Count^.Vi^[j - 1]); END; temp := 0; FOR j := 1 TO jlow DO temp := temp + Count^.Vi^[j] END; Count^.Vi^[jlow] := temp; temp := 0; FOR j := jhigh TO jmax DO temp := temp + Count^.Vi^[j] END; Count^.Vi^[jhigh] := temp; IF Detailed THEN FOR j := jlow TO jhigh DO IF j = jlow THEN WriteString (" <= ") ELSIF j = jhigh THEN WriteString (" >= ") ELSE WriteString (" = ") END; WriteInt (j, 3); WriteLR (NbExpected^.Vr^[j], 17, 2, 2); WriteInt (Count^.Vi^[j], 16); WriteLn; END; END; (* Calcul de la statistique du *) (* Khi-Deux pour une sequence. *) KhiSeq := ZeroLR; FOR j := jlow TO jhigh DO Diff := IToLR (Count^.Vi^[j]) - NbExpected^.Vr^[j]; KhiSeq := KhiSeq + Diff * Diff / NbExpected^.Vr^[j] END; Khi^.Vr^[Seq] := KhiSeq; IF Detailed THEN WriteLR (KhiSeq, 70, 2, 2); WriteLn; WriteLn; WriteString ('---------------------------------'); WriteLn; END; END; (* Test de K-S sur la repartition des N valeurs de KhiDeux obtenues; *) (* nombre de degres de liberte = nombre de categories - 1 = jhigh - jlow *) Stat := StatKS (FctnRepartKhiDeux, IToLR (jhigh - jlow) , Khi, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); DeleteTab (NbExpected); DeleteTab (Khi); DeleteTab (Count); END gTransform; (****************************************************************************) PROCEDURE BirthdaySpacings (N, k, m :INTEGER; n :LONGINT); CONST hache1 = 15901; (* Dimension de la table de hachage *) hache2 = 41; (* Hachage secondaire si collision *) jmax = 8000; (* Entier maximal pour la distribu- *) (* tion de Poisson *) VAR jlow, (* Indice de la plus basse categorie*) jhigh, (* Indice de la plus haute categorie*) Y, (* Nombre de collisions *) liberte, (* Nombre de degres liberte du Khi2 *) Seq, Rep, (* Compteurs ... *) i, j : INTEGER; (* ... *) tempLI : LONGINT; (* Variable temporaire *) nLR, (* n en LONGREAL *) mLR, (* m en LONGREAL *) kLR, (* k en LONGREAL *) lambda : LONGREAL; (* Parametre de Poisson *) Differences, (* Differences des dates *) Dates : FlexTab; (* "Dates de naissance" *) Prob : ARRAY [0..jmax] OF LONGREAL; (* Probabilite pour chaque categorie de Poisson *) (*****************) PROCEDURE Moniteur (); BEGIN GetOutput( fichierTemp ); (* SetOutput( StandardOutput );*) WriteInt(Seq, 4); WriteString(' / '); WriteInt(N, 3); WriteString( ' BirthdaySpacings (' ); WriteInt(N,1); WriteString(', '); WriteInt(k,1); WriteString(', '); WriteInt(m,1); WriteString(', '); WriteLongInt(n,1); WriteString(' )'); WriteLn; SetOutput( fichierTemp ); END Moniteur; (*****************) PROCEDURE ProcInit (); BEGIN nLR := LIToLR (n); mLR := IToLR (m); kLR := IToLR (k); lambda := mLR * mLR / nLR * mLR / FourLR ; IF lambda > IToLR(jmax) THEN MessageErreur(151) END; NewFlexTab (Dates, LI); NewFlexTab (Differences, LI); CreateTab (N+1, Khi); CreateTab (jmax + 100, NbExpected); CreateTab (jmax + 100, Count); CreateTab (m + 1, Dates); CreateTab (hache1, Differences); END ProcInit; (*****************) PROCEDURE QuickSort (VAR Table:FlexTab; mm :INTEGER); (* Trie les dates generees en ordre croissant *) PROCEDURE Sort (l, r : INTEGER); VAR i, j : INTEGER; x, w : LONGINT; BEGIN (* Sort *) i := l; j := r; x := Table^.Vli^[(l+r) DIV 2]; REPEAT WHILE (Table^.Vli^[i] < x) DO i := i + 1 END; WHILE (x < Table^.Vli^[j]) DO j := j - 1 END; IF (i <= j) THEN w := Table^.Vli^[i]; Table^.Vli^[i] := Table^.Vli^[j]; Table^.Vli^[j] := w; i := i + 1; j := j - 1 END; UNTIL (i > j); IF (l < j) THEN Sort (l, j) END; IF (i < r) THEN Sort (i, r) END END Sort; BEGIN (* QuickSort *) Sort (1, mm) END QuickSort; (*****************) PROCEDURE Poisson (lambda : LONGREAL); CONST limite = 0.0000001; VAR somme : LONGREAL; jmoy, (* Indice de la categorie avec pro- *) j : INTEGER; (* babilite maximale pour Poisson *) BEGIN jmoy := TRUNC (lambda); IF jmoy > jmax THEN MessageErreur (151) END; (* calculer d'abord la probabilite au maximum de Poisson *) somme := ZeroLR; FOR j := 2 TO jmoy DO (* Ln (Factorielle (jmoy)) *) somme := somme + Ln(IToLR(j)) END; somme := - lambda + IToLR(jmoy)*Ln(lambda) - somme; Prob[jmoy] := Exp(somme); (* Probabilite au maximum *) (* Trouver la categorie inferieure: soit j = 0, soit la premiere a *) (* gauche du maximum pour laquelle Probabilite < limite *) j := jmoy; WHILE (j > 0) AND (Prob[j] > limite) DO DEC(j); Prob[j] := Prob[j+1] * IToLR(j+1) / lambda; END; jlow := j; (* Proceder de meme pour les categories superieures *) j := jmoy; REPEAT INC(j); Prob[j] := Prob[j-1] * lambda / IToLR(j); UNTIL (j = jmax) OR (Prob[j] < limite); jhigh := j - 1; (* Calculer les nombres esperes dans chaque categorie *) FOR j := jlow TO jhigh DO NbExpected^.Vr^[j] := kLR * Prob[j] END; END Poisson; (******************) PROCEDURE Hachage (): INTEGER; (* Hachage double non ordonne: placer les differences de dates dans la table de hachage et retourner le nombre de collisions. *) VAR position, y, j, cle : INTEGER; BEGIN y := 0; WITH Differences^ DO FOR j := 1 TO m DO cle := Dates^.Vli^[j]; position := cle MOD hache1; LOOP IF Vli^[position] < 0 THEN Vli^[position] := cle; EXIT END; IF Vli^[position] = cle THEN INC(y); EXIT END; position := position - 23 - cle MOD hache2; IF position < 0 THEN position := position + hache1 END; END; END; END; RETURN y END Hachage; (******************) PROCEDURE CalculKhi2 (); (* Calculer une valeur de KhiDeux en comparant les nombres esperes avec *) (* ceux mesures pour chaque j; regrouper les categories adjacentes en une*) (* grosse categorie telle que son nombre espere soit >= MinExpected *) VAR mesure, mtemp, j : INTEGER; espere, etemp, Diff, KhiSeq : LONGREAL; BEGIN KhiSeq := ZeroLR; mesure := 0; espere := ZeroLR; j := jlow; liberte := -1; WHILE j <= jhigh DO etemp := espere; mtemp := mesure; mesure := 0; espere := ZeroLR; WHILE (j <= jhigh) AND (espere < MinExpected) DO espere := espere + NbExpected^.Vr^[j]; mesure := mesure + Count^.Vi^[j]; INC(j); END; (* Si la derniere categorie espere moins que MinExpected, la regrou-*) (* per avec l'avant-derniere. *) IF espere >= MinExpected THEN INC(liberte); ELSE KhiSeq := KhiSeq - Diff * Diff / etemp; espere := espere + etemp; mesure := mesure + mtemp END; Diff := IToLR(mesure) - espere; KhiSeq := KhiSeq + Diff * Diff / espere; END; Khi^.Vr^[Seq] := KhiSeq; (* une valeur de KhiDeux obtenue *) IF Detailed THEN WriteLR (KhiSeq, 7, 2, 2); WriteString (' '); IF (Seq MOD 10 = 0) THEN WriteLn END; END END CalculKhi2; (******************) BEGIN Entete ('Test BirthdaySpacings'); WriteString (' N = '); WriteInt (N, 1); WriteString ('; k = '); WriteInt (k, 1); WriteString ('; m = '); WriteInt(m, 1); WriteString ('; n = '); WriteLongInt (n, 1); IF m > 40000 THEN MessageErreur(150) END; IF k < 20 THEN MessageErreur(152) END; ProcInit (); WriteString ('; lambda = '); WriteLR(lambda, 6,4,2); WriteLn; WriteLn; IF Detailed THEN WriteString ('Valeur de la statistique pour chaque sequence :'); WriteLn; WriteLn END; Poisson (lambda); FOR Seq := 1 TO N DO (* Debut du test *) IF surMoniteur THEN Moniteur() END; FOR j := jlow TO jhigh DO (* Initialiser les compteurs de *) Count^.Vi^[j] := 0; (* chaque categorie *) END; FOR Rep := 1 TO k DO WITH Differences^ DO (* Initialiser la table de hachage *) FOR i := 0 TO hache1 - 1 DO Vli^[i] := -1; END; END; WITH Dates^ DO (* Generer les "dates de naissance" *) FOR i := 1 TO m DO Vli^[i] := TRUNC (nLR * Uniforme01 ()); END; END; QuickSort (Dates, m); (* Trier les dates *) WITH Dates^ DO (* Calculer la difference des dates *) tempLI := Vli^[1]; (* adjacentes *) FOR i := 1 TO m-1 DO Vli^[i] := Vli^[i+1] - Vli^[i]; END; Vli^[m] := tempLI + n - Vli^[m]; END; Y := Hachage (); (* Y = nombre de collisions *) IF Y <= jlow THEN INC (Count^.Vi^[jlow]) ELSIF Y >= jhigh THEN INC (Count^.Vi^[jhigh]) ELSE INC (Count^.Vi^[Y]) END; END; CalculKhi2; END; (* Test de K-S sur la repartition des N valeurs de KhiDeux obtenues *) Stat := StatKS (FctnRepartKhiDeux, IToLR (liberte) , Khi, N); WriteResults(Stat, OneLR - FctnRepartKS (IToLR(N), Stat)); DeleteTab (NbExpected); DeleteTab (Khi); DeleteTab (Count); DeleteTab (Dates); DeleteTab (Differences) END BirthdaySpacings; (****************************************************************************) BEGIN DixExpMoinsVingt := StringToLR ("10.0E-20"); MinExpected := StringToLR ("10.0"); Detailed := FALSE; surMoniteur := FALSE; LnDeux := Ln (TwoLR); UnSurLnDeux := OneLR / LnDeux; CorrCalculees := 0; Uniforme01 := NotInit; NewFlexTab (NbExpected, R); NewFlexTab (Khi, R); NewFlexTab (Count, I); NewFlexTab (Loca, I); ALLOCATE( DataFRNormale, SIZE(RecordFRNormale)); InitFRNormale( DataFRNormale); ALLOCATE( DataFRKS, SIZE(RecordFRKS)); InitFRKS( DataFRKS); ALLOCATE( DataFRKhiDeux, SIZE(RecordFRKhiDeux)); InitFRKhiDeux( DataFRKhiDeux); END STATUNIF.