C ALGORITHM 823, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 29,NO. 2, June, 2003, P. 95--109. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran77/ # Fortran77/Dp/ # Fortran77/Dp/Sfaur/ # Fortran77/Dp/Sfaur/Drivers/ # Fortran77/Dp/Sfaur/Drivers/Makefile # Fortran77/Dp/Sfaur/Drivers/data1 # Fortran77/Dp/Sfaur/Drivers/driver.f # Fortran77/Dp/Sfaur/Drivers/res1 # Fortran77/Dp/Sfaur/Src/ # Fortran77/Dp/Sfaur/Src/extra.f # Fortran77/Dp/Sfaur/Src/src.f # Fortran77/Dp/Snied/ # Fortran77/Dp/Snied/Drivers/ # Fortran77/Dp/Snied/Drivers/Makefile # Fortran77/Dp/Snied/Drivers/data1 # Fortran77/Dp/Snied/Drivers/driver.f # Fortran77/Dp/Snied/Drivers/res1 # Fortran77/Dp/Snied/Src/ # Fortran77/Dp/Snied/Src/extra.f # Fortran77/Dp/Snied/Src/src.f # Fortran77/Dp/Sniedx/ # Fortran77/Dp/Sniedx/Drivers/ # Fortran77/Dp/Sniedx/Drivers/Makefile # Fortran77/Dp/Sniedx/Drivers/driver.f # Fortran77/Dp/Sniedx/Drivers/nxs10m30 # Fortran77/Dp/Sniedx/Drivers/nxs11m30 # Fortran77/Dp/Sniedx/Drivers/nxs12m30 # Fortran77/Dp/Sniedx/Drivers/nxs13m30 # Fortran77/Dp/Sniedx/Drivers/nxs14m30 # Fortran77/Dp/Sniedx/Drivers/nxs15m30 # Fortran77/Dp/Sniedx/Drivers/nxs16m30 # Fortran77/Dp/Sniedx/Drivers/nxs4m30 # Fortran77/Dp/Sniedx/Drivers/nxs4m30.old # Fortran77/Dp/Sniedx/Drivers/nxs5m30 # Fortran77/Dp/Sniedx/Drivers/nxs6m30 # Fortran77/Dp/Sniedx/Drivers/nxs7m30 # Fortran77/Dp/Sniedx/Drivers/nxs8m30 # Fortran77/Dp/Sniedx/Drivers/nxs9m30 # Fortran77/Dp/Sniedx/Drivers/res1 # Fortran77/Dp/Sniedx/Src/ # Fortran77/Dp/Sniedx/Src/extra.f # Fortran77/Dp/Sniedx/Src/src.f # Fortran77/Dp/Ssobol/ # Fortran77/Dp/Ssobol/Drivers/ # Fortran77/Dp/Ssobol/Drivers/Makefile # Fortran77/Dp/Ssobol/Drivers/data1 # Fortran77/Dp/Ssobol/Drivers/driver.f # Fortran77/Dp/Ssobol/Drivers/res1 # Fortran77/Dp/Ssobol/Src/ # Fortran77/Dp/Ssobol/Src/extra.f # Fortran77/Dp/Ssobol/Src/src.f # This archive created: Mon Jul 7 15:01:40 2003 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Sfaur' then mkdir 'Sfaur' fi cd 'Sfaur' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' include makefile.inc all: Res1 SRCLIBS= $(LINPACK) $(EISPACK_DP) $(BLAS) $(PORT) src.o: $(INCSRC)/src.f $(F77) $(F77OPTS) -c $(INCSRC)/src.f extra.o: $(INCSRC)/extra.f $(F77) $(F77OPTS) -c $(INCSRC)/extra.f DRIVERS= driver RESULTS= Res1 SRCLIBS = $(TIMER) Objs1= driver.o src.o extra.o driver: $(Objs1) $(F77LINK) $(F77LINKOPTS) -o driver $(Objs1) $(SRCLIBS) Res1: driver ./driver >Res1 diffres: Res1 res1 echo "Differences in results from driver" $(DIFF) Res1 res1 clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << "SHAR_EOF" > 'data1' 683644482 135580830 114750887 803113694 903684270 72157854 59862132 505680715 779216195 193586105 329775408 1049470377 790869458 529177858 770596904 1072006936 140712293 266923202 844524433 114947652 963169856 269847095 96724448 724232372 729133779 872385703 803906910 638059470 526484481 843345462 431170622 1057600395 413165996 637333254 738214251 987594285 394318066 718140016 1064908418 636856468 856354958 404426257 6726334 916596837 284819911 905844483 331078851 514311810 280380871 612727606 184724875 578053121 197887820 50641495 423123265 227871988 36070792 465200006 187446371 99681590 983367218 414604440 243776932 990984772 289169669 72196638 522106132 18453143 1008259046 521211985 88531465 816692680 1001849242 432014653 124800720 554138134 166704786 634462798 1019392879 691725817 459507961 1043416360 536793729 1046579554 403993816 412804346 138344837 24248109 184604157 535890993 869278100 125797533 322976214 161488107 755812753 541238672 841393456 909350792 521628897 567331204 118313388 265208159 468752989 169896169 959104945 500814545 434091372 113024960 444725900 426474642 263586462 432120698 45230250 524531296 327681208 288240183 353741403 192041686 42986304 258035065 388012188 537416171 292523858 794840680 502322834 423539584 884204171 419181022 618098449 257684411 288733678 145066288 756812965 357038497 881037192 314634603 182012383 753118967 636654128 318546833 71382991 400554255 153755763 187747483 306343731 1044780 402486338 564947470 934172722 813806089 705671218 795868370 926341726 51211860 739601563 176287958 1016675820 950546833 267929114 275671035 666917527 431607592 177069271 570312219 125215260 886835363 75079442 767034988 879903592 525306349 560816857 973795128 774079162 143751759 137446578 209360240 1012892947 398130103 838233295 47600869 312065869 24489594 686091463 1064468817 610953508 593516621 736434891 331884657 295796744 991620199 197203105 694133440 834981040 813915471 1043396296 246522553 741863988 188767345 707349375 236730174 943797129 1073073438 917394477 81723820 598848363 17981413 322240492 766981617 282924362 643850936 914020733 586984616 128906358 739978294 543175409 351896331 1017052692 802102481 946877967 615436223 305014361 1038674293 372513724 902652345 552583238 507262995 36761397 753350193 116032862 328670174 103551339 1045088961 667010209 737933140 945361853 324875181 353779021 866388681 255506551 65493883 220170356 837239777 771797016 542598812 568193573 359763920 870736829 304824104 430870768 116534859 741311666 157845915 861348937 985097428 714048271 306167010 709889598 498574308 382576036 360956609 1005272948 680195561 845547054 1026217890 374002007 866217449 695335995 462692419 314161108 1048195907 709323042 604523286 581285813 463752482 919585493 518104348 905143615 884163267 99734344 332656669 398822849 973544107 459181032 1001255433 846689598 32962550 6372759 727573859 495342420 738892195 834835824 1047729137 529022179 1032177776 185764687 483299066 818974724 826550698 987239779 684232759 940422073 224992310 164068454 827234589 690302353 32197089 628566875 95166490 866947171 476745989 604118539 121738255 667653407 987681310 749263282 863528381 341372968 244006930 951984484 614173490 685305366 592527405 383027710 952965461 21271273 653806680 745437736 180332506 331728868 207307411 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' PROGRAM TSFAUR C C C THIS PROGRAM TESTS ACCURACY OF C NUMERICAL INTEGRATION USING "GNFAUR" C AND INTEGRAND (2) OF DAVIS AND C RABINOWITZ, PAGE 406 C C User Define: C DIMEN : dimension C ATMOST : sequence length C SAMS : Number of replications C MAXS : Maximum Digits of Scrambling Of Owen type Scrambling C IFLAG: User Choice of type Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,II,IFLAG,K,J,MAXS,SAM C .. C .. Array Arguments .. DOUBLE PRECISION OUTS(500,10000) REAL SECOND, T1 C .. C .. C .. External Subroutines .. EXTERNAL GNFAUR,TIME C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. SAM = 1 MAXS = 30 DIMEN = 2 IFLAG = 1 ATMOST = 2**12 DO 30 II = 1,SAM WRITE (*,FMT=*) 'I = ITERATION NUMBER' WRITE (*,FMT=*) 'EI = ESTIMATED INTEGRAL' T1 = SECOND() CALL GNFAUR(DIMEN,ATMOST,IFLAG,MAXS,OUTS) SUM = 0.0 K = 9 DO 20 I = 1,ATMOST F = 1.0 DO 10 J = 1,DIMEN F = F*ABS(4.0*OUTS(J,I)-2.0) 10 CONTINUE IF (MOD(I,2**K).EQ.0) THEN K = K + 1 WRITE (*,FMT=*) 'I = ',I WRITE (*,FMT=*) 'EI = ',SUM/I END IF SUM = SUM + F 20 CONTINUE T1 = SECOND() - T1 WRITE (*,FMT=*) 'Total time elapsed = ',T1 30 CONTINUE STOP END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' I = ITERATION NUMBER EI = ESTIMATED INTEGRAL I = 512 EI = 0.9977403526167160 I = 1024 EI = 0.9997283042920412 I = 2048 EI = 0.9991515530877914 I = 4096 EI = 0.9998810848129468 Total time elapsed = 0.4800000 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'extra.f' then echo shar: will not over-write existing file "'extra.f'" else cat << "SHAR_EOF" > 'extra.f' INTEGER FUNCTION EXOR(IIN,JIN) C C THIS FUNCTION CALCULATES THE EXCLUSIVE-OR OF ITS C TWO INPUT PARAMETERS C C .. Scalar Arguments .. INTEGER IIN,JIN C .. C .. Local Scalars .. INTEGER I,J,K,L C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. I = IIN J = JIN K = 0 L = 1 C 10 IF (I.EQ.J) THEN EXOR = K RETURN END IF C C CHECK THE CURRENT RIGHT-HAND BITS OF I AND J. C IF THEY DIFFER, SET THE APPROPRIATE BIT OF K. C IF (MOD(I,2).NE.MOD(J,2)) K = K + L I = I/2 J = J/2 L = 2*L GO TO 10 END DOUBLE PRECISION FUNCTION UNI() * * Random number generator, adapted from F. James * "A Review of Random Number Generators" * Comp. Phys. Comm. 60(1990), pp. 329-344. * C .. Parameters .. DOUBLE PRECISION TWOM24 PARAMETER (TWOM24=1D0/16777216.0) C .. C .. Local Scalars .. DOUBLE PRECISION CARRY INTEGER I,J C .. C .. Local Arrays .. DOUBLE PRECISION SEEDS(24) C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Save statement .. SAVE I,J,CARRY,SEEDS C .. C .. Data statements .. DATA I,J,CARRY/24,10,0.0d0/ DATA SEEDS/0.8804418d0,0.2694365d0,0.0367681d0,0.4068699d0, +0.4554052d0, + 0.2880635d0,0.1463408d0,0.2390333d0,0.6407298d0, +0.1755283d0,0.7132940d0, + 0.4913043d0,0.2979918d0,0.1396858d0,0.3589528d0, +0.5254809d0,0.9857749d0, + 0.4612127d0,0.2196441d0,0.7848351d0,0.4096100d0, +0.9807353d0,0.2689915d0, + 0.5140357d0/ C .. UNI = SEEDS(I) - SEEDS(J) - CARRY IF (UNI.LT.0) THEN UNI = UNI + 1 CARRY = TWOM24 ELSE CARRY = 0 END IF SEEDS(I) = UNI I = 24 - MOD(25-I,24) J = 24 - MOD(25-J,24) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE GNCRML(LSM,SHIFT,MAXL) C .. Parameters .. INTEGER MAXLN PARAMETER (MAXLN=30) C .. Scalar Arguments .. INTEGER MAXL C .. C .. Array Arguments .. INTEGER LSM(1:500,0:MAXLN,0:MAXLN),SHIFT(1:500,0:MAXLN) C .. C .. Scalars in Common .. DOUBLE PRECISION RQS INTEGER HISUM,NEXTN,QS,S,TESTN C .. C .. Arrays in Common .. INTEGER COEF(1:500,0:MAXLN,0:MAXLN),PGTEMP(0:MAXLN), + SCOEF(1:500,0:MAXLN,0:MAXLN),ZTEMP(1:500,0:MAXLN) C .. C .. Local Scalars .. INTEGER I,J,P,QSM1 C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Common blocks .. COMMON /FAURE/S,QS,HISUM,SCOEF,NEXTN,TESTN,RQS,PGTEMP,ZTEMP,COEF C .. C .. Save statement .. SAVE /FAURE/ C .. QSM1 = QS - 1 DO 30 P = 1,S DO 20 I = 0,MAXL SHIFT(P,I) = MOD((INT(UNI()*1000.0)),QS) DO 10 J = 0,HISUM IF (J.EQ.I) THEN LSM(P,I,J) = MOD((INT(UNI()*1000.0)),QSM1) + 1 ELSE IF (J.LT.I) THEN LSM(P,I,J) = MOD((INT(UNI()*1000.0)),QS) ELSE LSM(P,I,J) = 0 END IF 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GNCRMU(USM,USHIFT,MAXL) C GENERATING LOWER TRIANGULAR SCRMABLING MATRICES AND C SHIFT VECTORS. C INPUTS : C FROM INSOBL : MAXL C FROM BLOCK DATA "SOBOL" : S, MAXCOL, C C OUTPUTS : C TO INSOBL : USM, USHIFT C .. Parameters .. INTEGER MAXLN PARAMETER (MAXLN=30) C .. Scalar Arguments .. INTEGER MAXL C .. C .. Array Arguments .. INTEGER USHIFT(0:MAXLN),USM(0:MAXLN,0:MAXLN) C .. C .. Scalars in Common .. DOUBLE PRECISION RQS INTEGER HISUM,NEXTN,QS,S,TESTN C .. C .. Arrays in Common .. INTEGER COEF(1:500,0:MAXLN,0:MAXLN),PGTEMP(0:MAXLN), + SCOEF(1:500,0:MAXLN,0:MAXLN),ZTEMP(1:500,0:MAXLN) C .. C .. Local Scalars .. INTEGER I,J,QSM1,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Common blocks .. COMMON /FAURE/S,QS,HISUM,SCOEF,NEXTN,TESTN,RQS,PGTEMP,ZTEMP,COEF C .. C .. Save statement .. SAVE /FAURE/ C .. QSM1 = QS - 1 DO 20 I = 0,HISUM STEMP = MOD((INT(UNI()*1000.0)),QS) C STEMP = 0 USHIFT(I) = STEMP DO 10 J = 0,HISUM IF (J.EQ.I) THEN TEMP = MOD((INT(UNI()*1000.0)),QSM1) + 1 C TEMP = 1 ELSE IF (J.GT.I) THEN C TEMP = MOD((int(UNI()*1000.0)),QS) TEMP = 0 ELSE TEMP = 0 END IF USM(I,J) = TEMP 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE GNFAUR(DIMEN,ATMOST,IFLAG,MAXS,OUTS) C C User Define: C DIMEN : dimension C ATMOST : sequence length C MAXS : Maximum Digits of Scrambling Of Owen type Scrambling C IFLAG: User Choice of Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C C .. Scalar Arguments .. INTEGER ATMOST,DIMEN,IFLAG,MAXS C .. C .. Array Argument.. DOUBLE PRECISION OUTS(500,10000) C .. C .. Local Scalars .. INTEGER I,J,MAXL,MAXX C .. C .. Local Arrays .. DOUBLE PRECISION QUASI(500) LOGICAL FLAG(2) C .. C .. External Subroutines .. EXTERNAL GSFAUR,INFAUR C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. MAXL = MAXS - 1 CALL INFAUR(FLAG,DIMEN,ATMOST,QUASI,MAXL,IFLAG,MAXX) IF (.NOT.FLAG(2)) THEN WRITE (*,FMT=*) 'ATMOST = ',ATMOST WRITE (*,FMT=*) 'ATMOST IS NOT OK' STOP END IF DO 10 J = 1,DIMEN OUTS(J,1) = QUASI(J) 10 CONTINUE DO 30 I = 2,ATMOST CALL GSFAUR(QUASI,MAXX) DO 20 J = 1,DIMEN OUTS(J,I) = QUASI(J) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GSFAUR(QUASI,MAXL) C C THIS SUBROUTINE GENERATES A NEW C QUASIRANDOM VECTOR WITH EACH CALL BY USING THE ALPHA-ARY C GRAY CODE METHOD OF LICHTNER "SIAM J. DISCRETE MATH(1998), C 381-386". (SEE ESPECIALLY PAGE 381-382). C C THE USER MUST CALL "INFAUR" BEFORE C CALLING "GSFAUR". C AFTER CALLING "INFAUR", TEST FLAG(1) C AND FLAG(2); IF EITHER IS FALSE, DO C NOT CALL GOFAUR. READ THE COMMENTS AT C THE BEGINNING OF INFAUR AND THEN C THOSE BELOW. C C ALL INPUTS COME FROM "INFAUR" VIA C LABELLED COMMON "FAURE"; FOR THEIR C DEFINITIONS, SEE "INFAUR". C C INPUTS: C S,QS,COEF,NEXTN,TESTN,HISUM,RQS,MAXL C C OUTPUTS: C TO USER'S CALLING PROGRAM: C QUASI - A NEW SCRAMBLED QUASIRANDOM VECTOR C C C C C C NEXTN HAS A REPRESENTATION IN BASE C QS OF THE FORM: SUM OVER J FROM ZERO C TO HISUM OF YTEMP(J)*(QS**J) C C WE NOW COMPUTE THE YTEMP(J)'S SAME AS FAURE C C .. Parameters .. INTEGER MAXLN PARAMETER (MAXLN=30) C C .. Scalar Arguments .. INTEGER MAXL C .. C .. Array Arguments .. DOUBLE PRECISION QUASI(500) C .. C .. Scalars in Common .. DOUBLE PRECISION RQS INTEGER HISUM,NEXTN,QS,S,TESTN C .. C .. Arrays in Common .. INTEGER COEF(1:500,0:MAXLN,0:MAXLN),PGTEMP(0:MAXLN), + SCOEF(1:500,0:MAXLN,0:MAXLN),ZTEMP(1:500,0:MAXLN) C .. C .. Local Scalars .. DOUBLE PRECISION R INTEGER DIFF,I,K,KTEMP,LTEMP,MTEMP,POS,TEM C .. C .. Local Arrays .. INTEGER GTEMP(0:MAXLN),YTEMP(0:19) C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /FAURE/S,QS,HISUM,SCOEF,NEXTN,TESTN,RQS,PGTEMP,ZTEMP,COEF C .. C .. Save statement .. SAVE /FAURE/ C .. KTEMP = TESTN LTEMP = NEXTN DO 10 I = HISUM,0,-1 KTEMP = KTEMP/QS MTEMP = MOD(LTEMP,KTEMP) YTEMP(I) = (LTEMP-MTEMP)/KTEMP LTEMP = MTEMP 10 CONTINUE C C PERFORM THE CONVERT DIGIT FROM B-ARY TO B-ARY GRAY CODE C GTEMP(HISUM) = YTEMP(HISUM) DO 20 I = HISUM - 1,0,-1 TEM = YTEMP(I) - YTEMP(I+1) IF (TEM.LT.0) THEN GTEMP(I) = TEM + QS ELSE GTEMP(I) = MOD(TEM,QS) END IF 20 CONTINUE C C FINDING THE POSITION OF COLUMN THAT IT'S VALUE C HAS BEEN CHANGED. C DO 30 I = 0,HISUM IF (GTEMP(I).NE.PGTEMP(I)) THEN POS = I DIFF = GTEMP(I) - PGTEMP(I) END IF PGTEMP(I) = GTEMP(I) 30 CONTINUE C C UPDATE THE NEW VECTOR BY ADDING THE VALUE OF THE COLUMN THAT C HAS BEEN CHANGED. C DO 50 K = 1,S R = 0.0 DO 40 I = MAXL,0,-1 ZTEMP(K,I) = ZTEMP(K,I) + SCOEF(K,I,POS)*DIFF R = MOD(ZTEMP(K,I),QS) + RQS*R 40 CONTINUE QUASI(K) = R*RQS 50 CONTINUE C C UPDATE NEXTN AND, IF NEEDED, TESTN AND C HISUM C NEXTN = NEXTN + 1 IF (NEXTN.EQ.TESTN) THEN TESTN = TESTN*QS HISUM = HISUM + 1 C C SINCE FLAG(2) IS TRUE, C HISUM STAYS UNDER 20 C END IF C RETURN END C C THIS MODIFIED ALGORITHM 659, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 14, NO. 1, P.88. SUBROUTINE INFAUR(FLAG,DIMEN,ATMOST,QUASI,MAXL,IFLAG,MAXX) C C THIS SUBROUTINE FIRST CHECKS WHETHER C THE USER-SUPPLIED DIMENSION "DIMEN" OF THE C QUASIRANDOM VECTORS IS ACCEPTABLE C (STRICTLY BETWEEN 1 AND 41) : IF SO, C FLAG(1)=.TRUE. C C THEN IT CALCULATES AN UPPER SUMMATION C LIMIT "HISUM" BASED ON "DIMEN" AND THE C USER-SUPPLIED NUMBER "ATMOST" OF QUASIRANDOM C VECTORS REQUIRED. FLAG(2)=.TRUE. IF C ATMOST IS OK. C C IF FLAG(1) AND FLAG(2) ARE TRUE, C "INFAUR" NEXT PRODUCES THE OTHER C OUTPUTS LISTED BELOW PASSED TO C SUBROUTINE GOFAUR VIA LABELLED C COMMON "FAURE". THESE OUTPUTS ARE C IRRELEVANT TO THE USER. C C FIRST CALL INFAUR. IF FLAG(1) AND C FLAG(2) ARE TRUE, EACH (SUBSEQUENT) C CALL TO GSFAUR GENERATES A NEW C QUASIRANDOM VECTOR. C C INPUTS : DIMEN, ATMOST, MAXL C C OUTPUTS C TO USERS CALLING PROGRAM: C QUASI :FIRST SEQUENCE C FLAG C QSS : SAME AS QS - SEE BELOW C C C TO GSFAUR: C S :DIMENSION C QS :SMALLEST PRIME >=S C SCOEF :SCRAMBLING GENERATING MATRICES C NEXTN :THE NUMBER OF THE C NEXT QUASIRANDOM C VECTOR,INITIALIZED C TO TESTN-1 HERE. C TESTN :INITIALIZED TO 0 C HISUM :AFTER BEING USED TO C PRODUCE COEF, INITIALIZED C TO 3 FOR GOFAUR. C RQS :1.0/QS. C C C .. Parameters .. INTEGER MAXLN PARAMETER (MAXLN=30) C C C .. Scalar Arguments .. INTEGER ATMOST,DIMEN,IFLAG,MAXL,MAXX C .. C .. Local Arrays .. DOUBLE PRECISION QUASI(500) LOGICAL FLAG(2) C .. C .. Scalars in Common .. DOUBLE PRECISION RQS INTEGER HISUM,NEXTN,QS,S,TESTN C .. C .. Arrays in Common .. INTEGER COEF(1:500,0:MAXLN,0:MAXLN),PGTEMP(0:MAXLN), + SCOEF(1:500,0:MAXLN,0:MAXLN),ZTEMP(1:500,0:MAXLN) C .. C .. Local Scalars .. DOUBLE PRECISION R INTEGER I,IN,J,K,K1,K2,L,TEMS,TEMS2 C .. C .. Local Arrays .. INTEGER LSM(1:500,0:MAXLN,0:MAXLN),PRIMES(0:4,0:99), + SHIFT(1:500,0:MAXLN),TSCOEF(1:500,0:MAXLN,0:MAXLN), + USHIFT(0:MAXLN),USM(0:MAXLN,0:MAXLN) C .. C .. External Subroutines .. EXTERNAL GNCRML,GNCRMU C .. C .. Intrinsic Functions .. INTRINSIC LOG,MOD,NINT,REAL,INT C .. C .. Common blocks .. COMMON /FAURE/S,QS,HISUM,SCOEF,NEXTN,TESTN,RQS,PGTEMP,ZTEMP,COEF C .. C .. Save statement .. SAVE /FAURE/ C .. C .. Data statements .. C DATA (PRIMES(0,I),I=0,99)/0,1,2,3,5,5,7,7,11,11,11,11,13,13,17,17, + 17,17,19,19,23,23,23,23,29,29,29,29,29,29,31,31,37,37,37,37, + 37,37,41,41,41,41,43,43,47,47,47,47,53,53,53,53,53,53,59,59, + 59,59,59,59,61,61,67,67,67,67,67,67,71,71,71,71,73,73,79,79, + 79,79,79,79,83,83,83,83,89,89,89,89,89,89,97,97,97,97,97,97, + 97,97,101,101/ DATA (PRIMES(1,I),I=0,99)/101,101,103,103,107,107,107,107,109,109, + 113,113,113,113,127,127,127,127,127,127,127,127,127,127,127, + 127,127,127,131,131,131,131,137,137,137,137,137,137,139,139, + 149,149,149,149,149,149,149,149,149,149,151,151,157,157,157, + 157,157,157,163,163,163,163,163,163,167,167,167,167,173,173, + 173,173,173,173,179,179,179,179,179,179,181,181,191,191,191, + 191,191,191,191,191,191,191,193,193,197,197,197,197,199,199/ DATA (PRIMES(2,I),I=0,99)/211,211,211,211,211,211,211,211,211,211, + 211,211,223,223,223,223,223,223,223,223,223,223,223,223,227, + 227,227,227,229,229,233,233,233,233,239,239,239,239,239,239, + 241,241,251,251,251,251,251,251,251,251,251,251,257,257,257, + 257,257,257,263,263,263,263,263,263,269,269,269,269,269,269, + 271,271,277,277,277,277,277,277,281,281,281,281,283,283,293, + 293,293,293,293,293,293,293,293,293,307,307,307,307,307,307/ DATA (PRIMES(3,I),I=0,99)/307,307,307,307,307,307,307,307,311,311, + 311,311,313,313,317,317,317,317,331,331,331,331,331,331,331, + 331,331,331,331,331,331,331,337,337,337,337,337,337,347,347, + 347,347,347,347,347,347,347,347,349,349,353,353,353,353,359, + 359,359,359,359,359,367,367,367,367,367,367,367,367,373,373, + 373,373,373,373,379,379,379,379,379,379,383,383,383,383,389, + 389,389,389,389,389,397,397,397,397,397,397,397,397,401,401/ DATA (PRIMES(4,I),I=0,99)/401,401,409,409,409,409,409,409,409,409, + 419,419,419,419,419,419,419,419,419,419,421,421,431,431,431, + 431,431,431,431,431,431,431,433,433,439,439,439,439,439,439, + 443,443,443,443,449,449,449,449,449,449,457,457,457,457,457, + 457,457,457,461,461,461,461,463,463,467,467,467,467,479,479, + 479,479,479,479,479,479,479,479,479,479,487,487,487,487,487, + 487,487,487,491,491,491,491,499,499,499,499,499,499,499,499/ C C CHECK S C S = DIMEN FLAG(1) = S .GT. 1 .AND. S .LT. 41 IF (.NOT.FLAG(1)) RETURN C TESTN = 0 K1 = INT(S/100) K2 = (S-K1*100) QS = PRIMES(K1,K2) C C COMPUTE LOG(ATMOST+TESTN) IN BASE QS C USING A RATIO OF NATURAL LOGS TO GET C AN UPPER BOUND ON (THE NUMBER OF C DIGITS IN THE BASE QS REPRESENTATION C OF ATMOST+TESTN) MINUS ONE. C HISUM = NINT(LOG(REAL(ATMOST+TESTN))/LOG(REAL(QS))) FLAG(2) = HISUM .LT. 30 IF (.NOT.FLAG(2)) RETURN C C NOW FIND BINOMIAL COEFFICIENTS MOD QS C IN A UPPER-TRIANGULAR MATRIX "COEF" C USING RECURSION BINOM(I,J)=BINOM(I,J-1) C +BINOM(I-1,J-1) AND A=B+C IMPLIES MOD(A,D)= C MOD(MOD(B,D)+MOD(C,D),D) C C CONSTRUCTING UPPER-TIRANGULAR GENERATOR MATRICES C UP TO DIMENSION S DO 20 I = 0,HISUM DO 10 J = 0,HISUM IF (I.EQ.J) THEN COEF(1,I,J) = 1 ELSE COEF(1,I,J) = 0 END IF 10 CONTINUE 20 CONTINUE COEF(2,0,0) = 1 DO 30 J = 1,HISUM COEF(2,0,J) = 1 COEF(2,J,J) = 1 30 CONTINUE DO 50 I = 1,HISUM DO 40 J = I + 1,HISUM COEF(2,I,J) = MOD(COEF(2,I,J-1)+COEF(2,I-1,J-1),QS) 40 CONTINUE 50 CONTINUE DO 90 K = 3,S DO 80 J = 0,HISUM DO 70 I = 0,HISUM COEF(K,I,J) = 0 IF (J.LT.I) THEN GO TO 70 ELSE TEMS = 0 END IF DO 60 L = 0,J TEMS = TEMS + COEF(K-1,I,L)*COEF(2,L,J) 60 CONTINUE COEF(K,I,J) = MOD(TEMS,QS) 70 CONTINUE 80 CONTINUE 90 CONTINUE C C GENERATE USER CHOICE OF GENERATOR MATRICES C IF (IFLAG.EQ.0) THEN MAXX = HISUM DO 120 K = 1,S DO 110 J = 0,HISUM DO 100 I = 0,HISUM SCOEF(K,I,J) = COEF(K,I,J) 100 CONTINUE SHIFT(K,J) = 0 110 CONTINUE 120 CONTINUE ELSE IF ((IFLAG.EQ.1) .OR. (IFLAG.EQ.3)) THEN MAXX = MAXL CALL GNCRML(LSM,SHIFT,MAXL) DO 160 K = 1,S DO 150 J = 0,HISUM DO 140 I = 0,MAXL IF (K.EQ.1) THEN SCOEF(K,I,J) = LSM(K,I,J) IF (IFLAG.EQ.3) THEN TSCOEF(K,I,J) = SCOEF(K,I,J) END IF ELSE TEMS = 0 IF (I.GE.HISUM) THEN IN = HISUM ELSE IN = J END IF DO 130 L = 0,IN TEMS = TEMS + LSM(K,I,L)*COEF(K,L,J) 130 CONTINUE SCOEF(K,I,J) = MOD(TEMS,QS) IF (IFLAG.EQ.3) THEN TSCOEF(K,I,J) = SCOEF(K,I,J) END IF END IF 140 CONTINUE 150 CONTINUE 160 CONTINUE END IF IF ((IFLAG.EQ.2) .OR. (IFLAG.EQ.3)) THEN CALL GNCRMU(USM,USHIFT,HISUM) IF (IFLAG.EQ.2) THEN MAXX = HISUM END IF DO 200 K = 1,S DO 190 J = 0,HISUM DO 180 I = 0,MAXX TEMS = 0 TEMS2 = 0 IN = HISUM DO 170 L = 0,IN IF (IFLAG.EQ.2) THEN TEMS = TEMS + COEF(K,I,L)*USM(L,J) IF (J.EQ.0) THEN TEMS2 = TEMS2 + + COEF(K,I,L)*USHIFT(L) END IF END IF IF (IFLAG.EQ.3) THEN TEMS = TEMS + TSCOEF(K,I,L)*USM(L,J) IF (J.EQ.0) THEN TEMS2 = TEMS2 + + TSCOEF(K,I,L)*USHIFT(L) END IF END IF 170 CONTINUE SCOEF(K,I,J) = MOD(TEMS,QS) IF (J.EQ.0) THEN TEMS2 = MOD(TEMS2,QS) IF (IFLAG.EQ.3) THEN SHIFT(K,I) = MOD((TEMS2+SHIFT(K,I)), + QS) ELSE SHIFT(K,I) = TEMS2 END IF END IF 180 CONTINUE 190 CONTINUE 200 CONTINUE END IF END IF C CALCULATING THESE COEFFICIENTS C MOD QS AVOIDS POSSIBLE OVERFLOW C PROBLEMS WITH RAW BINOMIAL COEFFICIENTS C RQS = 1.0/REAL(QS) C GENERATE FIRST SCRAMBLED QUASI RANDOM VECTOR DO 220 K = 1,S R = 0.0 DO 210 I = MAXX,0,-1 ZTEMP(K,I) = SHIFT(K,I) R = ZTEMP(K,I) + RQS*R 210 CONTINUE QUASI(K) = R*RQS 220 CONTINUE C C SET UP FIRST VECTOR AND VALUES FOR GSFAUR C DO 230 I = 0,HISUM PGTEMP(I) = 0 230 CONTINUE TESTN = QS NEXTN = 1 HISUM = 0 C NOW COMPLETE INITIALIZATION C AS DESCRIBED IN SECTION 4. C NEXTN HAS 1 DIGITS IN BASE C QS, SO HISUM EQUALS 0. C RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Snied' then mkdir 'Snied' fi cd 'Snied' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' include makefile.inc all: Res1 SRCLIBS= $(LINPACK) $(EISPACK_DP) $(BLAS) $(PORT) src.o: $(INCSRC)/src.f $(F77) $(F77OPTS) -c $(INCSRC)/src.f extra.o: $(INCSRC)/extra.f $(F77) $(F77OPTS) -c $(INCSRC)/extra.f DRIVERS= driver RESULTS= Res1 SRCLIBS = $(TIMER) Objs1= driver.o src.o extra.o driver: $(Objs1) $(F77LINK) $(F77LINKOPTS) -o driver $(Objs1) $(SRCLIBS) Res1: driver ./driver >Res1 diffres: Res1 res1 echo "Differences in results from driver" $(DIFF) Res1 res1 clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << "SHAR_EOF" > 'data1' 683644482 135580830 114750887 803113694 903684270 72157854 59862132 505680715 779216195 193586105 329775408 1049470377 790869458 529177858 770596904 1072006936 140712293 266923202 844524433 114947652 963169856 269847095 96724448 724232372 729133779 872385703 803906910 638059470 526484481 843345462 431170622 1057600395 413165996 637333254 738214251 987594285 394318066 718140016 1064908418 636856468 856354958 404426257 6726334 916596837 284819911 905844483 331078851 514311810 280380871 612727606 184724875 578053121 197887820 50641495 423123265 227871988 36070792 465200006 187446371 99681590 983367218 414604440 243776932 990984772 289169669 72196638 522106132 18453143 1008259046 521211985 88531465 816692680 1001849242 432014653 124800720 554138134 166704786 634462798 1019392879 691725817 459507961 1043416360 536793729 1046579554 403993816 412804346 138344837 24248109 184604157 535890993 869278100 125797533 322976214 161488107 755812753 541238672 841393456 909350792 521628897 567331204 118313388 265208159 468752989 169896169 959104945 500814545 434091372 113024960 444725900 426474642 263586462 432120698 45230250 524531296 327681208 288240183 353741403 192041686 42986304 258035065 388012188 537416171 292523858 794840680 502322834 423539584 884204171 419181022 618098449 257684411 288733678 145066288 756812965 357038497 881037192 314634603 182012383 753118967 636654128 318546833 71382991 400554255 153755763 187747483 306343731 1044780 402486338 564947470 934172722 813806089 705671218 795868370 926341726 51211860 739601563 176287958 1016675820 950546833 267929114 275671035 666917527 431607592 177069271 570312219 125215260 886835363 75079442 767034988 879903592 525306349 560816857 973795128 774079162 143751759 137446578 209360240 1012892947 398130103 838233295 47600869 312065869 24489594 686091463 1064468817 610953508 593516621 736434891 331884657 295796744 991620199 197203105 694133440 834981040 813915471 1043396296 246522553 741863988 188767345 707349375 236730174 943797129 1073073438 917394477 81723820 598848363 17981413 322240492 766981617 282924362 643850936 914020733 586984616 128906358 739978294 543175409 351896331 1017052692 802102481 946877967 615436223 305014361 1038674293 372513724 902652345 552583238 507262995 36761397 753350193 116032862 328670174 103551339 1045088961 667010209 737933140 945361853 324875181 353779021 866388681 255506551 65493883 220170356 837239777 771797016 542598812 568193573 359763920 870736829 304824104 430870768 116534859 741311666 157845915 861348937 985097428 714048271 306167010 709889598 498574308 382576036 360956609 1005272948 680195561 845547054 1026217890 374002007 866217449 695335995 462692419 314161108 1048195907 709323042 604523286 581285813 463752482 919585493 518104348 905143615 884163267 99734344 332656669 398822849 973544107 459181032 1001255433 846689598 32962550 6372759 727573859 495342420 738892195 834835824 1047729137 529022179 1032177776 185764687 483299066 818974724 826550698 987239779 684232759 940422073 224992310 164068454 827234589 690302353 32197089 628566875 95166490 866947171 476745989 604118539 121738255 667653407 987681310 749263282 863528381 341372968 244006930 951984484 614173490 685305366 592527405 383027710 952965461 21271273 653806680 745437736 180332506 331728868 207307411 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' DOUBLE PRECISION FUNCTION TESTF(N,DIMEN,QUASI) C C This version : 4 Mar 1992 C C Provides a variety of test integrals for quasi-random C sequences. A call on TESTF computes an estimate of the C integral ; a call on EXACTF computes the exact value. C C .. Scalar Arguments .. INTEGER DIMEN,N C .. C .. Array Arguments .. DOUBLE PRECISION QUASI(*) C .. C .. Local Scalars .. DOUBLE PRECISION X INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,FLOAT,MOD,SIN C .. C .. Entry Points .. DOUBLE PRECISION EXACTF C .. GO TO (10,40,80,150) N C ENTRY EXACTF(N,DIMEN) GO TO (30,60,140,170) N C C Test integral 1 C 10 TESTF = 1.0 DO 20 I = 1,DIMEN TESTF = TESTF*ABS(4*QUASI(I)-2) 20 CONTINUE RETURN C 30 EXACTF = 1.0 RETURN C C Test integral 2 C 40 TESTF = 1.0 DO 50 I = 1,DIMEN TESTF = TESTF*I*COS(I*QUASI(I)) 50 CONTINUE RETURN C 60 EXACTF = 1.0 DO 70 I = 1,DIMEN EXACTF = EXACTF*SIN(FLOAT(I)) 70 CONTINUE RETURN C C Test integral 3 C 80 TESTF = 1.0 DO 130 I = 1,DIMEN X = 2*QUASI(I) - 1 GO TO (90,100,110,120) MOD(I,4) 90 TESTF = TESTF*X GO TO 130 100 TESTF = TESTF* (2*X*X-1) GO TO 130 110 TESTF = TESTF* (4*X*X-3)*X GO TO 130 120 X = X*X TESTF = TESTF* (8*X*X-8*X+1) 130 CONTINUE RETURN C 140 EXACTF = 0.0 RETURN C C Test integral 4 C 150 TESTF = 0 X = 1 DO 160 I = 1,DIMEN X = -X*QUASI(I) TESTF = TESTF + X 160 CONTINUE RETURN C C 170 X = 1.0/ (2** (DIMEN)) IF (MOD(DIMEN,2).EQ.0) THEN EXACTF = (X-1)/3 ELSE EXACTF = (X+1)/3 END IF RETURN C END PROGRAM TSNIED C C C THIS PROGRAM TESTS ACCURACY OF C NUMERICAL INTEGRATION USING "GNNIED" C AND INTEGRAND (2) OF DAVIS AND C RABINOWITZ, PAGE 406 C C User Define: C DIMEN : dimension C ATMOST : sequence length C SAMS : Number of replications C IFLAG: User Choice of type Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C C .. Parameters .. INTEGER MAXDIM PARAMETER (MAXDIM=318) C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,II,IFLAG,K,J,SAM C .. C .. Local Arrays .. DOUBLE PRECISION OUTS(MAXDIM,10000) REAL SECOND, T1 C .. C .. C .. External Subroutines .. EXTERNAL GNNIED,TIME C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. SAM = 1 DIMEN = 2 IFLAG = 1 ATMOST = 2**12 DO 30 II = 1,SAM WRITE (*,FMT=*) 'I = ITERATION NUMBER' WRITE (*,FMT=*) 'EI = ESTIMATED INTEGRAL' T1 = SECOND() CALL GNNIED(DIMEN,ATMOST,IFLAG,OUTS) SUM = 0.0 K = 9 DO 20 I = 1,ATMOST F = 1.0 DO 10 J = 1,DIMEN F = F*ABS(4.0*OUTS(J,I)-2.0) 10 CONTINUE IF (MOD(I,2**K).EQ.0) THEN K = K + 1 WRITE (*,FMT=*) 'I = ',I WRITE (*,FMT=*) 'EI = ',SUM/I END IF SUM = SUM + F 20 CONTINUE T1 = SECOND() - T1 WRITE (*,FMT=*) 'Total time elapsed = ',T1 30 CONTINUE STOP END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' I = ITERATION NUMBER EI = ESTIMATED INTEGRAL I = 512 EI = 0.9999562489621427 I = 1024 EI = 0.999713791146755 I = 2048 EI = 0.9999491497757775 I = 4096 EI = 0.9999468969847435 Total time elapsed = 0.13607088 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'extra.f' then echo shar: will not over-write existing file "'extra.f'" else cat << "SHAR_EOF" > 'extra.f' INTEGER FUNCTION EXOR(IIN,JIN) C C THIS FUNCTION CALCULATES THE EXCLUSIVE-OR OF ITS C TWO INPUT PARAMETERS C C .. Scalar Arguments .. INTEGER IIN,JIN C .. C .. Local Scalars .. INTEGER I,J,K,L C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. I = IIN J = JIN K = 0 L = 1 C 10 IF (I.EQ.J) THEN EXOR = K RETURN END IF C C CHECK THE CURRENT RIGHT-HAND BITS OF I AND J. C IF THEY DIFFER, SET THE APPROPRIATE BIT OF K. C IF (MOD(I,2).NE.MOD(J,2)) K = K + L I = I/2 J = J/2 L = 2*L GO TO 10 END DOUBLE PRECISION FUNCTION UNI() * * Random number generator, adapted from F. James * "A Review of Random Number Generators" * Comp. Phys. Comm. 60(1990), pp. 329-344. * C .. Parameters .. DOUBLE PRECISION TWOM24 PARAMETER (TWOM24=1D0/16777216.0) C .. C .. Local Scalars .. DOUBLE PRECISION CARRY INTEGER I,J C .. C .. Local Arrays .. DOUBLE PRECISION SEEDS(24) C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Save statement .. SAVE I,J,CARRY,SEEDS C .. C .. Data statements .. DATA I,J,CARRY/24,10,0.0d0/ DATA SEEDS/0.8804418d0,0.2694365d0,0.0367681d0,0.4068699d0, +0.4554052d0, + 0.2880635d0,0.1463408d0,0.2390333d0,0.6407298d0, +0.1755283d0,0.7132940d0, + 0.4913043d0,0.2979918d0,0.1396858d0,0.3589528d0, +0.5254809d0,0.9857749d0, + 0.4612127d0,0.2196441d0,0.7848351d0,0.4096100d0, +0.9807353d0,0.2689915d0, + 0.5140357d0/ C .. UNI = SEEDS(I) - SEEDS(J) - CARRY IF (UNI.LT.0) THEN UNI = UNI + 1 CARRY = TWOM24 ELSE CARRY = 0 END IF SEEDS(I) = UNI I = 24 - MOD(25-I,24) J = 24 - MOD(25-J,24) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE CALCV(PX,B,V,MAXV) C C This version : 12 February 1991 C C See the general comments on implementing Niederreiter's C low-discrepancy sequences. C C This program calculates the values of the constants V(J,R) as C described in BFN section 3.3. It is called from either CALCC or C SCLCC2. The values transmitted through common /FIELD/ determine C which field we are working in. C C INPUT : C PX is the appropriate irreducible polynomial for the dimension C currently being considered. The degree of PX will be called E. C B is the polynomial defined in section 2.3 of BFN. On entry to C the subroutine, the degree of B implicitly defines the parameter C J of section 3.3, by degree(B) = E*(J-1). C MAXV gives the dimension of the array V. C On entry, we suppose that the common block /FIELD/ has been set C up correctly (using SETFLD). C C OUTPUT : C On output B has been multiplied by PX, so its degree is now E*J. C V contains the values required. C C USES : C The subroutine PLYMUL is used to multiply polynomials. C C C ------------------------------------------------------------ C C The following COMMON block, used by many subroutines, C gives the order Q of a field, its characteristic P, and its C addition, multiplication, and subtraction tables. C The parameter MAXQ gives the order of the largest field to C be handled. C C C The following definitions concern the representation of C polynomials. C C C The parameter MAXDEG gives the highest degree of polynomial C to be handled. Polynomials stored as arrays have the C coefficient of degree n in POLY(N), and the degree of the C polynomial in POLY(-1). The parameter DEG is just to remind C us of this last fact. A polynomial which is identically 0 C is given degree -1. C C A polynomial can also be stored in an integer I, with C I = AN*Q**N + ... + A0. C Routines ITOP and PTOI convert between these two formats. C C ----------------------------------------------------------- C C C C C .. Parameters .. INTEGER MAXQ PARAMETER (MAXQ=50) INTEGER MAXDEG,DEG PARAMETER (MAXDEG=50,DEG=-1) C .. C .. Scalar Arguments .. INTEGER MAXV C .. C .. Array Arguments .. INTEGER B(-1:MAXDEG),PX(-1:MAXDEG),V(0:MAXV) C .. C .. Scalars in Common .. INTEGER P,Q C .. C .. Arrays in Common .. INTEGER ADD(0:MAXQ-1,0:MAXQ-1),MUL(0:MAXQ-1,0:MAXQ-1), + SUB(0:MAXQ-1,0:MAXQ-1) C .. C .. Local Scalars .. INTEGER BIGM,DUMMY,I,KJ,M,NONZER,R,TERM C .. C .. Local Arrays .. INTEGER H(-1:MAXDEG) C .. C .. External Subroutines .. EXTERNAL PLYMUL C .. C .. Common blocks .. COMMON /FIELD/P,Q,ADD,MUL,SUB C .. C .. Statement Functions .. INTEGER ARBIT C .. C .. Save statement .. SAVE /FIELD/ C .. C .. Statement Function definitions .. ARBIT(DUMMY) = 1 C .. C C We use ARBIT() to indicate where the user can place C an arbitrary element of the field of order Q, while NONZER C shows where he must put an arbitrary non-zero element C of the same field. For the code, C this means 0 <= ARBIT < Q and 0 < NONZER < Q. Within these C limits, the user can do what he likes. ARBIT is declared as C a function as a reminder that a different arbitrary value may C be returned each time ARBIT is referenced. C C BIGM is the M used in section 3.3. C It differs from the [little] m used in section 2.3, C denoted here by M. C NONZER = 1 C C E = PX(DEG) C C The poly H is PX**(J-1), which is the value of B on arrival. C In section 3.3, the values of Hi are defined with a minus sign : C don't forget this if you use them later ! C DO 10 I = -1,B(DEG) H(I) = B(I) 10 CONTINUE BIGM = H(DEG) C C Now multiply B by PX so B becomes PX**J. C In section 2.3, the values of Bi are defined with a minus sign : C don't forget this if you use them later ! C CALL PLYMUL(PX,B,B) M = B(DEG) C C We don't use J explicitly anywhere, but here it is just in case. C C J = M/E C C Now choose a value of Kj as defined in section 3.3. C We must have 0 <= Kj < E*J = M. C The limit condition on Kj does not seem very relevant C in this program. C KJ = BIGM C C Now choose values of V in accordance with the conditions in C section 3.3 C DO 20 R = 0,KJ - 1 V(R) = 0 20 CONTINUE V(KJ) = 1 C IF (KJ.LT.BIGM) THEN C TERM = SUB(0,H(KJ)) C DO 30 R = KJ + 1,BIGM - 1 V(R) = ARBIT(DUMMY) C C Check the condition of section 3.3, C remembering that the H's have the opposite sign. C TERM = SUB(TERM,MUL(H(R),V(R))) 30 CONTINUE C C Now V(BIGM) is anything but TERM C V(BIGM) = ADD(NONZER,TERM) C DO 40 R = BIGM + 1,M - 1 V(R) = ARBIT(DUMMY) 40 CONTINUE C ELSE C This is the case KJ .GE. BIGM C DO 50 R = KJ + 1,M - 1 V(R) = ARBIT(DUMMY) 50 CONTINUE C END IF C C Calculate the remaining V's using the recursion of section 2.3, C remembering that the B's have the opposite sign. C DO 70 R = 0,MAXV - M TERM = 0 DO 60 I = 0,M - 1 TERM = SUB(TERM,MUL(B(I),V(R+I))) 60 CONTINUE V(R+M) = TERM 70 CONTINUE C RETURN END C C ***** end of SUBROUTINE CALCV INTEGER FUNCTION CHARAC(QIN) C C This version : 12 December 1991 C C This function gives the characteristic for a field of C order QIN. If no such field exists, or if QIN is out of C the range we can handle, returns 0. C C C ------------------------------------------------------------ C C The following COMMON block, used by many subroutines, C gives the order Q of a field, its characteristic P, and its C addition, multiplication and subtraction tables. C The parameter MAXQ gives the order of the largest field to C be handled. C C C The following definitions concern the representation of C polynomials. C C C The parameter MAXDEG gives the highest degree of polynomial C to be handled. Polynomials stored as arrays have the C coefficient of degree n in POLY(N), and the degree of the C polynomial in POLY(-1). The parameter DEG is just to remind C us of this last fact. A polynomial which is identically 0 C is given degree -1. C C A polynomial can also be stored in an integer I, with C I = AN*Q**N + ... + A0. C Routines ITOP and PTOI convert between these two formats. C C ----------------------------------------------------------- C C C C .. Parameters .. INTEGER MAXQ PARAMETER (MAXQ=50) C .. C .. Scalar Arguments .. INTEGER QIN C .. C .. Local Arrays .. INTEGER CH(MAXQ) C .. C .. Save statement .. SAVE CH C .. C .. Data statements .. C DATA CH/0,2,3,2,5,0,7,2,3,0,11,0,13,0,0,2,17,0,19,0,0,0,23,0,5,0, + 3,0,29,0,31,2,0,0,0,0,37,0,0,0,41,0,43,0,0,0,47,0,7,0/ C .. C IF (QIN.LE.1 .OR. QIN.GT.MAXQ) THEN CHARAC = 0 ELSE CHARAC = CH(QIN) END IF C END C ***** end of SCLCC2 SUBROUTINE GNCRML(MAXDIM,NBITS,DIMEN,LSM,SHIFT) C GENERATING LOWER TRIAGULAR SCRMABLING MATRICES C AND SHIFT VECTORS. C INPUTS : C FROM SCLCC2 : MAXDIM,NBITS,DIMEN C C OUTPUTS : C TO SCLCC2 : LSM, SHIFT C .. Scalar Arguments .. INTEGER DIMEN,MAXDIM,NBITS C .. C .. Array Arguments .. INTEGER LSM(MAXDIM,0:NBITS),SHIFT(MAXDIM) C .. C .. Local Scalars .. INTEGER I,J,L,LL,P,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. DO 30 P = 1,DIMEN SHIFT(P) = 0 L = 1 DO 20 I = NBITS - 1,0,-1 LSM(P,I) = 0 STEMP = MOD((INT(UNI()*1000.0)),2) SHIFT(P) = SHIFT(P) + STEMP*L L = 2*L LL = 1 DO 10 J = NBITS - 1,0,-1 IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.LT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF LSM(P,I) = LSM(P,I) + TEMP*LL LL = 2*LL 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GNCRMU(NBITS,USM,USHIFT) C GENERATING UPPER TRIAGULAR SCRMABLING MATRIX C AND SHIFT VECTOR. C INPUTS : C FROM SCLCC2: NBITS C C OUTPUTS : C TO SCLCC2 : USM, USHIFT C .. Scalar Arguments .. INTEGER NBITS C .. C .. Array Arguments .. INTEGER USHIFT(0:NBITS-1),USM(0:NBITS-1,0:NBITS-1) C .. C .. Local Scalars .. INTEGER I,J,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. DO 20 I = 0,NBITS - 1 STEMP = MOD((INT(UNI()*1000.0)),2) USHIFT(I) = STEMP DO 10 J = 0,NBITS - 1 IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.GT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF USM(I,J) = TEMP 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE GNNIED(DIMEN,ATMOST,IFLAG,OUTS) C C This is Modified Program of Niederreiter Sequences Generator C It generates Various Scrambled Niederreiter Sequences. C SNIED and its associated subroutines are listed below. C C SNIED C SINLO2 C SGOLO2 C SCLCC2 C CALCV C CHARAC C SETFLD C PLYMUL C TESTF C C The suffix 2 indicates routines for use only by C the set of programs tailored for base 2. C C User Define: C DIMEN : dimension C ATMOST : sequence length C IFLAG: User Choice of Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C C C C .. Parameters .. INTEGER MAXDIM PARAMETER (MAXDIM=318) C .. C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,IFLAG,J,SKIP C .. C .. Local Arrays .. DOUBLE PRECISION QUASI(MAXDIM),OUTS(MAXDIM,10000) REAL TARRAY(2) C .. C .. External Subroutines .. EXTERNAL SGOLO2,SINLO2 C .. SKIP = 0 CALL SINLO2(DIMEN,SKIP,IFLAG) DO 20 I = 1,ATMOST CALL SGOLO2(QUASI) DO 10 J = 1,DIMEN OUTS(J,I) = QUASI(J) 10 CONTINUE 20 CONTINUE RETURN END C C ***** end of SUBROUTINE SINLO2 SUBROUTINE PLYMUL(PA,PB,PC) C C This version : 12 December 1991 C C C ------------------------------------------------------------ C C The following COMMON block, used by many subroutines, C gives the order Q of a field, its characteristic P, and its C addition, multiplication and subtraction tables. C The parameter MAXQ gives the order of the largest field to C be handled. C C C The following definitions concern the representation of C polynomials. C C C The parameter MAXDEG gives the highest degree of polynomial C to be handled. Polynomials stored as arrays have the C coefficient of degree n in POLY(N), and the degree of the C polynomial in POLY(-1). The parameter DEG is just to remind C us of this last fact. A polynomial which is identically 0 C is given degree -1. C C A polynomial can also be stored in an integer I, with C I = AN*Q**N + ... + A0. C Routines ITOP and PTOI convert between these two formats. C C ----------------------------------------------------------- C C C C C Multiplies polynomial PA by PB putting the result in PC. C Coefficients are elements of the field of order Q. C C .. Parameters .. INTEGER MAXQ PARAMETER (MAXQ=50) INTEGER MAXDEG,DEG PARAMETER (MAXDEG=50,DEG=-1) C .. C .. Array Arguments .. INTEGER PA(-1:MAXDEG),PB(-1:MAXDEG),PC(-1:MAXDEG) C .. C .. Scalars in Common .. INTEGER P,Q C .. C .. Arrays in Common .. INTEGER ADD(0:MAXQ-1,0:MAXQ-1),MUL(0:MAXQ-1,0:MAXQ-1), + SUB(0:MAXQ-1,0:MAXQ-1) C .. C .. Local Scalars .. INTEGER DEGA,DEGB,DEGC,I,J,TERM C .. C .. Local Arrays .. INTEGER PT(-1:MAXDEG) C .. C .. Intrinsic Functions .. INTRINSIC MAX,MIN C .. C .. Common blocks .. COMMON /FIELD/P,Q,ADD,MUL,SUB C .. C .. Save statement .. SAVE /FIELD/ C .. DEGA = PA(DEG) DEGB = PB(DEG) IF (DEGA.EQ.-1 .OR. DEGB.EQ.-1) THEN DEGC = -1 ELSE DEGC = DEGA + DEGB END IF IF (DEGC.GT.MAXDEG) THEN WRITE (*,FMT=*) ' PLYMUL : Degree of product exceeds MAXDEG' STOP END IF C DO 20 I = 0,DEGC TERM = 0 DO 10 J = MAX(0,I-DEGA),MIN(DEGB,I) TERM = ADD(TERM,MUL(PA(I-J),PB(J))) 10 CONTINUE PT(I) = TERM 20 CONTINUE C PC(DEG) = DEGC DO 30 I = 0,DEGC PC(I) = PT(I) 30 CONTINUE DO 40 I = DEGC + 1,MAXDEG PC(I) = 0 40 CONTINUE RETURN END C SUBROUTINE SCLCC2(IFLAG) C C This is modified routine of SCLCC2. C C This program calculates the values of the constants C(I,J,R). C As far as possible, we use Niederreiter's notation. C For each value of I, we first calculate all the corresponding C values of C : these are held in the array CI. All these C values are either 0 or 1. Next we pack the values into the C array CJ, in such a way that CJ(I,R) holds the values of C C for the indicated values of I and R and for every value of C J from 1 to NBITS. The most significant bit of CJ(I,R) C (not counting the sign bit) is C(I,1,R) and the least C significant bit is C(I,NBITS,R). C When all the values of CJ have been calculated, we return C this array to the calling program. C C -------------------------------------------------------------- C C We define the common block /COMM2/ and some C associated parameters. These are for use in the base 2 C version of the generator. C C C The parameter MAXDIM is the maximum dimension that will be used. C NBITS is the number of bits (not counting the sign) in a C fixed-point integer. C C C The common block /COMM2/ : C CJ - Contains the packed values of Niederreiter's C(I,J,R) C SCJ - Contains the packed values of User choice of scrambled C Niederreiter's C(I,J,R) C DIMEN - The dimension of the sequence to be generated C COUNT - The index of the current item in the sequence, C expressed as an array of bits. COUNT(R) is the same C as Niederreiter's AR(N) (page 54) except that N is C implicit. C NEXTQ - The numerators of the next item in the series. These C are like Niederreiter's XI(N) (page 54) except that C N is implicit, and the NEXTQ are integers. To obtain C the values of XI(N), multiply by RECIP (see GOLO2). C C Array SCJ of the common block is set up by subroutine SCLCC2. C The other items in the common block are set up by INLO2. C C -------------------------------------------------------------- C C The following COMMON block, used by many subroutines, C gives the order Q of a field, its characteristic P, and its C addition, multiplication and subtraction tables. C The parameter MAXQ gives the order of the largest field to C be handled. C C C The following definitions concern the representation of C polynomials. C C C The parameter MAXDEG gives the highest degree of polynomial C to be handled. Polynomials stored as arrays have the C coefficient of degree n in POLY(N), and the degree of the C polynomial in POLY(-1). The parameter DEG is just to remind C us of this last fact. A polynomial which is identically 0 C is given degree -1. C C A polynomial can also be stored in an integer I, with C I = AN*Q**N + ... + A0. C Routines ITOP and PTOI convert between these two formats. C C --------------------------------------------------------------- C C C C MAXE - We need MAXDIM irreducible polynomials over Z2. C MAXE is the highest degree among these. C MAXV - The maximum possible index used in V. C C C INPUT : C The array SCJ to be initialised, and DIMEN the number of C dimensions we are using, are transmitted through /COMM2/. C C OUTPUT : C The array SCJ is returned to the calling program. C C USES : C Subroutine SETFLD is used to set up field arithmetic tables. C (Although this is a little heavy-handed for the field of C order 2, it makes for uniformity with the general program.) C Subroutine CALCV is used to for the auxiliary calculation C of the values V needed to get the Cs. C C C Non-Standard Intrinsic Funtion for f77 C But Standard Intrinsic Fuction for f90 IBITS IS USED. C C .. Parameters .. INTEGER MAXDIM,NBITS PARAMETER (MAXDIM=318,NBITS=31) INTEGER MAXQ PARAMETER (MAXQ=50) INTEGER MAXDEG,DEG PARAMETER (MAXDEG=50,DEG=-1) INTEGER MAXE,MAXV PARAMETER (MAXE=11,MAXV=NBITS+MAXE) C .. C .. Scalar Arguments .. INTEGER IFLAG C .. C .. Scalars in Common .. INTEGER COUNT,DIMEN,P,Q C .. C .. Arrays in Common .. INTEGER ADD(0:MAXQ-1,0:MAXQ-1),MUL(0:MAXQ-1,0:MAXQ-1), + NEXTQ(MAXDIM),SCJ(MAXDIM,0:NBITS-1),SUB(0:MAXQ-1,0:MAXQ-1) C .. C .. Local Scalars .. INTEGER E,I,J,K,L,PP,R,TEMP1,TEMP2,TEMP3,TEMP4,TERM,U C .. C .. Local Arrays .. INTEGER B(-1:MAXDEG),CI(NBITS,0:NBITS-1),CJ(MAXDIM,0:NBITS-1), + IRRED(MAXDIM,-1:MAXE),LSM(MAXDIM,0:NBITS-1),PX(-1:MAXDEG), + SHIFT(MAXDIM),TV(MAXDIM,0:NBITS-1,0:NBITS-1), + USHIFT(0:NBITS-1),USM(0:NBITS-1,0:NBITS-1),V(0:MAXV) C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. External Subroutines .. EXTERNAL CALCV,GNCRML,GNCRMU,SETFLD C .. C .. Intrinsic Functions .. INTRINSIC IBITS,MOD C .. C .. Common blocks .. COMMON /COMM2/SCJ,DIMEN,COUNT,NEXTQ COMMON /FIELD/P,Q,ADD,MUL,SUB C .. C .. Save statement .. SAVE /COMM2/,/FIELD/,IRRED C .. C .. Data statements .. C C This DATA statement supplies the coefficients and the C degrees of the first 12 irreducible polynomials over Z2. C They are taken from Lidl and Niederreiter, FINITE FIELDS, C Cambridge University Press (1984), page 553. C The order of these polynomials is the same as the order in C file 'irrtabs.dat' used by the general program. C C In this block PX(I, -1) is the degree of the Ith polynomial, C and PX(I, N) is the coefficient of x**n in the Ith polynomial. C DATA (IRRED(1,I),I=-1,1)/1,0,1/ DATA (IRRED(2,I),I=-1,1)/1,1,1/ DATA (IRRED(3,I),I=-1,2)/2,1,1,1/ DATA (IRRED(4,I),I=-1,3)/3,1,1,0,1/ DATA (IRRED(5,I),I=-1,3)/3,1,0,1,1/ DATA (IRRED(6,I),I=-1,4)/4,1,1,0,0,1/ DATA (IRRED(7,I),I=-1,4)/4,1,0,0,1,1/ DATA (IRRED(8,I),I=-1,4)/4,1,1,1,1,1/ DATA (IRRED(9,I),I=-1,5)/5,1,0,1,0,0,1/ DATA (IRRED(10,I),I=-1,5)/5,1,0,0,1,0,1/ DATA (IRRED(11,I),I=-1,5)/5,1,1,1,1,0,1/ DATA (IRRED(12,I),I=-1,5)/5,1,1,1,0,1,1/ DATA (IRRED(13,I),I=-1,5)/5,1,1,0,1,1,1/ DATA (IRRED(14,I),I=-1,5)/5,1,0,1,1,1,1/ DATA (IRRED(15,I),I=-1,6)/6,1,1,0,0,0,0,1/ DATA (IRRED(16,I),I=-1,6)/6,1,0,0,1,0,0,1/ DATA (IRRED(17,I),I=-1,6)/6,1,1,1,0,1,0,1/ DATA (IRRED(18,I),I=-1,6)/6,1,1,0,1,1,0,1/ DATA (IRRED(19,I),I=-1,6)/6,1,0,0,0,0,1,1/ DATA (IRRED(20,I),I=-1,6)/6,1,1,1,0,0,1,1/ DATA (IRRED(21,I),I=-1,6)/6,1,0,1,1,0,1,1/ DATA (IRRED(22,I),I=-1,6)/6,1,1,0,0,1,1,1/ DATA (IRRED(23,I),I=-1,6)/6,1,0,1,0,1,1,1/ DATA (IRRED(24,I),I=-1,7)/7,1,1,0,0,0,0,0,1/ DATA (IRRED(25,I),I=-1,7)/7,1,0,0,1,0,0,0,1/ DATA (IRRED(26,I),I=-1,7)/7,1,1,1,1,0,0,0,1/ DATA (IRRED(27,I),I=-1,7)/7,1,0,0,0,1,0,0,1/ DATA (IRRED(28,I),I=-1,7)/7,1,0,1,1,1,0,0,1/ DATA (IRRED(29,I),I=-1,7)/7,1,1,1,0,0,1,0,1/ DATA (IRRED(30,I),I=-1,7)/7,1,1,0,1,0,1,0,1/ DATA (IRRED(31,I),I=-1,7)/7,1,0,0,1,1,1,0,1/ DATA (IRRED(32,I),I=-1,7)/7,1,1,1,1,1,1,0,1/ DATA (IRRED(33,I),I=-1,7)/7,1,0,0,0,0,0,1,1/ DATA (IRRED(34,I),I=-1,7)/7,1,1,0,1,0,0,1,1/ DATA (IRRED(35,I),I=-1,7)/7,1,1,0,0,1,0,1,1/ DATA (IRRED(36,I),I=-1,7)/7,1,0,1,0,1,0,1,1/ DATA (IRRED(37,I),I=-1,7)/7,1,0,1,0,0,1,1,1/ DATA (IRRED(38,I),I=-1,7)/7,1,1,1,1,0,1,1,1/ DATA (IRRED(39,I),I=-1,7)/7,1,0,0,0,1,1,1,1/ DATA (IRRED(40,I),I=-1,7)/7,1,1,1,0,1,1,1,1/ DATA (IRRED(41,I),I=-1,7)/7,1,0,1,1,1,1,1,1/ DATA (IRRED(42,I),I=-1,8)/8,1,1,0,1,1,0,0,0,1/ DATA (IRRED(43,I),I=-1,8)/8,1,0,1,1,1,0,0,0,1/ DATA (IRRED(44,I),I=-1,8)/8,1,1,0,1,0,1,0,0,1/ DATA (IRRED(45,I),I=-1,8)/8,1,0,1,1,0,1,0,0,1/ DATA (IRRED(46,I),I=-1,8)/8,1,0,0,1,1,1,0,0,1/ DATA (IRRED(47,I),I=-1,8)/8,1,1,1,1,1,1,0,0,1/ DATA (IRRED(48,I),I=-1,8)/8,1,0,1,1,0,0,1,0,1/ DATA (IRRED(49,I),I=-1,8)/8,1,1,1,1,1,0,1,0,1/ DATA (IRRED(50,I),I=-1,8)/8,1,1,0,0,0,1,1,0,1/ DATA (IRRED(51,I),I=-1,8)/8,1,0,1,0,0,1,1,0,1/ DATA (IRRED(52,I),I=-1,8)/8,1,0,0,1,0,1,1,0,1/ DATA (IRRED(53,I),I=-1,8)/8,1,0,0,0,1,1,1,0,1/ DATA (IRRED(54,I),I=-1,8)/8,1,1,1,0,1,1,1,0,1/ DATA (IRRED(55,I),I=-1,8)/8,1,1,0,1,1,1,1,0,1/ DATA (IRRED(56,I),I=-1,8)/8,1,1,1,0,0,0,0,1,1/ DATA (IRRED(57,I),I=-1,8)/8,1,1,0,1,0,0,0,1,1/ DATA (IRRED(58,I),I=-1,8)/8,1,0,1,1,0,0,0,1,1/ DATA (IRRED(59,I),I=-1,8)/8,1,1,1,1,1,0,0,1,1/ DATA (IRRED(60,I),I=-1,8)/8,1,1,0,0,0,1,0,1,1/ DATA (IRRED(61,I),I=-1,8)/8,1,0,0,1,0,1,0,1,1/ DATA (IRRED(62,I),I=-1,8)/8,1,0,0,0,1,1,0,1,1/ DATA (IRRED(63,I),I=-1,8)/8,1,0,1,1,1,1,0,1,1/ DATA (IRRED(64,I),I=-1,8)/8,1,1,0,0,0,0,1,1,1/ DATA (IRRED(65,I),I=-1,8)/8,1,1,1,1,0,0,1,1,1/ DATA (IRRED(66,I),I=-1,8)/8,1,1,1,0,1,0,1,1,1/ DATA (IRRED(67,I),I=-1,8)/8,1,0,1,1,1,0,1,1,1/ DATA (IRRED(68,I),I=-1,8)/8,1,1,1,0,0,1,1,1,1/ DATA (IRRED(69,I),I=-1,8)/8,1,1,0,0,1,1,1,1,1/ DATA (IRRED(70,I),I=-1,8)/8,1,0,1,0,1,1,1,1,1/ DATA (IRRED(71,I),I=-1,8)/8,1,0,0,1,1,1,1,1,1/ DATA (IRRED(72,I),I=-1,9)/9,1,1,0,0,0,0,0,0,0,1/ DATA (IRRED(73,I),I=-1,9)/9,1,0,0,0,1,0,0,0,0,1/ DATA (IRRED(74,I),I=-1,9)/9,1,1,1,0,1,0,0,0,0,1/ DATA (IRRED(75,I),I=-1,9)/9,1,1,0,1,1,0,0,0,0,1/ DATA (IRRED(76,I),I=-1,9)/9,1,0,0,0,0,1,0,0,0,1/ DATA (IRRED(77,I),I=-1,9)/9,1,0,1,1,0,1,0,0,0,1/ DATA (IRRED(78,I),I=-1,9)/9,1,1,0,0,1,1,0,0,0,1/ DATA (IRRED(79,I),I=-1,9)/9,1,1,0,1,0,0,1,0,0,1/ DATA (IRRED(80,I),I=-1,9)/9,1,0,0,1,1,0,1,0,0,1/ DATA (IRRED(81,I),I=-1,9)/9,1,1,1,1,1,0,1,0,0,1/ DATA (IRRED(82,I),I=-1,9)/9,1,0,1,0,0,1,1,0,0,1/ DATA (IRRED(83,I),I=-1,9)/9,1,0,0,1,0,1,1,0,0,1/ DATA (IRRED(84,I),I=-1,9)/9,1,1,1,1,0,1,1,0,0,1/ DATA (IRRED(85,I),I=-1,9)/9,1,1,1,0,1,1,1,0,0,1/ DATA (IRRED(86,I),I=-1,9)/9,1,0,1,1,1,1,1,0,0,1/ DATA (IRRED(87,I),I=-1,9)/9,1,1,1,0,0,0,0,1,0,1/ DATA (IRRED(88,I),I=-1,9)/9,1,0,1,0,1,0,0,1,0,1/ DATA (IRRED(89,I),I=-1,9)/9,1,0,0,1,1,0,0,1,0,1/ DATA (IRRED(90,I),I=-1,9)/9,1,1,0,0,0,1,0,1,0,1/ DATA (IRRED(91,I),I=-1,9)/9,1,0,1,0,0,1,0,1,0,1/ DATA (IRRED(92,I),I=-1,9)/9,1,1,1,1,0,1,0,1,0,1/ DATA (IRRED(93,I),I=-1,9)/9,1,1,1,0,1,1,0,1,0,1/ DATA (IRRED(94,I),I=-1,9)/9,1,0,1,1,1,1,0,1,0,1/ DATA (IRRED(95,I),I=-1,9)/9,1,1,1,1,0,0,1,1,0,1/ DATA (IRRED(96,I),I=-1,9)/9,1,0,0,0,1,0,1,1,0,1/ DATA (IRRED(97,I),I=-1,9)/9,1,1,0,1,1,0,1,1,0,1/ DATA (IRRED(98,I),I=-1,9)/9,1,0,1,0,1,1,1,1,0,1/ DATA (IRRED(99,I),I=-1,9)/9,1,0,0,1,1,1,1,1,0,1/ DATA (IRRED(100,I),I=-1,9)/9,1,0,0,0,0,0,0,0,1,1/ DATA (IRRED(101,I),I=-1,9)/9,1,1,0,0,1,0,0,0,1,1/ DATA (IRRED(102,I),I=-1,9)/9,1,0,1,0,1,0,0,0,1,1/ DATA (IRRED(103,I),I=-1,9)/9,1,1,1,1,1,0,0,0,1,1/ DATA (IRRED(104,I),I=-1,9)/9,1,1,0,0,0,1,0,0,1,1/ DATA (IRRED(105,I),I=-1,9)/9,1,0,0,0,1,1,0,0,1,1/ DATA (IRRED(106,I),I=-1,9)/9,1,1,0,1,1,1,0,0,1,1/ DATA (IRRED(107,I),I=-1,9)/9,1,0,0,1,0,0,1,0,1,1/ DATA (IRRED(108,I),I=-1,9)/9,1,1,1,1,0,0,1,0,1,1/ DATA (IRRED(109,I),I=-1,9)/9,1,1,0,1,1,0,1,0,1,1/ DATA (IRRED(110,I),I=-1,9)/9,1,0,0,0,0,1,1,0,1,1/ DATA (IRRED(111,I),I=-1,9)/9,1,1,0,1,0,1,1,0,1,1/ DATA (IRRED(112,I),I=-1,9)/9,1,0,1,1,0,1,1,0,1,1/ DATA (IRRED(113,I),I=-1,9)/9,1,1,0,0,1,1,1,0,1,1/ DATA (IRRED(114,I),I=-1,9)/9,1,1,1,1,1,1,1,0,1,1/ DATA (IRRED(115,I),I=-1,9)/9,1,0,1,0,0,0,0,1,1,1/ DATA (IRRED(116,I),I=-1,9)/9,1,1,1,1,0,0,0,1,1,1/ DATA (IRRED(117,I),I=-1,9)/9,1,0,0,0,0,1,0,1,1,1/ DATA (IRRED(118,I),I=-1,9)/9,1,0,1,0,1,1,0,1,1,1/ DATA (IRRED(119,I),I=-1,9)/9,1,0,0,1,1,1,0,1,1,1/ DATA (IRRED(120,I),I=-1,9)/9,1,1,1,0,0,0,1,1,1,1/ DATA (IRRED(121,I),I=-1,9)/9,1,1,0,1,0,0,1,1,1,1/ DATA (IRRED(122,I),I=-1,9)/9,1,0,1,1,0,0,1,1,1,1/ DATA (IRRED(123,I),I=-1,9)/9,1,0,1,0,1,0,1,1,1,1/ DATA (IRRED(124,I),I=-1,9)/9,1,0,0,1,1,0,1,1,1,1/ DATA (IRRED(125,I),I=-1,9)/9,1,1,0,0,0,1,1,1,1,1/ DATA (IRRED(126,I),I=-1,9)/9,1,0,0,1,0,1,1,1,1,1/ DATA (IRRED(127,I),I=-1,9)/9,1,1,0,1,1,1,1,1,1,1/ DATA (IRRED(128,I),I=-1,10)/10,1,0,0,1,0,0,0,0,0,0,1/ DATA (IRRED(129,I),I=-1,10)/10,1,1,1,1,0,0,0,0,0,0,1/ DATA (IRRED(130,I),I=-1,10)/10,1,1,0,1,1,0,0,0,0,0,1/ DATA (IRRED(131,I),I=-1,10)/10,1,0,1,1,1,0,0,0,0,0,1/ DATA (IRRED(132,I),I=-1,10)/10,1,1,1,0,0,1,0,0,0,0,1/ DATA (IRRED(133,I),I=-1,10)/10,1,0,1,1,0,1,0,0,0,0,1/ DATA (IRRED(134,I),I=-1,10)/10,1,0,1,0,1,1,0,0,0,0,1/ DATA (IRRED(135,I),I=-1,10)/10,1,1,1,0,0,0,1,0,0,0,1/ DATA (IRRED(136,I),I=-1,10)/10,1,1,0,0,1,0,1,0,0,0,1/ DATA (IRRED(137,I),I=-1,10)/10,1,1,0,0,0,1,1,0,0,0,1/ DATA (IRRED(138,I),I=-1,10)/10,1,0,1,0,0,1,1,0,0,0,1/ DATA (IRRED(139,I),I=-1,10)/10,1,1,1,1,0,1,1,0,0,0,1/ DATA (IRRED(140,I),I=-1,10)/10,1,0,0,0,0,0,0,1,0,0,1/ DATA (IRRED(141,I),I=-1,10)/10,1,1,0,1,0,0,0,1,0,0,1/ DATA (IRRED(142,I),I=-1,10)/10,1,0,0,1,1,0,0,1,0,0,1/ DATA (IRRED(143,I),I=-1,10)/10,1,0,0,1,0,1,0,1,0,0,1/ DATA (IRRED(144,I),I=-1,10)/10,1,1,1,1,0,1,0,1,0,0,1/ DATA (IRRED(145,I),I=-1,10)/10,1,0,1,0,0,0,1,1,0,0,1/ DATA (IRRED(146,I),I=-1,10)/10,1,0,0,1,0,0,1,1,0,0,1/ DATA (IRRED(147,I),I=-1,10)/10,1,1,1,0,1,0,1,1,0,0,1/ DATA (IRRED(148,I),I=-1,10)/10,1,1,1,0,0,1,1,1,0,0,1/ DATA (IRRED(149,I),I=-1,10)/10,1,0,1,1,0,1,1,1,0,0,1/ DATA (IRRED(150,I),I=-1,10)/10,1,1,0,0,1,1,1,1,0,0,1/ DATA (IRRED(151,I),I=-1,10)/10,1,1,1,1,1,1,1,1,0,0,1/ DATA (IRRED(152,I),I=-1,10)/10,1,1,0,1,0,0,0,0,1,0,1/ DATA (IRRED(153,I),I=-1,10)/10,1,0,1,1,0,0,0,0,1,0,1/ DATA (IRRED(154,I),I=-1,10)/10,1,0,0,1,1,0,0,0,1,0,1/ DATA (IRRED(155,I),I=-1,10)/10,1,1,1,1,1,0,0,0,1,0,1/ DATA (IRRED(156,I),I=-1,10)/10,1,1,0,0,0,1,0,0,1,0,1/ DATA (IRRED(157,I),I=-1,10)/10,1,0,0,0,1,1,0,0,1,0,1/ DATA (IRRED(158,I),I=-1,10)/10,1,0,1,1,1,1,0,0,1,0,1/ DATA (IRRED(159,I),I=-1,10)/10,1,1,0,0,0,0,1,0,1,0,1/ DATA (IRRED(160,I),I=-1,10)/10,1,1,1,0,1,0,1,0,1,0,1/ DATA (IRRED(161,I),I=-1,10)/10,1,0,0,0,0,1,1,0,1,0,1/ DATA (IRRED(162,I),I=-1,10)/10,1,1,1,0,0,1,1,0,1,0,1/ DATA (IRRED(163,I),I=-1,10)/10,1,1,0,1,0,1,1,0,1,0,1/ DATA (IRRED(164,I),I=-1,10)/10,1,0,1,0,0,0,0,1,1,0,1/ DATA (IRRED(165,I),I=-1,10)/10,1,1,1,1,0,0,0,1,1,0,1/ DATA (IRRED(166,I),I=-1,10)/10,1,1,1,0,1,0,0,1,1,0,1/ DATA (IRRED(167,I),I=-1,10)/10,1,1,0,1,1,0,0,1,1,0,1/ DATA (IRRED(168,I),I=-1,10)/10,1,0,0,0,0,1,0,1,1,0,1/ DATA (IRRED(169,I),I=-1,10)/10,1,1,0,1,0,1,0,1,1,0,1/ DATA (IRRED(170,I),I=-1,10)/10,1,0,0,1,1,1,0,1,1,0,1/ DATA (IRRED(171,I),I=-1,10)/10,1,0,0,0,0,0,1,1,1,0,1/ DATA (IRRED(172,I),I=-1,10)/10,1,1,1,0,0,0,1,1,1,0,1/ DATA (IRRED(173,I),I=-1,10)/10,1,0,1,0,0,1,1,1,1,0,1/ DATA (IRRED(174,I),I=-1,10)/10,1,1,1,0,1,1,1,1,1,0,1/ DATA (IRRED(175,I),I=-1,10)/10,1,1,0,1,1,1,1,1,1,0,1/ DATA (IRRED(176,I),I=-1,10)/10,1,1,0,0,1,0,0,0,0,1,1/ DATA (IRRED(177,I),I=-1,10)/10,1,0,1,0,1,0,0,0,0,1,1/ DATA (IRRED(178,I),I=-1,10)/10,1,1,0,0,0,1,0,0,0,1,1/ DATA (IRRED(179,I),I=-1,10)/10,1,0,1,0,0,1,0,0,0,1,1/ DATA (IRRED(180,I),I=-1,10)/10,1,0,0,0,1,1,0,0,0,1,1/ DATA (IRRED(181,I),I=-1,10)/10,1,1,1,0,1,1,0,0,0,1,1/ DATA (IRRED(182,I),I=-1,10)/10,1,1,0,0,0,0,1,0,0,1,1/ DATA (IRRED(183,I),I=-1,10)/10,1,1,1,1,0,0,1,0,0,1,1/ DATA (IRRED(184,I),I=-1,10)/10,1,0,0,0,1,0,1,0,0,1,1/ DATA (IRRED(185,I),I=-1,10)/10,1,1,0,1,1,0,1,0,0,1,1/ DATA (IRRED(186,I),I=-1,10)/10,1,0,0,1,1,1,1,0,0,1,1/ DATA (IRRED(187,I),I=-1,10)/10,1,1,1,1,1,1,1,0,0,1,1/ DATA (IRRED(188,I),I=-1,10)/10,1,0,1,0,0,0,0,1,0,1,1/ DATA (IRRED(189,I),I=-1,10)/10,1,0,0,1,0,0,0,1,0,1,1/ DATA (IRRED(190,I),I=-1,10)/10,1,1,1,0,0,1,0,1,0,1,1/ DATA (IRRED(191,I),I=-1,10)/10,1,0,1,1,0,1,0,1,0,1,1/ DATA (IRRED(192,I),I=-1,10)/10,1,0,1,0,1,1,0,1,0,1,1/ DATA (IRRED(193,I),I=-1,10)/10,1,1,1,1,1,1,0,1,0,1,1/ DATA (IRRED(194,I),I=-1,10)/10,1,0,0,0,0,0,1,1,0,1,1/ DATA (IRRED(195,I),I=-1,10)/10,1,0,1,1,0,0,1,1,0,1,1/ DATA (IRRED(196,I),I=-1,10)/10,1,1,0,0,1,0,1,1,0,1,1/ DATA (IRRED(197,I),I=-1,10)/10,1,1,1,1,1,0,1,1,0,1,1/ DATA (IRRED(198,I),I=-1,10)/10,1,1,1,0,1,1,1,1,0,1,1/ DATA (IRRED(199,I),I=-1,10)/10,1,0,1,1,1,1,1,1,0,1,1/ DATA (IRRED(200,I),I=-1,10)/10,1,1,1,1,0,0,0,0,1,1,1/ DATA (IRRED(201,I),I=-1,10)/10,1,0,0,0,1,0,0,0,1,1,1/ DATA (IRRED(202,I),I=-1,10)/10,1,1,1,0,1,0,0,0,1,1,1/ DATA (IRRED(203,I),I=-1,10)/10,1,0,1,1,1,0,0,0,1,1,1/ DATA (IRRED(204,I),I=-1,10)/10,1,0,0,0,0,1,0,0,1,1,1/ DATA (IRRED(205,I),I=-1,10)/10,1,1,0,1,0,1,0,0,1,1,1/ DATA (IRRED(206,I),I=-1,10)/10,1,0,1,0,1,1,0,0,1,1,1/ DATA (IRRED(207,I),I=-1,10)/10,1,0,0,1,1,1,0,0,1,1,1/ DATA (IRRED(208,I),I=-1,10)/10,1,1,1,0,0,0,1,0,1,1,1/ DATA (IRRED(209,I),I=-1,10)/10,1,0,1,1,0,0,1,0,1,1,1/ DATA (IRRED(210,I),I=-1,10)/10,1,0,1,0,1,0,1,0,1,1,1/ DATA (IRRED(211,I),I=-1,10)/10,1,0,0,1,1,0,1,0,1,1,1/ DATA (IRRED(212,I),I=-1,10)/10,1,1,0,0,0,1,1,0,1,1,1/ DATA (IRRED(213,I),I=-1,10)/10,1,1,0,1,1,1,1,0,1,1,1/ DATA (IRRED(214,I),I=-1,10)/10,1,0,1,1,1,1,1,0,1,1,1/ DATA (IRRED(215,I),I=-1,10)/10,1,0,0,0,0,0,0,1,1,1,1/ DATA (IRRED(216,I),I=-1,10)/10,1,1,1,0,0,0,0,1,1,1,1/ DATA (IRRED(217,I),I=-1,10)/10,1,0,1,1,0,0,0,1,1,1,1/ DATA (IRRED(218,I),I=-1,10)/10,1,1,0,0,1,0,0,1,1,1,1/ DATA (IRRED(219,I),I=-1,10)/10,1,0,0,1,0,1,0,1,1,1,1/ DATA (IRRED(220,I),I=-1,10)/10,1,0,0,0,1,1,0,1,1,1,1/ DATA (IRRED(221,I),I=-1,10)/10,1,0,1,0,0,0,1,1,1,1,1/ DATA (IRRED(222,I),I=-1,10)/10,1,1,0,1,1,0,1,1,1,1,1/ DATA (IRRED(223,I),I=-1,10)/10,1,1,0,1,0,1,1,1,1,1,1/ DATA (IRRED(224,I),I=-1,10)/10,1,1,0,0,1,1,1,1,1,1,1/ DATA (IRRED(225,I),I=-1,10)/10,1,0,0,1,1,1,1,1,1,1,1/ DATA (IRRED(226,I),I=-1,10)/10,1,1,1,1,1,1,1,1,1,1,1/ DATA (IRRED(227,I),I=-1,11)/11,1,0,1,0,0,0,0,0,0,0,0,1/ DATA (IRRED(228,I),I=-1,11)/11,1,1,1,0,1,0,0,0,0,0,0,1/ DATA (IRRED(229,I),I=-1,11)/11,1,1,0,1,0,1,0,0,0,0,0,1/ DATA (IRRED(230,I),I=-1,11)/11,1,0,1,1,0,1,0,0,0,0,0,1/ DATA (IRRED(231,I),I=-1,11)/11,1,1,1,0,0,0,1,0,0,0,0,1/ DATA (IRRED(232,I),I=-1,11)/11,1,1,0,0,0,1,1,0,0,0,0,1/ DATA (IRRED(233,I),I=-1,11)/11,1,0,1,0,0,1,1,0,0,0,0,1/ DATA (IRRED(234,I),I=-1,11)/11,1,0,0,0,1,1,1,0,0,0,0,1/ DATA (IRRED(235,I),I=-1,11)/11,1,1,0,1,1,1,1,0,0,0,0,1/ DATA (IRRED(236,I),I=-1,11)/11,1,0,1,1,0,0,0,1,0,0,0,1/ DATA (IRRED(237,I),I=-1,11)/11,1,0,1,0,1,0,0,1,0,0,0,1/ DATA (IRRED(238,I),I=-1,11)/11,1,1,1,1,1,0,0,1,0,0,0,1/ DATA (IRRED(239,I),I=-1,11)/11,1,0,0,1,0,1,0,1,0,0,0,1/ DATA (IRRED(240,I),I=-1,11)/11,1,0,0,0,1,1,0,1,0,0,0,1/ DATA (IRRED(241,I),I=-1,11)/11,1,1,0,0,0,0,1,1,0,0,0,1/ DATA (IRRED(242,I),I=-1,11)/11,1,1,1,1,0,0,1,1,0,0,0,1/ DATA (IRRED(243,I),I=-1,11)/11,1,0,0,0,1,0,1,1,0,0,0,1/ DATA (IRRED(244,I),I=-1,11)/11,1,0,0,0,0,1,1,1,0,0,0,1/ DATA (IRRED(245,I),I=-1,11)/11,1,1,1,0,0,1,1,1,0,0,0,1/ DATA (IRRED(246,I),I=-1,11)/11,1,1,0,1,0,1,1,1,0,0,0,1/ DATA (IRRED(247,I),I=-1,11)/11,1,0,1,0,1,1,1,1,0,0,0,1/ DATA (IRRED(248,I),I=-1,11)/11,1,0,1,1,0,0,0,0,1,0,0,1/ DATA (IRRED(249,I),I=-1,11)/11,1,1,0,0,1,0,0,0,1,0,0,1/ DATA (IRRED(250,I),I=-1,11)/11,1,0,1,0,0,1,0,0,1,0,0,1/ DATA (IRRED(251,I),I=-1,11)/11,1,0,0,1,0,1,0,0,1,0,0,1/ DATA (IRRED(252,I),I=-1,11)/11,1,1,1,0,1,1,0,0,1,0,0,1/ DATA (IRRED(253,I),I=-1,11)/11,1,1,0,1,1,1,0,0,1,0,0,1/ DATA (IRRED(254,I),I=-1,11)/11,1,0,1,1,1,1,0,0,1,0,0,1/ DATA (IRRED(255,I),I=-1,11)/11,1,0,1,0,0,0,1,0,1,0,0,1/ DATA (IRRED(256,I),I=-1,11)/11,1,0,0,1,0,0,1,0,1,0,0,1/ DATA (IRRED(257,I),I=-1,11)/11,1,0,0,0,1,0,1,0,1,0,0,1/ DATA (IRRED(258,I),I=-1,11)/11,1,1,0,1,1,0,1,0,1,0,0,1/ DATA (IRRED(259,I),I=-1,11)/11,1,1,0,0,1,1,1,0,1,0,0,1/ DATA (IRRED(260,I),I=-1,11)/11,1,0,1,0,1,1,1,0,1,0,0,1/ DATA (IRRED(261,I),I=-1,11)/11,1,1,1,1,1,1,1,0,1,0,0,1/ DATA (IRRED(262,I),I=-1,11)/11,1,1,0,0,0,0,0,1,1,0,0,1/ DATA (IRRED(263,I),I=-1,11)/11,1,1,1,1,0,0,0,1,1,0,0,1/ DATA (IRRED(264,I),I=-1,11)/11,1,1,0,1,0,1,0,1,1,0,0,1/ DATA (IRRED(265,I),I=-1,11)/11,1,0,1,1,0,1,0,1,1,0,0,1/ DATA (IRRED(266,I),I=-1,11)/11,1,0,0,1,1,1,0,1,1,0,0,1/ DATA (IRRED(267,I),I=-1,11)/11,1,1,1,0,0,0,1,1,1,0,0,1/ DATA (IRRED(268,I),I=-1,11)/11,1,0,0,1,1,0,1,1,1,0,0,1/ DATA (IRRED(269,I),I=-1,11)/11,1,0,1,0,0,1,1,1,1,0,0,1/ DATA (IRRED(270,I),I=-1,11)/11,1,1,1,1,0,1,1,1,1,0,0,1/ DATA (IRRED(271,I),I=-1,11)/11,1,1,1,0,1,1,1,1,1,0,0,1/ DATA (IRRED(272,I),I=-1,11)/11,1,0,0,0,0,0,0,0,0,1,0,1/ DATA (IRRED(273,I),I=-1,11)/11,1,1,1,0,0,0,0,0,0,1,0,1/ DATA (IRRED(274,I),I=-1,11)/11,1,1,0,0,1,0,0,0,0,1,0,1/ DATA (IRRED(275,I),I=-1,11)/11,1,0,1,0,1,0,0,0,0,1,0,1/ DATA (IRRED(276,I),I=-1,11)/11,1,0,0,1,0,1,0,0,0,1,0,1/ DATA (IRRED(277,I),I=-1,11)/11,1,0,0,1,0,0,1,0,0,1,0,1/ DATA (IRRED(278,I),I=-1,11)/11,1,0,0,0,0,1,1,0,0,1,0,1/ DATA (IRRED(279,I),I=-1,11)/11,1,0,1,1,0,1,1,0,0,1,0,1/ DATA (IRRED(280,I),I=-1,11)/11,1,0,0,1,1,1,1,0,0,1,0,1/ DATA (IRRED(281,I),I=-1,11)/11,1,1,1,1,1,1,1,0,0,1,0,1/ DATA (IRRED(282,I),I=-1,11)/11,1,0,1,0,0,0,0,1,0,1,0,1/ DATA (IRRED(283,I),I=-1,11)/11,1,0,0,0,1,0,0,1,0,1,0,1/ DATA (IRRED(284,I),I=-1,11)/11,1,0,1,1,1,0,0,1,0,1,0,1/ DATA (IRRED(285,I),I=-1,11)/11,1,1,1,0,0,1,0,1,0,1,0,1/ DATA (IRRED(286,I),I=-1,11)/11,1,1,0,1,0,1,0,1,0,1,0,1/ DATA (IRRED(287,I),I=-1,11)/11,1,1,0,0,1,1,0,1,0,1,0,1/ DATA (IRRED(288,I),I=-1,11)/11,1,0,1,0,1,1,0,1,0,1,0,1/ DATA (IRRED(289,I),I=-1,11)/11,1,0,1,0,1,0,1,1,0,1,0,1/ DATA (IRRED(290,I),I=-1,11)/11,1,1,1,1,1,0,1,1,0,1,0,1/ DATA (IRRED(291,I),I=-1,11)/11,1,1,0,0,0,1,1,1,0,1,0,1/ DATA (IRRED(292,I),I=-1,11)/11,1,0,0,1,0,1,1,1,0,1,0,1/ DATA (IRRED(293,I),I=-1,11)/11,1,1,1,1,0,1,1,1,0,1,0,1/ DATA (IRRED(294,I),I=-1,11)/11,1,0,0,0,1,1,1,1,0,1,0,1/ DATA (IRRED(295,I),I=-1,11)/11,1,1,0,1,1,1,1,1,0,1,0,1/ DATA (IRRED(296,I),I=-1,11)/11,1,1,0,0,0,0,0,0,1,1,0,1/ DATA (IRRED(297,I),I=-1,11)/11,1,0,0,1,0,0,0,0,1,1,0,1/ DATA (IRRED(298,I),I=-1,11)/11,1,0,0,0,1,0,0,0,1,1,0,1/ DATA (IRRED(299,I),I=-1,11)/11,1,1,0,0,1,1,0,0,1,1,0,1/ DATA (IRRED(300,I),I=-1,11)/11,1,1,1,1,1,1,0,0,1,1,0,1/ DATA (IRRED(301,I),I=-1,11)/11,1,0,0,0,0,0,1,0,1,1,0,1/ DATA (IRRED(302,I),I=-1,11)/11,1,1,0,1,0,0,1,0,1,1,0,1/ DATA (IRRED(303,I),I=-1,11)/11,1,0,0,1,1,0,1,0,1,1,0,1/ DATA (IRRED(304,I),I=-1,11)/11,1,1,1,1,1,0,1,0,1,1,0,1/ DATA (IRRED(305,I),I=-1,11)/11,1,0,1,0,0,1,1,0,1,1,0,1/ DATA (IRRED(306,I),I=-1,11)/11,1,1,1,1,0,1,1,0,1,1,0,1/ DATA (IRRED(307,I),I=-1,11)/11,1,0,1,1,1,1,1,0,1,1,0,1/ DATA (IRRED(308,I),I=-1,11)/11,1,1,1,0,0,0,0,1,1,1,0,1/ DATA (IRRED(309,I),I=-1,11)/11,1,1,0,1,0,0,0,1,1,1,0,1/ DATA (IRRED(310,I),I=-1,11)/11,1,1,0,0,1,0,0,1,1,1,0,1/ DATA (IRRED(311,I),I=-1,11)/11,1,0,1,0,1,0,0,1,1,1,0,1/ DATA (IRRED(312,I),I=-1,11)/11,1,1,1,1,0,1,0,1,1,1,0,1/ DATA (IRRED(313,I),I=-1,11)/11,1,1,1,0,1,1,0,1,1,1,0,1/ DATA (IRRED(314,I),I=-1,11)/11,1,0,1,1,1,1,0,1,1,1,0,1/ DATA (IRRED(315,I),I=-1,11)/11,1,0,0,1,0,0,1,1,1,1,0,1/ DATA (IRRED(316,I),I=-1,11)/11,1,1,0,1,1,0,1,1,1,1,0,1/ DATA (IRRED(317,I),I=-1,11)/11,1,0,1,1,1,0,1,1,1,1,0,1/ DATA (IRRED(318,I),I=-1,11)/11,1,1,1,0,0,1,1,1,1,1,0,1/ C .. C C all the 10 th degree polynomials for mod 2 are used; n=10 C the first 92, 11 th degree polyn. are also included; n=11. C niederreiter-lidl book, pg 385 - the first 3 columns are done C C Prepare to work in Z2 C CALL SETFLD(2) C DO 60 I = 1,DIMEN C C For each dimension, we need to calculate powers of an C appropriate irreducible polynomial : see Niederreiter C page 65, just below equation (19). C Copy the appropriate irreducible polynomial into PX, C and its degree into E. Set polynomial B = PX ** 0 = 1. C M is the degree of B. Subsequently B will hold higher C powers of PX. C E = IRRED(I,DEG) DO 10 J = -1,E PX(J) = IRRED(I,J) 10 CONTINUE B(DEG) = 0 B(0) = 1 C C Niederreiter (page 56, after equation (7), defines two C variables Q and U. We do not need Q explicitly, but we C do need U. C U = 0 C DO 30 J = 1,NBITS C C If U = 0, we need to set B to the next power of PX C and recalculate V. This is done by subroutine CALCV. C IF (U.EQ.0) CALL CALCV(PX,B,V,MAXV) C C Now C is obtained from V. Niederreiter C obtains A from V (page 65, near the bottom), and then gets C C from A (page 56, equation (7)). However this can be done C in one step. Here CI(J,R) corresponds to C Niederreiter's C(I,J,R). C DO 20 R = 0,NBITS - 1 CI(J,R) = V(R+U) 20 CONTINUE C C Increment U. If U = E, then U = 0 and in Niederreiter's C paper Q = Q + 1. Here, however, Q is not used explicitly. C U = U + 1 IF (U.EQ.E) U = 0 30 CONTINUE C C The array CI now holds the values of C(I,J,R) for this value C of I. We pack them into array CJ so that CJ(I,R) holds all C the values of C(I,J,R) for J from 1 to NBITS. C DO 50 R = 0,NBITS - 1 TERM = 0 DO 40 J = 1,NBITS TERM = 2*TERM + CI(J,R) 40 CONTINUE CJ(I,R) = TERM 50 CONTINUE C 60 CONTINUE C C COMPUTING GENERATOR MATRICES OF USER CHOICE C IF (IFLAG.EQ.0) THEN DO 80 I = 1,DIMEN DO 70 J = 0,NBITS - 1 SCJ(I,J) = CJ(I,J) 70 CONTINUE SHIFT(I) = 0 80 CONTINUE ELSE IF ((IFLAG.EQ.1) .OR. (IFLAG.EQ.3)) THEN CALL GNCRML(MAXDIM,NBITS,DIMEN,LSM,SHIFT) DO 120 I = 1,DIMEN DO 110 J = 0,NBITS - 1 L = 1 TEMP2 = 0 DO 100 P = NBITS - 1,0,-1 TEMP1 = 0 DO 90 K = 0,NBITS - 1 TEMP1 = TEMP1 + (IBITS(LSM(I,P),K,1)* + IBITS(CJ(I,J),K,1)) 90 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L L = 2*L 100 CONTINUE SCJ(I,J) = TEMP2 110 CONTINUE 120 CONTINUE END IF IF ((IFLAG.EQ.2) .OR. (IFLAG.EQ.3)) THEN CALL GNCRMU(NBITS,USM,USHIFT) DO 180 I = 1,DIMEN DO 140 J = 0,NBITS - 1 P = NBITS - 1 DO 130 K = 0,NBITS - 1 IF (IFLAG.EQ.2) THEN TV(I,P,J) = IBITS(CJ(I,J),K,1) ELSE TV(I,P,J) = IBITS(SCJ(I,J),K,1) END IF P = P - 1 130 CONTINUE 140 CONTINUE DO 170 PP = 0,NBITS - 1 TEMP2 = 0 TEMP4 = 0 L = 1 DO 160 J = NBITS - 1,0,-1 TEMP1 = 0 TEMP3 = 0 DO 150 P = 0,NBITS - 1 TEMP1 = TEMP1 + TV(I,J,P)*USM(P,PP) IF (PP.EQ.0) THEN TEMP3 = TEMP3 + TV(I,J,P)*USHIFT(P) END IF 150 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L IF (PP.EQ.1) THEN TEMP3 = MOD(TEMP1,2) TEMP4 = TEMP4 + TEMP3*L END IF L = 2*L 160 CONTINUE SCJ(I,PP) = TEMP2 IF (PP.EQ.0) THEN IF (IFLAG.EQ.3) THEN SHIFT(I) = EXOR(TEMP4,SHIFT(I)) ELSE SHIFT(I) = TEMP4 END IF END IF 170 CONTINUE 180 CONTINUE END IF END IF DO 190 I = 1,DIMEN NEXTQ(I) = SHIFT(I) 190 CONTINUE RETURN END SUBROUTINE SETFLD(QIN) C C This version : 12 December 1991 C C This subroutine sets up addition, multiplication, and C subtraction tables for the finite field of order QIN. C If necessary, it reads precalculated tables from the file C 'gftabs.dat' using unit 1. These precalculated tables C are supposed to have been put there by GFARIT. C C ***** For the base-2 programs, these precalculated C ***** tables are not needed and, therefore, neither C ***** is GFARIT. C C C Unit 1 is closed both before and after the call of SETFLD. C C USES C Integer function CHARAC gets the characteristic of a field. C C C ------------------------------------------------------------ C C The following COMMON block, used by many subroutines, C gives the order Q of a field, its characteristic P, and its C addition, multiplication and subtraction tables. C The parameter MAXQ gives the order of the largest field to C be handled. C C C The following definitions concern the representation of C polynomials. C C C The parameter MAXDEG gives the highest degree of polynomial C to be handled. Polynomials stored as arrays have the C coefficient of degree n in POLY(N), and the degree of the C polynomial in POLY(-1). The parameter DEG is just to remind C us of this last fact. A polynomial which is identically 0 C is given degree -1. C C A polynomial can also be stored in an integer I, with C I = AN*Q**N + ... + A0. C Routines ITOP and PTOI convert between these two formats. C C ----------------------------------------------------------- C C C C C .. Parameters .. INTEGER MAXQ PARAMETER (MAXQ=50) C .. C .. Scalar Arguments .. INTEGER QIN C .. C .. Scalars in Common .. INTEGER P,Q C .. C .. Arrays in Common .. INTEGER ADD(0:MAXQ-1,0:MAXQ-1),MUL(0:MAXQ-1,0:MAXQ-1), + SUB(0:MAXQ-1,0:MAXQ-1) C .. C .. Local Scalars .. INTEGER I,J,N C .. C .. External Functions .. INTEGER CHARAC EXTERNAL CHARAC C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /FIELD/P,Q,ADD,MUL,SUB C .. C .. Save statement .. SAVE /FIELD/ C .. IF (QIN.LE.1 .OR. QIN.GT.MAXQ) THEN WRITE (*,FMT=*) ' SETFLD : Bad value of Q' STOP END IF C Q = QIN P = CHARAC(Q) C IF (P.EQ.0) THEN WRITE (*,FMT=*) ' SETFLD : There is no field of order',Q STOP END IF C C Set up to handle a field of prime order : calculate ADD and MUL. C IF (P.EQ.Q) THEN DO 20 I = 0,Q - 1 DO 10 J = 0,Q - 1 ADD(I,J) = MOD(I+J,P) MUL(I,J) = MOD(I*J,P) 10 CONTINUE 20 CONTINUE C C Set up to handle a field of prime-power order : tables for C ADD and MUL are in the file 'gftabs.dat'. C ELSE OPEN (UNIT=1,FILE='gftabs.dat',STATUS='old') C C ***** OPEN statements are system dependent. C 30 READ (1,FMT=9000,END=80) N DO 40 I = 0,N - 1 READ (1,FMT=9000) (ADD(I,J),J=0,N-1) 40 CONTINUE DO 50 I = 0,N - 1 READ (1,FMT=9000) (MUL(I,J),J=0,N-1) 50 CONTINUE IF (N.NE.Q) GO TO 30 CLOSE (1) END IF C C Now use the addition table to set the subtraction table. C DO 70 I = 0,Q - 1 DO 60 J = 0,Q - 1 SUB(ADD(I,J),I) = J 60 CONTINUE 70 CONTINUE RETURN C 80 WRITE (*,FMT=*) ' SETFLD : Tables for q =',Q,' not found' STOP C 9000 FORMAT (20I3) END C C ***** end of SUBROUTINE PLYMUL SUBROUTINE SGOLO2(QUASI) C C This modified routine of GOLO2 C C C This subroutine generates a new quasi-random vector C on each call. C C INPUT C From SINLO2, labelled common /COMM2/, properly initialized. C C OUTPUT C To the caller, the next vector in the sequence in the C array QUASI. C C ------------------------------------------------------------ C C C This file defines the common block /COMM2/ and some C associated parameters. These are for use in the base 2 C version of the generator. C C C The parameter MAXDIM is the maximum dimension that will be used. C NBITS is the number of bits (not counting the sign) in a C fixed-point integer. C C C The common block /COMM2/ : C SCJ - Contains the packed values of Niederreiter's C(I,J,R) C DIMEN - The dimension of the sequence to be generated C COUNT - The index of the current item in the sequence, C expressed as an array of bits. COUNT(R) is the same C as Niederreiter's AR(N) (page 54) except that N is C implicit. C NEXTQ - The numerators of the next item in the series. These C are like Niederreiter's XI(N) (page 54) except that C N is implicit, and the NEXTQ are integers. To obtain C the values of XI(N), multiply by RECIP (see GOLO2). C C Array SCJ of the common block is set up by subroutine SCLCC2. C The other items in the common block are set up by INLO2. C C ------------------------------------------------------------ C C C C C The parameter RECIP is the multiplier which changes the C integers in NEXTQ into the required real values in QUASI. C C C Multiply the numerators in NEXTQ by RECIP to get the next C quasi-random vector C C .. Parameters .. INTEGER MAXDIM,NBITS PARAMETER (MAXDIM=318,NBITS=31) DOUBLE PRECISION RECIP PARAMETER (RECIP=2.0** (-NBITS)) C .. C .. Array Arguments .. DOUBLE PRECISION QUASI(*) C .. C .. Scalars in Common .. INTEGER COUNT,DIMEN C .. C .. Arrays in Common .. INTEGER NEXTQ(MAXDIM),SCJ(MAXDIM,0:NBITS-1) C .. C .. Local Scalars .. INTEGER I,R C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /COMM2/SCJ,DIMEN,COUNT,NEXTQ C .. C .. Save statement .. SAVE /COMM2/ C .. DO 10 I = 1,DIMEN QUASI(I) = NEXTQ(I)*RECIP 10 CONTINUE C C Find the position of the right-hand zero in COUNT. This C is the bit that changes in the Gray-code representation as C we go from COUNT to COUNT+1. C R = 0 I = COUNT 20 IF (MOD(I,2).NE.0) THEN R = R + 1 I = I/2 GO TO 20 END IF C C Check that we have not passed 2**NBITS calls on GOLO2 C IF (R.GE.NBITS) THEN WRITE (*,FMT=*) ' SGOLO2 : Too many calls' STOP END IF C C Compute the new numerators in vector NEXTQ C DO 30 I = 1,DIMEN NEXTQ(I) = EXOR(NEXTQ(I),SCJ(I,R)) C QUASI(I) = NEXTQ(I)*RECIP 30 CONTINUE C COUNT = COUNT + 1 RETURN END C C ***** end of INTEGER FUNCTION CHARAC SUBROUTINE SINLO2(DIM,SKIP,IFLAG) C C This Modified Routine of INLO2 C C C This subroutine calculates the values of Niederreiter's C C(I,J,R) and Various Scrambled Niederreiter's C(I,J,R) C performs other initialisation necessary C before calling GOLO2. C C INPUT : C DIMEN - The dimension of the sequence to be generated. C {DIMEN is called DIM in the argument of INLO2, C because DIMEN is subsequently passed via COMMON C and is called DIMEN there.} C C SKIP - The number of values to throw away at the beginning C of the sequence. C C IFLAG -User choice of Genertor matrices. C C C OUTPUT : C To SGOLO2, labelled common /COMM2/. C C USES : C Calls SCLCC2 to calculate the values of SCJ. C C C ------------------------------------------------------------ C C C This file defines the common block /COMM2/ and some C associated parameters. These are for use in the base 2 C version of the generator. C C C The parameter MAXDIM is the maximum dimension that will be used. C NBITS is the number of bits (not counting the sign) in a C fixed-point integer. C C C The common block /COMM2/ : C SCJ - Contains the packed values of User choice of C Niederreiter's C(I,J,R) C DIMEN - The dimension of the sequence to be generated C COUNT - The index of the current item in the sequence, C expressed as an array of bits. COUNT(R) is the same C as Niederreiter's AR(N) (page 54) except that N is C implicit. C NEXTQ - The numerators of the next item in the series. These C are like Niederreiter's XI(N) (page 54) except that C N is implicit, and the NEXTQ are integers. To obtain C the values of XI(N), multiply by RECIP (see GOLO2). C C Array SCJ of the common block is set up by subroutine SCLCC2. C The other items in the common block are set up by INLO2. C C ------------------------------------------------------------ C C C C C .. Parameters .. INTEGER MAXDIM,NBITS PARAMETER (MAXDIM=318,NBITS=31) C .. C .. Scalar Arguments .. INTEGER DIM,IFLAG,SKIP C .. C .. Scalars in Common .. INTEGER COUNT,DIMEN C .. C .. Arrays in Common .. INTEGER NEXTQ(MAXDIM),SCJ(MAXDIM,0:NBITS-1) C .. C .. Local Scalars .. INTEGER GRAY,I,R C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. External Subroutines .. EXTERNAL SCLCC2 C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /COMM2/SCJ,DIMEN,COUNT,NEXTQ C .. C .. Save statement .. SAVE /COMM2/ C .. DIMEN = DIM C C This assignment just relabels the variable for C subsequent use. C C IF ((DIMEN .LE. 0) .OR. C * (DIMEN .GT. MAXDIM)) THEN C WRITE (*,*) ' SINLO2 : Bad dimension' C STOP C ENDIF C CALL SCLCC2(IFLAG) C C Translate SKIP into Gray code C GRAY = EXOR(SKIP,SKIP/2) C C Now set up NEXTQ appropriately for this value of the Gray code C C R = 0 10 IF (GRAY.NE.0) THEN IF (MOD(GRAY,2).NE.0) THEN DO 20 I = 1,DIMEN NEXTQ(I) = EXOR(NEXTQ(I),SCJ(I,R)) 20 CONTINUE END IF GRAY = GRAY/2 R = R + 1 GO TO 10 END IF C COUNT = SKIP RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Sniedx' then mkdir 'Sniedx' fi cd 'Sniedx' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' include makefile.inc all: Res1 SRCLIBS= $(LINPACK) $(EISPACK_DP) $(BLAS) $(PORT) src.o: $(INCSRC)/src.f $(F77) $(F77OPTS) -c $(INCSRC)/src.f extra.o: $(INCSRC)/extra.f $(F77) $(F77OPTS) -c $(INCSRC)/extra.f DRIVERS= driver RESULTS= Res1 SRCLIBS = $(TIMER) Objs1= driver.o src.o extra.o driver: $(Objs1) $(F77LINK) $(F77LINKOPTS) -o driver $(Objs1) $(SRCLIBS) Res1: driver ./driver >Res1 diffres: Res1 res1 echo "Differences in results from driver" $(DIFF) Res1 res1 clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' PROGRAM TNIEDX C C C THIS PROGRAM TESTS ACCURACY OF C NUMERICAL INTEGRATION USING "GNNIEDX" C AND INTEGRAND (2) OF DAVIS AND C RABINOWITZ, PAGE 406 C C User Define: C DIMEN : dimension C ATMOST : sequence length C SAMS : Number of replications C MAXS : Maximum Digits of Scrambling Of Owen type Scrambling C IFLAG: User Choice of type Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,II,IFLAG,K,J,MAXS,SAM C .. C .. Local Arrays .. DOUBLE PRECISION OUTS(40,10000) REAL SECOND, T1 C .. C .. C .. External Subroutines .. EXTERNAL GNNEDX,TIME C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. SAM = 2 MAXS = 31 DIMEN = 2 IFLAG = 1 ATMOST = 2**12 DO 30 II = 1,SAM WRITE (*,FMT=*) 'I = ITERATION NUMBER' WRITE (*,FMT=*) 'EI = ESTIMATED INTEGRAL' T1 = SECOND() CALL GNNEDX(DIMEN,ATMOST,IFLAG,MAXS,OUTS) SUM = 0.0 K = 9 DO 20 I = 1,ATMOST F = 1.0 DO 10 J = 1,DIMEN F = F*ABS(4.0*OUTS(J,I)-2.0) 10 CONTINUE IF (MOD(I,2**K).EQ.0) THEN K = K + 1 WRITE (*,FMT=*) 'I = ',I WRITE (*,FMT=*) 'EI = ',SUM/I END IF SUM = SUM + F 20 CONTINUE T1 = SECOND() - T1 WRITE (*,FMT=*) 'Total time elapsed = ',T1 30 CONTINUE STOP END SHAR_EOF fi # end of overwriting check if test -f 'nxs10m30' then echo shar: will not over-write existing file "'nxs10m30'" else cat << "SHAR_EOF" > 'nxs10m30' 601100350 67680573 1003196763 53112190 492752882 781048748 1069541657 224897837 26974858 517275145 451634153 30468723 313436293 752607609 187954075 724531163 999249265 40633782 478349820 293281792 216175174 144616791 33199907 130589174 80853660 13346823 44263879 57650027 26225812 21324732 687408335 1030296690 175105697 700504437 756791551 4206536 49860717 833303642 713701768 216343475 632636211 347953333 907692016 71399327 824918801 830143602 449297368 970005117 1018642708 577263874 751927066 258968587 123507468 478028610 203228312 177582754 155992224 54811421 19937300 113152282 921329673 649217276 727081678 547133168 767606725 631686577 710872245 866003693 37583364 635204116 204885920 478116368 164232231 199943409 615983887 989479985 711448216 841924349 271236280 833057230 25760385 389611716 580794291 44232556 920750979 962104500 185550548 730922590 1054857325 836052499 28934698 1068535527 755323706 526909559 804493913 1014525536 764397968 619975748 948657094 429469934 362457301 916436499 568092370 848255946 160246088 765495053 685927399 188505087 300899219 418292680 958128278 886337310 387653334 555779 1008436211 568562020 991258129 749968623 942282279 253779389 375007715 124858755 414730981 48823893 651890047 409231005 539706186 687786607 41086666 627879067 514503092 788547937 107586243 447614530 386025311 103422703 881307080 730530903 861344152 989276983 404865006 393338132 126108916 664019830 718337097 253850012 936179030 113777181 857262754 958455437 712826990 946979290 193472891 401936841 91455682 840790241 438488014 467139983 505031665 894492005 1042693246 1009011283 394533495 93427647 486374176 398103100 880362108 1073355548 512954968 844823800 352879549 108809984 217619968 962645870 252262189 1057997208 504537080 426089156 667552385 701857563 423829955 741237157 460497274 858536270 550366782 107064721 595468149 566030992 246604953 312884224 88913615 758861086 781240458 417558862 177460770 549410222 210428889 90596570 233828267 487312547 469458463 238901236 501346906 107338005 783050872 757526602 184581810 574671503 1014985482 990799856 611844349 137738448 519272294 475815320 673907190 100083641 568553807 718554306 956710208 732127731 1022304853 537564689 612416799 36505570 94186904 125473587 410062645 963718952 554345898 138179562 567629417 592623580 638067409 902136638 291608000 3517876 130676446 936083264 1011580017 573503066 227878475 958446878 593604780 495026943 915398759 743517565 507335141 102089007 306996234 419433675 693311399 679609836 314930875 338971979 537304327 462280239 113432193 909480871 284588106 701962134 1013431625 497903278 99876062 346028649 337394225 966537213 871038360 520762684 789290483 740734599 811972729 956690834 52357130 789379354 597858417 1034654136 780511120 326550583 478918767 98795812 456478730 258162657 19360625 683666161 805999686 233824978 440505476 875932528 662991655 320940586 64417554 252603479 118324798 237418699 264231894 369726699 380463530 386299779 834720337 555983499 SHAR_EOF fi # end of overwriting check if test -f 'nxs11m30' then echo shar: will not over-write existing file "'nxs11m30'" else cat << "SHAR_EOF" > 'nxs11m30' 683644482 135580830 114750887 803113694 903684270 72157854 59862132 505680715 779216195 193586105 329775408 1049470377 790869458 529177858 770596904 1072006936 140712293 266923202 844524433 114947652 963169856 269847095 96724448 724232372 729133779 872385703 803906910 638059470 526484481 843345462 431170622 1057600395 413165996 637333254 738214251 987594285 394318066 718140016 1064908418 636856468 856354958 404426257 6726334 916596837 284819911 905844483 331078851 514311810 280380871 612727606 184724875 578053121 197887820 50641495 423123265 227871988 36070792 465200006 187446371 99681590 983367218 414604440 243776932 990984772 289169669 72196638 522106132 18453143 1008259046 521211985 88531465 816692680 1001849242 432014653 124800720 554138134 166704786 634462798 1019392879 691725817 459507961 1043416360 536793729 1046579554 403993816 412804346 138344837 24248109 184604157 535890993 869278100 125797533 322976214 161488107 755812753 541238672 841393456 909350792 521628897 567331204 118313388 265208159 468752989 169896169 959104945 500814545 434091372 113024960 444725900 426474642 263586462 432120698 45230250 524531296 327681208 288240183 353741403 192041686 42986304 258035065 388012188 537416171 292523858 794840680 502322834 423539584 884204171 419181022 618098449 257684411 288733678 145066288 756812965 357038497 881037192 314634603 182012383 753118967 636654128 318546833 71382991 400554255 153755763 187747483 306343731 1044780 402486338 564947470 934172722 813806089 705671218 795868370 926341726 51211860 739601563 176287958 1016675820 950546833 267929114 275671035 666917527 431607592 177069271 570312219 125215260 886835363 75079442 767034988 879903592 525306349 560816857 973795128 774079162 143751759 137446578 209360240 1012892947 398130103 838233295 47600869 312065869 24489594 686091463 1064468817 610953508 593516621 736434891 331884657 295796744 991620199 197203105 694133440 834981040 813915471 1043396296 246522553 741863988 188767345 707349375 236730174 943797129 1073073438 917394477 81723820 598848363 17981413 322240492 766981617 282924362 643850936 914020733 586984616 128906358 739978294 543175409 351896331 1017052692 802102481 946877967 615436223 305014361 1038674293 372513724 902652345 552583238 507262995 36761397 753350193 116032862 328670174 103551339 1045088961 667010209 737933140 945361853 324875181 353779021 866388681 255506551 65493883 220170356 837239777 771797016 542598812 568193573 359763920 870736829 304824104 430870768 116534859 741311666 157845915 861348937 985097428 714048271 306167010 709889598 498574308 382576036 360956609 1005272948 680195561 845547054 1026217890 374002007 866217449 695335995 462692419 314161108 1048195907 709323042 604523286 581285813 463752482 919585493 518104348 905143615 884163267 99734344 332656669 398822849 973544107 459181032 1001255433 846689598 32962550 6372759 727573859 495342420 738892195 834835824 1047729137 529022179 1032177776 185764687 483299066 818974724 826550698 987239779 684232759 940422073 224992310 164068454 827234589 690302353 32197089 628566875 95166490 866947171 476745989 604118539 121738255 667653407 987681310 749263282 863528381 341372968 244006930 951984484 614173490 685305366 592527405 383027710 952965461 21271273 653806680 745437736 180332506 331728868 207307411 SHAR_EOF fi # end of overwriting check if test -f 'nxs12m30' then echo shar: will not over-write existing file "'nxs12m30'" else cat << "SHAR_EOF" > 'nxs12m30' 455734274 770944497 251274392 943939844 634593383 142647816 518405655 1035112060 686290424 65615251 693694295 260630152 43423443 881037832 267195896 910495556 605034094 965336736 151131045 582450365 288640393 124889081 201348992 232395086 964441802 612863572 1058289724 400689983 769068112 870006638 442034405 43938580 530724371 1060644604 1057605153 312158786 847892251 1022077706 18501972 383515470 371298986 47901519 323627800 867544471 192786854 203960857 720589590 545513545 464140138 547595136 376987492 501322404 192506954 737887979 949009184 729312279 192921452 991695031 102006938 776205479 388727334 696774466 350675311 1063266075 378755246 1019011080 540197062 820200324 366443609 41446382 903249144 329275036 1052721345 498980004 58772674 999024019 969764777 1051794123 99786007 727840180 172373356 204871234 384889035 806403840 53695037 13616293 66161898 186684471 98164394 355219262 217386272 693439660 965615244 154762973 436447974 799925933 668492905 842607298 791036533 626067281 22239456 678586727 806410990 935957395 13378684 139417906 542947431 866282064 955470920 476415717 29319813 738723019 310863222 837714036 747917136 328347967 851029220 767666096 964745988 401734446 124426506 534017541 504807159 389066462 356357537 725316735 659277475 756136726 958461663 655058111 114257772 133795835 1023812401 378615015 543820138 371991839 810602804 101893785 247079695 198575477 217168727 930751227 875740551 1021343329 521456293 239470050 84958534 620743129 123271774 544238382 669572624 333843660 444346220 356140799 190081678 512353278 796921871 600415539 453475202 848657261 603109415 745283747 151883079 1060772111 714582242 271029343 735760924 35818366 208189017 742934659 723194823 754382856 220628662 157700063 118306598 536381889 565112788 101605440 100486671 530359790 719846920 52922347 936119517 608310191 510418084 1033430645 567827140 661160957 258328750 718946670 236310941 133803377 168935681 1043464914 987196881 459837085 823469307 851754223 184441250 596703707 290780171 59741935 815598374 416271082 334275066 488096772 482629079 1019633822 257691498 612146801 915787995 171742109 134941161 748302528 401228388 175745030 524781059 599865139 22024809 736041817 37791189 53443125 1022208028 618558541 973936607 680773347 231100310 188099243 387492997 997909507 712762618 864453370 971076324 422303104 887951713 70908379 321106277 37257531 377525304 993066805 541243164 109743129 792736603 187956981 26958383 853015617 329470240 172383273 874347828 235666399 427079847 514437634 633169581 837815932 208607260 888137353 588084793 689649957 913247571 957569656 82894851 272748879 517167675 597600805 14301272 421502435 535550497 489600072 588119866 486064591 229677813 200942584 1056966489 703305147 118583330 634757241 367117372 850863584 691542163 625261051 975595510 687522116 186898165 247234909 363416034 1035231735 75052476 868760869 950694368 433949613 253713205 898716701 270518328 164766745 50123223 106155344 691291959 83686480 349811315 867635353 762318664 1073642543 564967196 486146330 994152282 480436438 40635107 901402239 767712015 403085167 913686502 676736454 283632149 207423788 1032456183 556764953 573432388 539147584 1053096806 246348695 866546419 318131241 409479478 19314747 142040277 669577250 670828194 766440062 1037065698 815986243 1068147734 982381558 190909672 227921657 374739209 345354308 301243259 490874673 152439690 424307886 256544451 446409440 88939318 927493863 540479486 605696652 105999616 700777601 839609487 642661344 269620503 649934339 297064636 632059693 671701847 257012421 270575278 928304922 853023247 878153656 SHAR_EOF fi # end of overwriting check if test -f 'nxs13m30' then echo shar: will not over-write existing file "'nxs13m30'" else cat << "SHAR_EOF" > 'nxs13m30' 455734274 770944497 251274392 943939844 634593383 142647816 518405655 1035112060 686290424 65615251 693694295 260630152 43423443 881037832 267195896 910495556 605034094 965336736 151131045 582450365 288640393 124889081 201348992 232395086 964441802 612863572 1058289724 400689983 769068112 870006638 442034405 43938580 530724371 1060644604 1057605153 312158786 847892251 1022077706 18501972 383515470 371298986 47901519 323627800 867544471 192786854 203960857 720589590 545513545 464140138 547595136 376987492 501322404 192506954 737887979 949009184 729312279 192921452 991695031 102006938 776205479 388727334 696774466 350675311 1063266075 378755246 1019011080 540197062 820200324 366443609 41446382 903249144 329275036 1052721345 498980004 58772674 999024019 969764777 1051794123 99786007 727840180 172373356 204871234 384889035 806403840 53695037 13616293 66161898 186684471 98164394 355219262 217386272 693439660 965615244 154762973 436447974 799925933 668492905 842607298 791036533 626067281 22239456 678586727 806410990 935957395 13378684 139417906 542947431 866282064 955470920 476415717 29319813 738723019 310863222 837714036 747917136 328347967 851029220 767666096 964745988 401734446 124426506 534017541 504807159 389066462 356357537 725316735 659277475 756136726 958461663 655058111 114257772 133795835 1023812401 378615015 543820138 371991839 810602804 101893785 247079695 198575477 217168727 930751227 875740551 1021343329 521456293 239470050 84958534 620743129 123271774 544238382 669572624 333843660 444346220 356140799 190081678 512353278 796921871 600415539 453475202 848657261 603109415 745283747 151883079 1060772111 714582242 271029343 735760924 35818366 208189017 742934659 723194823 754382856 220628662 157700063 118306598 536381889 565112788 101605440 100486671 530359790 719846920 52922347 936119517 608310191 510418084 1033430645 567827140 661160957 258328750 718946670 236310941 133803377 168935681 1043464914 987196881 459837085 823469307 851754223 184441250 596703707 290780171 59741935 815598374 416271082 334275066 488096772 482629079 1019633822 257691498 612146801 915787995 171742109 134941161 748302528 401228388 175745030 524781059 599865139 22024809 736041817 37791189 53443125 1022208028 618558541 973936607 680773347 231100310 188099243 387492997 997909507 712762618 864453370 971076324 422303104 887951713 70908379 321106277 37257531 377525304 993066805 541243164 109743129 792736603 187956981 26958383 853015617 329470240 172383273 874347828 235666399 427079847 514437634 633169581 837815932 208607260 888137353 588084793 689649957 913247571 957569656 82894851 272748879 517167675 597600805 14301272 421502435 535550497 489600072 588119866 486064591 229677813 200942584 1056966489 703305147 118583330 634757241 367117372 850863584 691542163 625261051 975595510 687522116 186898165 247234909 363416034 1035231735 75052476 868760869 950694368 433949613 253713205 898716701 270518328 164766745 50123223 106155344 691291959 83686480 349811315 867635353 762318664 1073642543 564967196 486146330 994152282 480436438 40635107 901402239 767712015 403085167 913686502 676736454 283632149 207423788 1032456183 556764953 573432388 539147584 1053096806 246348695 866546419 318131241 409479478 19314747 142040277 669577250 670828194 766440062 1037065698 815986243 1068147734 982381558 190909672 227921657 374739209 345354308 301243259 490874673 152439690 424307886 256544451 446409440 88939318 927493863 540479486 605696652 105999616 700777601 839609487 642661344 269620503 649934339 297064636 632059693 671701847 257012421 270575278 928304922 853023247 878153656 908127231 456920491 575103776 1014645311 649190824 471861745 303688659 191446425 1059703468 735083227 1073560050 399977794 353339158 897373374 901344073 326945933 972224406 541542931 797731354 787796494 617219174 750951417 1027765865 918702517 147665819 56969400 191051689 648505617 1020627885 515999460 SHAR_EOF fi # end of overwriting check if test -f 'nxs14m30' then echo shar: will not over-write existing file "'nxs14m30'" else cat << "SHAR_EOF" > 'nxs14m30' 363719563 126351411 333577928 1056266471 863617271 545956392 264984465 681943159 78085838 303365306 714461760 878096179 498252671 361921399 1063052908 70946244 104581208 524819162 118770356 796438038 248234408 966943067 684422354 949155117 14743351 682186614 276383621 954908515 908779080 188650849 726717595 594712286 189627205 962626877 339264590 919793712 551920957 1071032665 966676047 389380967 907398637 1006152404 844213372 488904161 922394447 158553023 92371877 636183969 966656831 508946316 1007314493 706515892 896386477 22341291 841210574 287818538 583447977 19696092 215190261 728128467 409766165 80949732 958979209 201976339 297121774 214729864 1043245631 270423146 301241045 679963780 566269527 422403774 165302426 914777783 144300990 437527082 419837406 885130654 774249209 824807568 172485538 384520341 995290890 522922015 381287522 732993502 420132767 554262480 306479536 641326327 467561434 930123715 597698854 648058317 947215434 275322696 189119691 925373570 342298247 854593552 306693696 427297494 565611478 689496359 407839947 480181175 474306049 623831144 573116812 1025353306 646127332 20717872 692237780 790982509 780585624 520913570 542647738 114251809 908490564 323313502 598016715 36879633 413357082 688134786 529390887 836525405 101207742 490664799 11494833 522299314 551475773 292805888 497724752 588769302 730496039 386917419 875293549 309251161 897721384 773076175 599095939 647256214 285476502 674944667 56159762 276099032 568461366 844675922 37541153 455402722 540203051 758670719 741617654 37819314 140948332 49760762 21597375 526385708 1061539985 1045701947 603427637 786977027 254345255 1026578814 188180846 403679758 46785960 556911344 817653155 129910636 383529639 702519067 1008955264 816285592 389874521 1025108190 370052760 361843901 782388792 727300723 874780780 1032318137 786775357 593391885 441843415 878415332 42938588 326487870 130723815 191703994 387341724 771600788 710973497 154624564 994157100 978001805 434779620 1055151620 881172851 146590316 623010517 791412177 915915557 915210661 1072755931 510341731 93253334 795920433 574195950 611866662 990826533 203952029 870492277 1052220963 186404076 793778659 184746614 408220858 79873149 150093288 58150406 326849960 839453518 505919819 120084119 686127014 994898089 776213883 736633924 615787916 197990472 271596250 910133156 81863853 334711432 104742150 599466742 411643997 30005230 473224083 89878894 581365075 871966718 936443977 868585761 59510824 211223195 741257226 566789837 324012150 1000420904 390294924 559828109 547969020 899816526 44034976 402610207 879319035 751749594 750457816 896109678 200664812 356839039 377943788 75871134 988191444 612525298 428181965 79060962 245808923 659438054 976230761 621901013 990671526 779364542 1008101719 1048836861 753543499 531728282 630806886 407682815 527132550 731180892 196920828 663869657 1063558824 652317284 14622439 877620946 1008507639 904200452 71792432 893965205 336930819 744920405 884732703 593736443 383950873 713171048 212736844 251296924 866846431 882585823 726089887 570430933 124674441 1014322098 191282236 401072059 257700805 623630253 721097945 314891164 733331637 161550360 502391973 816418068 184735944 880183527 833093731 766925329 862523522 766437102 314856969 693813092 577743131 6892938 693316693 513249816 715150528 902131452 536908016 703024026 255588483 797459039 642799176 125289328 34940810 814984727 154130419 954568665 364941906 776843713 523478561 621488164 344765525 209934804 296960280 4686125 323886097 274739907 52824308 374139161 216143239 1040631726 250585515 1072159745 1070231107 295878124 356423628 707647585 668470041 846895146 890530169 397463630 310838951 249187095 826097215 49220519 335241096 538058126 720431917 553546564 475440125 359663258 627836411 525395389 479379746 139983909 166538312 767722472 307846915 155931062 286500052 986356837 394275280 135296169 815275901 202685428 647738757 531267531 63904746 194425588 529308760 139050208 550214416 841938907 851674637 474033192 137868193 934809352 867254475 714397996 188211231 53025880 906031050 966318476 6712754 406760659 922452103 522341224 570275994 81275582 926019935 447282841 562574040 394653639 794947205 668724913 686109321 379752491 650190375 323192837 734194683 38112520 58093839 389562127 875285157 787936082 760183554 29226022 532719339 344312085 142617226 588163644 105426877 177250683 162868885 574641189 7669603 193724285 185612450 217098721 135392219 1033477436 821412987 847784694 641540976 35355503 447014323 SHAR_EOF fi # end of overwriting check if test -f 'nxs15m30' then echo shar: will not over-write existing file "'nxs15m30'" else cat << "SHAR_EOF" > 'nxs15m30' 363719563 126351411 333577928 1056266471 863617271 545956392 264984465 681943159 78085838 303365306 714461760 878096179 498252671 361921399 1063052908 70946244 104581208 524819162 118770356 796438038 248234408 966943067 684422354 949155117 14743351 682186614 276383621 954908515 908779080 188650849 726717595 594712286 189627205 962626877 339264590 919793712 551920957 1071032665 966676047 389380967 907398637 1006152404 844213372 488904161 922394447 158553023 92371877 636183969 966656831 508946316 1007314493 706515892 896386477 22341291 841210574 287818538 583447977 19696092 215190261 728128467 409766165 80949732 958979209 201976339 297121774 214729864 1043245631 270423146 301241045 679963780 566269527 422403774 165302426 914777783 144300990 437527082 419837406 885130654 774249209 824807568 172485538 384520341 995290890 522922015 381287522 732993502 420132767 554262480 306479536 641326327 467561434 930123715 597698854 648058317 947215434 275322696 189119691 925373570 342298247 854593552 306693696 427297494 565611478 689496359 407839947 480181175 474306049 623831144 573116812 1025353306 646127332 20717872 692237780 790982509 780585624 520913570 542647738 114251809 908490564 323313502 598016715 36879633 413357082 688134786 529390887 836525405 101207742 490664799 11494833 522299314 551475773 292805888 497724752 588769302 730496039 386917419 875293549 309251161 897721384 773076175 599095939 647256214 285476502 674944667 56159762 276099032 568461366 844675922 37541153 455402722 540203051 758670719 741617654 37819314 140948332 49760762 21597375 526385708 1061539985 1045701947 603427637 786977027 254345255 1026578814 188180846 403679758 46785960 556911344 817653155 129910636 383529639 702519067 1008955264 816285592 389874521 1025108190 370052760 361843901 782388792 727300723 874780780 1032318137 786775357 593391885 441843415 878415332 42938588 326487870 130723815 191703994 387341724 771600788 710973497 154624564 994157100 978001805 434779620 1055151620 881172851 146590316 623010517 791412177 915915557 915210661 1072755931 510341731 93253334 795920433 574195950 611866662 990826533 203952029 870492277 1052220963 186404076 793778659 184746614 408220858 79873149 150093288 58150406 326849960 839453518 505919819 120084119 686127014 994898089 776213883 736633924 615787916 197990472 271596250 910133156 81863853 334711432 104742150 599466742 411643997 30005230 473224083 89878894 581365075 871966718 936443977 868585761 59510824 211223195 741257226 566789837 324012150 1000420904 390294924 559828109 547969020 899816526 44034976 402610207 879319035 751749594 750457816 896109678 200664812 356839039 377943788 75871134 988191444 612525298 428181965 79060962 245808923 659438054 976230761 621901013 990671526 779364542 1008101719 1048836861 753543499 531728282 630806886 407682815 527132550 731180892 196920828 663869657 1063558824 652317284 14622439 877620946 1008507639 904200452 71792432 893965205 336930819 744920405 884732703 593736443 383950873 713171048 212736844 251296924 866846431 882585823 726089887 570430933 124674441 1014322098 191282236 401072059 257700805 623630253 721097945 314891164 733331637 161550360 502391973 816418068 184735944 880183527 833093731 766925329 862523522 766437102 314856969 693813092 577743131 6892938 693316693 513249816 715150528 902131452 536908016 703024026 255588483 797459039 642799176 125289328 34940810 814984727 154130419 954568665 364941906 776843713 523478561 621488164 344765525 209934804 296960280 4686125 323886097 274739907 52824308 374139161 216143239 1040631726 250585515 1072159745 1070231107 295878124 356423628 707647585 668470041 846895146 890530169 397463630 310838951 249187095 826097215 49220519 335241096 538058126 720431917 553546564 475440125 359663258 627836411 525395389 479379746 139983909 166538312 767722472 307846915 155931062 286500052 986356837 394275280 135296169 815275901 202685428 647738757 531267531 63904746 194425588 529308760 139050208 550214416 841938907 851674637 474033192 137868193 934809352 867254475 714397996 188211231 53025880 906031050 966318476 6712754 406760659 922452103 522341224 570275994 81275582 926019935 447282841 562574040 394653639 794947205 668724913 686109321 379752491 650190375 323192837 734194683 38112520 58093839 389562127 875285157 787936082 760183554 29226022 532719339 344312085 142617226 588163644 105426877 177250683 162868885 574641189 7669603 193724285 185612450 217098721 135392219 1033477436 821412987 847784694 641540976 35355503 447014323 SHAR_EOF fi # end of overwriting check if test -f 'nxs16m30' then echo shar: will not over-write existing file "'nxs16m30'" else cat << "SHAR_EOF" > 'nxs16m30' 363719563 126351411 333577928 1056266471 863617271 545956392 264984465 681943159 52463357 303365306 341871271 878096179 1027086679 361921399 39103888 62748151 338638050 524819162 3876208 267521086 458865375 966943067 802258923 713427351 531164653 682186614 37351620 762889236 260620051 188650849 726717595 594712286 189627205 962626877 339264590 919793712 551920957 1071032665 966676047 1054115635 907398637 1006152404 162072233 344789835 922394440 158553011 92371900 873005565 640169584 237806147 641542344 519118056 896386491 536714453 859608546 634014818 618417251 3393265 156430672 666226982 409766165 80949732 958979209 201976339 297121774 214729864 1043245631 270423146 301241045 679963780 566269527 955107561 284626468 914777783 144300991 437527080 419837403 442171142 965324889 1065357291 410833689 561063529 575986391 237595697 210154797 521379252 53870289 554262484 709618621 927211987 467561434 930123715 597698854 648058317 947215434 275322696 189119691 925373570 342298247 854593552 396346066 427297494 1024886163 336487607 407839949 480181179 474306072 65831972 730212082 687327018 644792536 1051908400 889709511 801681759 1023408433 855329054 1008563125 341106972 941232885 939036962 598016715 36879633 413357082 688134786 529390887 836525405 101207742 490664799 11494833 522299314 583110352 292805888 247866556 1052591429 730496033 386917415 875293557 99782499 920602639 425898890 156162358 107809551 623898110 46232889 637144328 372203353 365552457 910212709 280354127 260652191 874780780 1032318137 786775357 593391885 441843415 878415332 42938588 326487870 130723815 785144252 387341724 771600788 127740332 592968206 994157099 978001793 434779645 729282287 354410120 791171421 491923658 351481648 915915554 241023036 649736868 632545424 469673852 909774405 1003732247 531359009 540203051 758670719 741617654 37819314 140948332 49760762 21597375 526385708 1061539985 579527250 603427637 786977027 37143057 1026578814 188180847 403679756 46785965 667156179 891312700 471918364 294250646 432700999 1008955269 938630575 36531043 218890627 413601250 8463491 782388798 91934927 990826533 203952029 870492277 1052220963 186404076 793778659 184746614 408220858 79873149 150093288 58150406 269228974 570757344 505919819 120084118 686127012 994898093 298863363 484781311 687199643 650162558 61321752 225670925 285598097 274126380 759224780 538885718 411643995 571217690 114403460 89878894 581365075 871966718 936443977 868585761 59510824 211223195 741257226 566789837 324012150 207363495 390294924 559828109 377003347 877987647 44034976 402610207 448122428 925503731 282184435 603842202 427105648 29326592 862664874 1062493067 275667799 802501187 198359584 67791473 940026761 659438054 976230761 621901013 990671526 779364542 1008101719 1048836861 753543499 877275007 241367427 289394386 527132550 731180892 196920828 741154085 1063558824 73339515 671306328 177839438 280603233 69598435 553636047 21758641 679392417 1040597771 639040741 867477365 195068519 904480750 430967557 251296924 866846431 882585823 726089887 570430933 124674441 1014322098 191282236 584440131 257700805 623630253 199904397 422483729 733331637 818940721 502391973 600211113 832718859 880183527 105319506 911207421 75898035 546969523 569282401 997955765 943189515 6892938 392476249 519515299 878147683 902131452 536908016 703024026 255588483 797459039 642799176 125289328 34940810 814984727 154130419 954568665 757407115 380124696 523478561 343228356 344765525 209934804 733893040 4686125 157645789 1063363845 420310840 421343735 843471422 128937584 38733292 181218315 860436375 398555761 534645125 403099131 1036462677 892066497 341684534 223489259 234965009 19905077 931842084 1026230795 1011059645 746733161 358619256 432009135 983738138 24518058 297205887 347292813 58075623 134401733 539020927 770756513 301173416 363639933 789096518 418407949 242383984 662950093 1055777508 79433174 640498181 707647585 668470041 846895146 890530169 397463630 310838951 249187095 826097215 49220519 335241096 782929433 720431917 553546564 591356063 380704065 627836411 525395389 795160220 747434339 134050143 928864726 711080970 160565023 698130484 346507239 762040514 929275449 425334989 968144228 422912892 379752491 650190375 323192837 734194683 38112520 58093839 389562127 875285157 787936082 760183554 29226022 511529677 356007219 142617226 1052163461 105426877 177250683 265775255 574641189 509636365 109461310 352669900 797298167 858221260 605131882 721505885 563118811 749568523 926871602 621924245 531267531 63904746 194425588 529308760 139050208 550214416 841938907 851674637 967506763 137868193 934809352 1056799343 421894631 188211231 289105312 906031050 857200977 369465542 406760659 409343592 23228780 257961589 866153129 766964599 636314676 365872358 394653639 553244867 894090625 255777848 SHAR_EOF fi # end of overwriting check if test -f 'nxs4m30' then echo shar: will not over-write existing file "'nxs4m30'" else cat << "SHAR_EOF" > 'nxs4m30' 982024320 536870912 1027883072 268435456 513941536 134217728 256970768 67108864 128485384 33554432 64242692 16777216 32121346 8388608 16060673 4194304 8030336 2097152 4015168 1048576 2007584 524288 1003792 262144 501896 131072 250948 65536 125474 32768 890306817 706871810 1046873606 680003592 134094366 518794284 447507564 142639232 985003508 610986728 893304504 581478560 189187800 333241584 329030896 8421376 904193297 717802018 1054400102 673196168 121186302 506601196 439769772 134776960 978303668 613311592 898315576 571122336 193809240 319804400 753991424 144747008 817941504 545294336 123770368 977331200 201375744 134250496 641686528 69888000 1070346240 713564160 253106176 904138752 12582912 8388608 744230911 141036202 808455372 538970248 132642454 985310436 214745280 143163520 647586652 82723432 1060371504 706914336 261688856 890369040 319750399 571080874 1007419596 671613064 301334678 514064612 213909696 142606464 90439516 105032296 66898992 44599328 17733144 32826384 12632064 8421376 5463055 6449162 3944460 2629640 1173494 2013860 838848 559232 353948 410856 261936 174624 69528 128272 SHAR_EOF fi # end of overwriting check if test -f 'nxs4m30.old' then echo shar: will not over-write existing file "'nxs4m30.old'" else cat << "SHAR_EOF" > 'nxs4m30.old' SHAR_EOF fi # end of overwriting check if test -f 'nxs5m30' then echo shar: will not over-write existing file "'nxs5m30'" else cat << "SHAR_EOF" > 'nxs5m30' 936807843 468403921 828984403 356582721 595421402 297710701 648208792 527539158 263769579 470128177 172395292 86197646 192816965 105493743 52746871 112025243 35045196 17522598 38715510 32079951 16039975 31293608 11436052 5718026 12329167 6467624 3233812 7083215 2317878 1158939 189424280 960093815 536870912 691658359 345829179 172914589 142740098 71370049 35685024 45744832 22872416 11436208 8421384 4210692 2105346 2691301 1345650 672825 559752 279876 139938 178360 89180 44590 32768 16384 8192 10553 5276 2638 378848561 757697122 594201266 168969361 957154112 186898679 393620873 642667409 900036949 815822848 930586077 399959500 148592212 804995055 651141544 447900379 983427706 201460354 246707548 606769380 829226943 607693803 872469364 288913566 393224777 755278824 590387622 186144468 959712835 189907697 686339538 1027381875 213116256 554325474 200129167 153587669 177289635 431207461 1054283608 406761944 315351808 155168968 577870202 802378588 446674946 753819964 361512087 277004059 34461209 249231638 738561728 376489647 181133977 393332398 676050087 1006672852 231152336 567934834 169205014 144990292 627373762 875882387 710212082 1041558963 283353473 315523601 258432656 98890848 601641341 852033833 147134524 706550037 817216966 660987536 936040432 584028852 155262085 66212886 528536171 306157518 961418166 692786855 393225133 1060844555 607693495 889246273 730652526 1071156015 293079512 333204864 SHAR_EOF fi # end of overwriting check if test -f 'nxs6m30' then echo shar: will not over-write existing file "'nxs6m30'" else cat << "SHAR_EOF" > 'nxs6m30' 606480775 465803779 181587458 568432641 871533762 290784517 90793729 536870912 615257700 516719748 592784257 268435456 962944305 954326085 559569344 94739461 568852125 453491749 279784672 134217728 438382924 183798805 99814005 67108864 11863203 262703112 137055032 33554432 176704595 131351556 853996292 526861827 728436739 11570176 325294091 493386514 782564896 605115473 170921120 551607552 425266499 933836653 384112693 525740423 1010005209 196499649 915852042 541723895 923324257 608453636 710752575 55412690 1023321049 392372758 775661454 139030565 1016207989 991651735 435488018 85197121 699255567 511050258 191574051 558978113 439001226 1004174746 138036802 821882480 828971145 554762513 41095779 398529653 271830963 703243775 663834619 579951279 54389364 721959544 301998112 603996224 866216640 1028990854 652422655 1043378552 705154807 719738807 846089353 9681424 700527139 98879553 281303814 985894407 950700034 85005313 860600331 371193886 1056048168 59563588 734736515 614763796 922903331 626957236 670320823 928085694 511862620 1063432009 773068452 352314605 480397148 71319552 364298516 948946717 1004560439 61684053 917725859 313931751 1014484092 6707055 741982243 548460804 640002308 1046696450 576271878 426455041 912785406 891996329 39800934 78835236 9150916 105674258 299118841 825150637 122507264 93169063 164308668 427645603 1024535069 846396261 747694451 423456867 736365516 429426358 580263466 922096900 414940055 932227189 1062437416 201658913 820188420 1072646548 169454328 999738199 908931378 682273312 499382810 710890364 810766637 646637595 486301003 377255267 262577752 247127853 231815664 26982614 3437264 24227907 34545098 17521923 56762569 54526153 3331460 7447619 3793274 4311084 4633376 322319 700549 353514 930320 954426 SHAR_EOF fi # end of overwriting check if test -f 'nxs7m30' then echo shar: will not over-write existing file "'nxs7m30'" else cat << "SHAR_EOF" > 'nxs7m30' 606480775 465803779 181587458 568432641 871533762 290784517 90793729 536870912 615257700 516719748 592784257 268435456 962944305 954326085 559569344 94739461 568852125 453491749 279784672 134217728 438382924 183798805 99814005 67108864 11863203 262703112 137055032 33554432 176704595 131351556 853996292 526861827 728436739 11570176 325294091 493386514 782564896 605115473 170921120 551607552 425266499 933836653 384112693 525740423 1010005209 196499649 915852042 541723895 923324257 608453636 710752575 55412690 1023321049 392372758 775661454 139030565 1016207989 991651735 435488018 85197121 699255567 511050258 191574051 558978113 439001226 1004174746 138036802 821882480 828971145 554762513 41095779 398529653 271830963 703243775 663834619 579951279 54389364 721959544 301998112 603996224 866216640 1028990854 652422655 1043378552 705154807 719738807 846089353 9681424 700527139 98879553 281303814 985894407 950700034 85005313 860600331 371193886 1056048168 59563588 734736515 614763796 922903331 626957236 670320823 928085694 511862620 1063432009 773068452 352314605 480397148 71319552 364298516 948946717 1004560439 61684053 917725859 313931751 1014484092 6707055 741982243 548460804 640002308 1046696450 576271878 426455041 912785406 891996329 39800934 78835236 9150916 105674258 299118841 825150637 122507264 93169063 164308668 427645603 1024535069 846396261 747694451 423456867 736365516 429426358 580263466 922096900 414940055 932227189 1062437416 201658913 820188420 1072646548 169454328 999738199 908931378 682273312 499382810 710890364 810766637 646637595 486301003 377255267 262577752 247127853 231815664 26982614 3437264 24227907 34545098 17521923 56762569 54526153 3331460 7447619 3793274 4311084 4633376 322319 700549 353514 930320 954426 388460551 931657735 421192443 621105157 1037106001 669513735 647982689 67387396 221483079 513879041 350950732 105063427 929604210 805699590 290999399 537133060 753493059 341801984 259503533 290668550 453613872 438666243 1068517918 568630272 1019013139 1024809985 239569137 1052841990 367304536 100663302 SHAR_EOF fi # end of overwriting check if test -f 'nxs8m30' then echo shar: will not over-write existing file "'nxs8m30'" else cat << "SHAR_EOF" > 'nxs8m30' 971165186 508331020 732674767 1046873606 87632649 918850564 203264367 559219973 483398024 1020005388 406998203 813695112 396379340 536870912 195423833 25709644 901129312 268435456 39206696 1050812448 450564656 134217728 646272410 525406224 225282328 67108864 757378251 262703112 112641164 33554432 386036107 312339215 707937815 43697717 967607410 570951840 484762168 771247179 943823356 695597499 468984043 516686852 322278341 242177886 986861169 398251091 283855535 164576526 697032771 880616056 798901250 493511676 1065294580 143130632 508870870 322823590 728708934 214657483 1001687952 691120683 731609491 619933475 112872027 945702653 208699895 384599718 1032831408 505693482 685312077 35212225 740868458 537397376 74445050 246491299 450793213 692604740 857703665 851904404 396914158 822478136 925665122 289259643 345391241 8388616 723813487 609897559 114017584 941232093 210567856 377228463 403431811 806863623 544899611 551191097 1031342207 888046758 25494524 924942566 700979533 474591531 770139276 374080644 617990995 793401641 986842673 550719552 324394601 1043405112 472893038 530334688 806545654 900193647 1069820069 142608392 82061105 74650251 184121264 187751457 147170845 51907337 599612282 744088995 896805021 233032466 205913572 1045309068 926717027 256109464 20067032 805127486 251664498 545259520 427529290 977131340 52090010 855955499 752626937 1045572096 302263456 767751821 837744084 737998121 1025014392 910459908 145474983 1041363978 504172586 232612596 712421995 782115630 649735037 278383021 469385363 397086490 537295335 679520898 129309802 889794450 344282194 1017227189 501838591 516197388 676419649 195482432 42046614 749820973 412209920 688554157 181370167 945770665 372476172 281343112 291909000 1045203972 531692812 1068918520 607143639 679295073 830151498 53677572 852307082 644957963 500220665 133106526 339376740 518270092 605778238 329885136 232120047 325946539 593330004 536873096 272495248 538599187 897129537 959045204 667190 139207123 323571381 658897106 536468757 390033626 625726741 910461956 189848208 144186599 918437726 130754967 535474913 70659206 477210486 114380203 703905949 623871774 393638376 35695240 591360090 325659536 730726691 517012203 110572048 680036488 508120792 986513555 623281929 496378704 1002367817 460715509 1002469322 397683064 350220086 1018491818 778980373 910984204 369842672 810128731 167405666 339609997 276968643 641303465 SHAR_EOF fi # end of overwriting check if test -f 'nxs9m30' then echo shar: will not over-write existing file "'nxs9m30'" else cat << "SHAR_EOF" > 'nxs9m30' 601100350 67680573 1003196763 53112190 492752882 781048748 1069541657 224897837 26974858 517275145 451634153 30468723 313436293 752607609 187954075 724531163 999249265 40633782 478349820 293281792 216175174 144616791 33199907 130589174 80853660 13346823 44263879 57650027 26225812 21324732 687408335 1030296690 175105697 700504437 756791551 4206536 49860717 833303642 713701768 216343475 632636211 347953333 907692016 71399327 824918801 830143602 449297368 970005117 1018642708 577263874 751927066 258968587 123507468 478028610 203228312 177582754 155992224 54811421 19937300 113152282 921329673 649217276 727081678 547133168 767606725 631686577 710872245 866003693 37583364 635204116 204885920 478116368 164232231 199943409 615983887 989479985 711448216 841924349 271236280 833057230 25760385 389611716 580794291 44232556 920750979 962104500 185550548 730922590 1054857325 836052499 28934698 1068535527 755323706 526909559 804493913 1014525536 764397968 619975748 948657094 429469934 362457301 916436499 568092370 848255946 160246088 765495053 685927399 188505087 300899219 418292680 958128278 886337310 387653334 555779 1008436211 568562020 991258129 749968623 942282279 253779389 375007715 124858755 414730981 48823893 651890047 409231005 539706186 687786607 41086666 627879067 514503092 788547937 107586243 447614530 386025311 103422703 881307080 730530903 861344152 989276983 404865006 393338132 126108916 664019830 718337097 253850012 936179030 113777181 857262754 958455437 712826990 946979290 193472891 401936841 91455682 840790241 438488014 467139983 505031665 894492005 1042693246 1009011283 394533495 93427647 486374176 398103100 880362108 1073355548 512954968 844823800 352879549 108809984 217619968 962645870 252262189 1057997208 504537080 426089156 667552385 701857563 423829955 741237157 460497274 858536270 550366782 107064721 595468149 566030992 246604953 312884224 88913615 758861086 781240458 417558862 177460770 549410222 210428889 90596570 233828267 487312547 469458463 238901236 501346906 107338005 783050872 757526602 184581810 574671503 1014985482 990799856 611844349 137738448 519272294 475815320 673907190 100083641 568553807 718554306 956710208 732127731 1022304853 537564689 612416799 36505570 94186904 125473587 410062645 963718952 554345898 138179562 567629417 592623580 638067409 902136638 291608000 3517876 130676446 936083264 1011580017 573503066 227878475 958446878 593604780 495026943 915398759 743517565 507335141 102089007 306996234 419433675 693311399 679609836 314930875 338971979 537304327 462280239 113432193 909480871 284588106 701962134 1013431625 497903278 99876062 346028649 337394225 966537213 871038360 520762684 789290483 740734599 SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' I = ITERATION NUMBER EI = ESTIMATED INTEGRAL I = 512 EI = 0.9999238584613173 I = 1024 EI = 0.9999746664717595 I = 2048 EI = 0.9996899431116691 I = 4096 EI = 0.9999899950893638 Total time elapsed = 0.05282024 I = ITERATION NUMBER EI = ESTIMATED INTEGRAL I = 512 EI = 0.9948311093778697 I = 1024 EI = 0.9981168231523748 I = 2048 EI = 0.9999985867572326 I = 4096 EI = 0.9998295181301297 Total time elapsed = 0.0166993 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'extra.f' then echo shar: will not over-write existing file "'extra.f'" else cat << "SHAR_EOF" > 'extra.f' INTEGER FUNCTION EXOR(IIN,JIN) C C THIS FUNCTION CALCULATES THE EXCLUSIVE-OR OF ITS C TWO INPUT PARAMETERS C C .. Scalar Arguments .. INTEGER IIN,JIN C .. C .. Local Scalars .. INTEGER I,J,K,L C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. I = IIN J = JIN K = 0 L = 1 C 10 IF (I.EQ.J) THEN EXOR = K RETURN END IF C C CHECK THE CURRENT RIGHT-HAND BITS OF I AND J. C IF THEY DIFFER, SET THE APPROPRIATE BIT OF K. C IF (MOD(I,2).NE.MOD(J,2)) K = K + L I = I/2 J = J/2 L = 2*L GO TO 10 END DOUBLE PRECISION FUNCTION UNI() * * Random number generator, adapted from F. James * "A Review of Random Number Generators" * Comp. Phys. Comm. 60(1990), pp. 329-344. * C .. Parameters .. DOUBLE PRECISION TWOM24 PARAMETER (TWOM24=1D0/16777216.0) C .. C .. Local Scalars .. DOUBLE PRECISION CARRY INTEGER I,J C .. C .. Local Arrays .. DOUBLE PRECISION SEEDS(24) C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Save statement .. SAVE I,J,CARRY,SEEDS C .. C .. Data statements .. DATA I,J,CARRY/24,10,0.0d0/ DATA SEEDS/0.8804418d0,0.2694365d0,0.0367681d0,0.4068699d0, +0.4554052d0, + 0.2880635d0,0.1463408d0,0.2390333d0,0.6407298d0, +0.1755283d0,0.7132940d0, + 0.4913043d0,0.2979918d0,0.1396858d0,0.3589528d0, +0.5254809d0,0.9857749d0, + 0.4612127d0,0.2196441d0,0.7848351d0,0.4096100d0, +0.9807353d0,0.2689915d0, + 0.5140357d0/ C .. UNI = SEEDS(I) - SEEDS(J) - CARRY IF (UNI.LT.0) THEN UNI = UNI + 1 CARRY = TWOM24 ELSE CARRY = 0 END IF SEEDS(I) = UNI I = 24 - MOD(25-I,24) J = 24 - MOD(25-J,24) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE GNCRML(MAXS,LSM,SHIFT) C GENERATING LOWER TRIAGULAR SCRMABLING MATRICES C AND SHIFT VECTORS. C INPUTS : C FROM INNEDX : MAXS C FROM BLOCK DATA "NIEDX" : S, MAXCOL C C OUTPUTS : C TO INNEDX : LSM, SHIFT C .. Scalar Arguments .. INTEGER MAXS C .. C .. Array Arguments .. INTEGER LSM(40,31),SHIFT(40) C .. C .. Local Scalars .. INTEGER I,J,L,LL,P,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Common blocks .. COMMON /NIEDX/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /NIEDX/ C .. DO 30 P = 1,S SHIFT(P) = 0 L = 1 DO 20 I = MAXS,1,-1 LSM(P,I) = 0 STEMP = MOD((INT(UNI()*1000.0)),2) SHIFT(P) = SHIFT(P) + STEMP*L L = 2*L LL = 1 DO 10 J = 30,1,-1 IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.LT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF LSM(P,I) = LSM(P,I) + TEMP*LL LL = 2*LL 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GNCRMU(USM,USHIFT) C GENERATING UPPER TRIANGULAR SCRMABLING MATRICES AND C SHIFT VECTORS. C INPUTS : C FROM BLOCK DATA "NIEDX" : S, MAXCOL, C C OUTPUTS : C TO INSOBL : USM, USHIFT C .. Array Arguments .. INTEGER USHIFT(30),USM(30,30) C .. C .. Local Scalars .. INTEGER I,J,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Common blocks .. COMMON /NIEDX/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. SAVE /NIEDX/ C .. DO 20 I = 1,MAXCOL STEMP = MOD((INT(UNI()*1000.0)),2) USHIFT(I) = STEMP DO 10 J = 1,MAXCOL IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.GT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF USM(I,J) = TEMP 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE GNNEDX(DIMEN,ATMOST,IFLAG,MAXS,OUTS) C C User Define: C DIMEN : dimension C ATMOST : sequence length C MAXS : Maximum Digits of Scrambling C IFLAG: User Choice of Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,IFLAG,J,MAXS C .. C .. Local Arrays .. DOUBLE PRECISION QUASI(40),OUTS(40,10000) LOGICAL FLAG(2) C .. External Subroutines .. EXTERNAL GONEDX,INNEDX C .. CALL INNEDX(FLAG,DIMEN,ATMOST,QUASI,MAXS,IFLAG) DO 10 J = 1,DIMEN OUTS(J,1) = QUASI(J) 10 CONTINUE DO 30 I = 2,ATMOST CALL GONEDX(QUASI) DO 20 J = 1,DIMEN OUTS(J,I) = QUASI(J) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GONEDX(QUASI) C C THIS SUBROUTINE GENERATES A NEW C QUASIRANDOM VECTOR WITH EACH CALL C C IT ADAPTS THE IDEAS OF ANTONOV AND SALEEV, C USSR COMPUT. MATHS. MATH. PHYS. 19 (1980), C 252 - 256 C C THE USER MUST CALL "INNEDX" BEFORE CALLING C "GONEDX". AFTER CALLING "INNEDX", TEST C FLAG(1) AND FLAG(2); IF EITHER IS FALSE, C DO NOT CALL "GONEDX". "GOGONEDXCHECKS C THAT THE USER DOES NOT MAKE MORE CALLS C THAN HE SAID HE WOULD : SEE THE COMMENTS C TO "INNEDX". C C INPUTS: C FROM USER'S CALLING PROGRAM: C NONE C C FROM LABELLED COMMON /NIEDX/: C SV TABLE OF DIRECTION NUMBERS C S DIMENSION C MAXCOL LAST COLUMN OF V TO BE USED C COUNT SEQUENCE NUMBER OF THIS CALL C LASTQ NUMERATORS FOR LAST VECTOR GENERATED C RECIPD (1/DENOMINATOR) FOR THESE NUMERATORS C C C FIND THE POSITION OF THE RIGHT-HAND ZERO IN COUNT C C .. Array Arguments .. DOUBLE PRECISION QUASI(40) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Local Scalars .. INTEGER I,L C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /NIEDX/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /NIEDX/ C .. L = 0 I = COUNT 10 L = L + 1 IF (MOD(I,2).EQ.1) THEN I = I/2 GO TO 10 END IF C C CHECK THAT THE USER IS NOT CHEATING ! C IF (L.GT.MAXCOL) STOP ' TOO MANY CALLS ON GOSOBL' C C CALCULATE THE NEW COMPONENTS OF QUASI, C FIRST THE NUMERATORS, THEN NORMALIZED C DO 20 I = 1,S LASTQ(I) = EXOR(LASTQ(I),SV(I,L)) C C IF A FULL-WORD EXCLUSIVE-OR, SAY .XOR., IS AVAILABLE C THEN REPLACE THE PRECEDING STATEMENT BY C C LASTQ(I) = LASTQ(I) .XOR. SV(I,L) C C TO GET A FASTER, EXTENDED FORTRAN PROGRAM C QUASI(I) = LASTQ(I)*RECIPD 20 CONTINUE C COUNT = COUNT + 1 C RETURN END SUBROUTINE INNEDX(FLAG,DIMEN,ATMOST,QUASI,MAXS,IFLAG) C C THIS INITIALIZE NIED-XING SEQUENCE. C FIRST CHECK WHETHER THE USER-SUPPLIED C DIMENSION "DIMEN" OF THE QUASI-RANDOM C VECTORS IS STRICTLY BETWEEN 1 AND 16. C IF SO, FLAG(1) = .TRUE. C C NEXT CHECK "ATMOST", AN UPPER BOUND ON THE NUMBER C OF CALLS THE USER INTENDS TO MAKE ON "GOSOBL". IF C THIS IS POSITIVE AND LESS THAN 2**30, THEN FLAG(2) = .TRUE. C (WE ASSUME WE ARE WORKING ON A COMPUTER WITH C WORD LENGTH AT LEAST 31 BITS EXCLUDING SIGN.) C THE NUMBER OF COLUMNS OF THE ARRAY V WHICH C ARE INITIALIZED IS C MAXCOL = NUMBER OF BITS IN ATMOST. C IN "GONEDX" WE CHECK THAT THIS IS NOT EXCEEDED. C C NX STORED NIEDERREiITER-XING SEQUENCE GENERATORS C WHICH IS BINARY FRACTIONS WITH RESPECT TO ROWS. C THE GENERATORS CAN BE FOUND URL OF WEBSITE C http://www.dismat.oeaw.ac.at/pirs/niedxing.html C C C THE NUMBERS IN V ARE ACTUALLY BINARY FRACTIONS WITH C RESPECT TO COLUMNS. C C LSM ARE LOWER TRIANGULAR SCRAMBLING MATRICES. C USM ARE UPPER TRIANGULAR SCRAMBLING MATRICES. C SV ARE SCRAMBLING GENERATOR MATRICES AND THE NUMBERS C ARE BINARY FRACTIONS. C "RECIPD" HOLDS 1/(THE COMMON DENOMINATOR OF ALL C OF THEM). C C C "INNEDX" IMPLICITLY COMPUTES THE FIRST SHIFTED C VECTOR "LASTQ", AND RETURN IT TO THE CALLING C PROGRAM. SUBSEQUENT VECTORS COME FROM "GONEDX". C "LASTQ" HOLDS NUMERATORS OF THE LAST VECTOR GENERATED. C C C Non-Standard Intrinsic Funtion for f77 C But Standard Intrinsic Fuction for f90 IBITS IS USED. C C INPUTS : C FROM USER'S PROGRAM : DIMEN, ATMOST, MAXS, C IFLAG, QUASI C C OUTPUTS : C TO USER'S PROGRAM : FLAG, LASTQ C TO "GONEDX" VIA /NIEDX/ : C SV, S, MAXCOL, COUNT, LASTQ, RECIPD C C C C CHECK PARAMETERS C C .. Scalar Arguments .. INTEGER ATMOST,DIMEN,IFLAG,MAXS C .. C .. Array Arguments .. DOUBLE PRECISION QUASI(40) LOGICAL FLAG(2) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Local Scalars .. DOUBLE PRECISION LL INTEGER I,J,K,L,P,PP,TEMP1,TEMP2,TEMP3,TEMP4 CHARACTER*20 NAME C .. C .. Local Arrays .. INTEGER LSM(40,31),NX(16,30),SHIFT(40),TV(40,31,31),USHIFT(30), + USM(30,30),V(40,30) C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. External Subroutines .. EXTERNAL GNCRML,GNCRMU C .. C .. Intrinsic Functions .. INTRINSIC IBITS,MOD C .. C .. Common blocks .. COMMON /NIEDX/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /NIEDX/ C .. S = DIMEN FLAG(1) = (S.GE.1 .AND. S.LE.16) FLAG(2) = (ATMOST.GT.0 .AND. ATMOST.LT.2**30) IF (.NOT. (FLAG(1).AND.FLAG(2))) RETURN I = ATMOST MAXCOL = 0 10 MAXCOL = MAXCOL + 1 I = I/2 IF (I.GT.0) GO TO 10 C C GET GENERATOR MATRICES C IF (S.LE.4) THEN NAME = 'nxs4m30' ELSE IF (S.EQ.5) THEN NAME = 'nxs5m30' ELSE IF (S.EQ.6) THEN NAME = 'nxs6m30' ELSE IF (S.LE.7) THEN NAME = 'nxs7m30' ELSE IF (S.EQ.8) THEN NAME = 'nxs8m30' ELSE IF (S.EQ.9) THEN NAME = 'nxs9m30' ELSE IF (S.EQ.10) THEN NAME = 'nxs10m30' ELSE IF (S.LE.11) THEN NAME = 'nxs11m30' ELSE IF (S.EQ.12) THEN NAME = 'nxs12m30' ELSE IF (S.EQ.13) THEN NAME = 'nxs13m30' ELSE IF (S.LE.14) THEN NAME = 'nxs14m30' ELSE IF (S.EQ.15) THEN NAME = 'nxs15m30' ELSE IF (S.EQ.16) THEN NAME = 'nxs16m30' END IF OPEN (UNIT=1,FILE=NAME,STATUS='UNKNOWN') DO 20 I = 1,S READ (1,FMT=*) (NX(I,J),J=1,30) 20 CONTINUE CLOSE (UNIT=1,STATUS='KEEP') DO 40 I = 1,S DO 30 J = 1,MAXCOL V(I,J) = 0 30 CONTINUE 40 CONTINUE L = 1 DO 70 J = 30,1,-1 DO 60 I = 1,S DO 50 K = 1,MAXCOL V(I,K) = V(I,K) + IBITS(NX(I,J),30-K,1)*L 50 CONTINUE 60 CONTINUE L = 2*L 70 CONTINUE C C COMPUTING GENERATOR MATRICES OF USER CHOICE C IF (IFLAG.EQ.0) THEN DO 90 I = 1,S DO 80 J = 1,MAXCOL SV(I,J) = V(I,J) 80 CONTINUE SHIFT(I) = 0 90 CONTINUE LL = 2.0** (30) ELSE IF ((IFLAG.EQ.1) .OR. (IFLAG.EQ.3)) THEN CALL GNCRML(MAXS,LSM,SHIFT) DO 130 I = 1,S DO 120 J = 1,MAXCOL L = 1 TEMP2 = 0 DO 110 P = MAXS,1,-1 TEMP1 = 0 DO 100 K = 1,30 TEMP1 = TEMP1 + (IBITS(LSM(I,P),K-1,1)* + IBITS(V(I,J),K-1,1)) 100 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L L = 2*L 110 CONTINUE SV(I,J) = TEMP2 120 CONTINUE 130 CONTINUE LL = 2.0** (MAXS) END IF IF ((IFLAG.EQ.2) .OR. (IFLAG.EQ.3)) THEN CALL GNCRMU(USM,USHIFT) IF (IFLAG.EQ.2) THEN MAXS = 30 END IF DO 190 I = 1,S DO 150 J = 1,MAXCOL P = MAXS DO 140 K = 1,MAXS IF (IFLAG.EQ.2) THEN TV(I,P,J) = IBITS(V(I,J),K-1,1) ELSE TV(I,P,J) = IBITS(SV(I,J),K-1,1) END IF P = P - 1 140 CONTINUE 150 CONTINUE DO 180 PP = 1,MAXCOL TEMP2 = 0 TEMP4 = 0 L = 1 DO 170 J = MAXS,1,-1 TEMP1 = 0 TEMP3 = 0 DO 160 P = 1,MAXCOL TEMP1 = TEMP1 + TV(I,J,P)*USM(P,PP) IF (PP.EQ.1) THEN TEMP3 = TEMP3 + TV(I,J,P)*USHIFT(P) END IF 160 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L IF (PP.EQ.1) THEN TEMP3 = MOD(TEMP3,2) TEMP4 = TEMP4 + TEMP3*L END IF L = 2*L 170 CONTINUE SV(I,PP) = TEMP2 IF (PP.EQ.1) THEN IF (IFLAG.EQ.3) THEN SHIFT(I) = EXOR(TEMP4,SHIFT(I)) ELSE SHIFT(I) = TEMP4 END IF END IF 180 CONTINUE 190 CONTINUE LL = 2.0** (MAXS) END IF END IF C C RECIPD IS 1/(COMMON DENOMINATOR OF THE ELEMENTS IN V) C RECIPD = 1.0/LL C C SET UP FIRST VECTOR AND VALUES FOR "GONEDX" C AND RETURN THE FIRST SEQUENCE C COUNT = 0 DO 200 I = 1,S LASTQ(I) = SHIFT(I) QUASI(I) = LASTQ(I)*RECIPD 200 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Ssobol' then mkdir 'Ssobol' fi cd 'Ssobol' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' include makefile.inc all: Res1 SRCLIBS= $(LINPACK) $(EISPACK_DP) $(BLAS) $(PORT) src.o: $(INCSRC)/src.f $(F77) $(F77OPTS) -c $(INCSRC)/src.f extra.o: $(INCSRC)/extra.f $(F77) $(F77OPTS) -c $(INCSRC)/extra.f DRIVERS= driver RESULTS= Res1 SRCLIBS = $(TIMER) Objs1= driver.o src.o extra.o driver: $(Objs1) $(F77LINK) $(F77LINKOPTS) -o driver $(Objs1) $(SRCLIBS) Res1: driver ./driver >Res1 diffres: Res1 res1 echo "Differences in results from driver" $(DIFF) Res1 res1 clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << "SHAR_EOF" > 'data1' 683644482 135580830 114750887 803113694 903684270 72157854 59862132 505680715 779216195 193586105 329775408 1049470377 790869458 529177858 770596904 1072006936 140712293 266923202 844524433 114947652 963169856 269847095 96724448 724232372 729133779 872385703 803906910 638059470 526484481 843345462 431170622 1057600395 413165996 637333254 738214251 987594285 394318066 718140016 1064908418 636856468 856354958 404426257 6726334 916596837 284819911 905844483 331078851 514311810 280380871 612727606 184724875 578053121 197887820 50641495 423123265 227871988 36070792 465200006 187446371 99681590 983367218 414604440 243776932 990984772 289169669 72196638 522106132 18453143 1008259046 521211985 88531465 816692680 1001849242 432014653 124800720 554138134 166704786 634462798 1019392879 691725817 459507961 1043416360 536793729 1046579554 403993816 412804346 138344837 24248109 184604157 535890993 869278100 125797533 322976214 161488107 755812753 541238672 841393456 909350792 521628897 567331204 118313388 265208159 468752989 169896169 959104945 500814545 434091372 113024960 444725900 426474642 263586462 432120698 45230250 524531296 327681208 288240183 353741403 192041686 42986304 258035065 388012188 537416171 292523858 794840680 502322834 423539584 884204171 419181022 618098449 257684411 288733678 145066288 756812965 357038497 881037192 314634603 182012383 753118967 636654128 318546833 71382991 400554255 153755763 187747483 306343731 1044780 402486338 564947470 934172722 813806089 705671218 795868370 926341726 51211860 739601563 176287958 1016675820 950546833 267929114 275671035 666917527 431607592 177069271 570312219 125215260 886835363 75079442 767034988 879903592 525306349 560816857 973795128 774079162 143751759 137446578 209360240 1012892947 398130103 838233295 47600869 312065869 24489594 686091463 1064468817 610953508 593516621 736434891 331884657 295796744 991620199 197203105 694133440 834981040 813915471 1043396296 246522553 741863988 188767345 707349375 236730174 943797129 1073073438 917394477 81723820 598848363 17981413 322240492 766981617 282924362 643850936 914020733 586984616 128906358 739978294 543175409 351896331 1017052692 802102481 946877967 615436223 305014361 1038674293 372513724 902652345 552583238 507262995 36761397 753350193 116032862 328670174 103551339 1045088961 667010209 737933140 945361853 324875181 353779021 866388681 255506551 65493883 220170356 837239777 771797016 542598812 568193573 359763920 870736829 304824104 430870768 116534859 741311666 157845915 861348937 985097428 714048271 306167010 709889598 498574308 382576036 360956609 1005272948 680195561 845547054 1026217890 374002007 866217449 695335995 462692419 314161108 1048195907 709323042 604523286 581285813 463752482 919585493 518104348 905143615 884163267 99734344 332656669 398822849 973544107 459181032 1001255433 846689598 32962550 6372759 727573859 495342420 738892195 834835824 1047729137 529022179 1032177776 185764687 483299066 818974724 826550698 987239779 684232759 940422073 224992310 164068454 827234589 690302353 32197089 628566875 95166490 866947171 476745989 604118539 121738255 667653407 987681310 749263282 863528381 341372968 244006930 951984484 614173490 685305366 592527405 383027710 952965461 21271273 653806680 745437736 180332506 331728868 207307411 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' PROGRAM TSSOB C C C THIS PROGRAM TESTS ACCURACY OF C NUMERICAL INTEGRATION USING "GNSSOB" C AND INTEGRAND (2) OF DAVIS AND C RABINOWITZ, PAGE 406 C C User Define: C DIMEN : dimension C ATMOST : sequence length C SAMS : Number of replications C MAXS : Maximum Digits of Scrambling Of Owen type Scrambling C IFLAG: User Choice of type Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C .. Local Scalars .. DOUBLE PRECISION F,SUM INTEGER ATMOST,DIMEN,I,II,IFLAG,K,J,MAXS,SAM C .. C .. Array Arguments .. DOUBLE PRECISION OUTS(40,10000) REAL SECOND, T1 C .. C .. C .. External Subroutines .. EXTERNAL GNSSOB,TIME C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. SAM = 1 MAXS = 31 DIMEN = 2 IFLAG = 1 ATMOST = 2**12 DO 30 II = 1,SAM WRITE (*,FMT=*) 'I = ITERATION NUMBER' WRITE (*,FMT=*) 'EI = ESTIMATED INTEGRAL' T1 = SECOND() CALL GNSSOB(DIMEN,ATMOST,IFLAG,MAXS,OUTS) SUM = 0.0 K = 9 DO 20 I = 1,ATMOST F = 1.0 DO 10 J = 1,DIMEN F = F*ABS(4.0*OUTS(J,I)-2.0) 10 CONTINUE IF (MOD(I,2**K).EQ.0) THEN K = K + 1 WRITE (*,FMT=*) 'I = ',I WRITE (*,FMT=*) 'EI = ',SUM/I END IF SUM = SUM + F 20 CONTINUE T1 = SECOND() - T1 WRITE (*,FMT=*) 'Total time elapsed = ',T1 30 CONTINUE STOP END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' I = ITERATION NUMBER EI = ESTIMATED INTEGRAL I = 512 EI = 0.9957062293734682 I = 1024 EI = 0.998864951179457 I = 2048 EI = 0.9991789653087108 I = 4096 EI = 0.9999722196422461 Total time elapsed = 0.029989865 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'extra.f' then echo shar: will not over-write existing file "'extra.f'" else cat << "SHAR_EOF" > 'extra.f' INTEGER FUNCTION EXOR(IIN,JIN) C C THIS FUNCTION CALCULATES THE EXCLUSIVE-OR OF ITS C TWO INPUT PARAMETERS C C .. Scalar Arguments .. INTEGER IIN,JIN C .. C .. Local Scalars .. INTEGER I,J,K,L C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. I = IIN J = JIN K = 0 L = 1 C 10 IF (I.EQ.J) THEN EXOR = K RETURN END IF C C CHECK THE CURRENT RIGHT-HAND BITS OF I AND J. C IF THEY DIFFER, SET THE APPROPRIATE BIT OF K. C IF (MOD(I,2).NE.MOD(J,2)) K = K + L I = I/2 J = J/2 L = 2*L GO TO 10 END DOUBLE PRECISION FUNCTION UNI() * * Random number generator, adapted from F. James * "A Review of Random Number Generators" * Comp. Phys. Comm. 60(1990), pp. 329-344. * C .. Parameters .. DOUBLE PRECISION TWOM24 PARAMETER (TWOM24=1D0/16777216.0) C .. C .. Local Scalars .. DOUBLE PRECISION CARRY INTEGER I,J C .. C .. Local Arrays .. DOUBLE PRECISION SEEDS(24) C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Save statement .. SAVE I,J,CARRY,SEEDS C .. C .. Data statements .. DATA I,J,CARRY/24,10,0.0d0/ DATA SEEDS/0.8804418d0,0.2694365d0,0.0367681d0,0.4068699d0, +0.4554052d0, + 0.2880635d0,0.1463408d0,0.2390333d0,0.6407298d0, +0.1755283d0,0.7132940d0, + 0.4913043d0,0.2979918d0,0.1396858d0,0.3589528d0, +0.5254809d0,0.9857749d0, + 0.4612127d0,0.2196441d0,0.7848351d0,0.4096100d0, +0.9807353d0,0.2689915d0, + 0.5140357d0/ C .. UNI = SEEDS(I) - SEEDS(J) - CARRY IF (UNI.LT.0) THEN UNI = UNI + 1 CARRY = TWOM24 ELSE CARRY = 0 END IF SEEDS(I) = UNI I = 24 - MOD(25-I,24) J = 24 - MOD(25-J,24) RETURN END SHAR_EOF fi # end of overwriting check if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' BLOCK DATA BDSOBL C C INITIALIZES LABELLED COMMON /SOBDAT/ C FOR "INSOBL". C C THE ARRAY POLY GIVES SUCCESSIVE PRIMITIVE C POLYNOMIALS CODED IN BINARY, E.G. C 45 = 100101 C HAS BITS 5, 2, AND 0 SET (COUNTING FROM THE C RIGHT) AND THEREFORE REPRESENTS C X**5 + X**2 + X**0 C C THESE POLYNOMIALS ARE IN THE ORDER USED BY C SOBOL IN USSR COMPUT. MATHS. MATH. PHYS. 16 (1977), C 236-242. A MORE COMPLETE TABLE IS GIVEN IN SOBOL AND C LEVITAN, THE PRODUCTION OF POINTS UNIFORMLY C DISTRIBUTED IN A MULTIDIMENSIONAL CUBE (IN RUSSIAN), C PREPRINT IPM AKAD. NAUK SSSR, NO. 40, MOSCOW 1976. C C THE INITIALIZATION OF THE ARRAY VINIT IS FROM THE C LATTER PAPER. FOR A POLYNOMIAL OF DEGREE M, M INITIAL C VALUES ARE NEEDED : THESE ARE THE VALUES GIVEN HERE. C SUBSEQUENT VALUES ARE CALCULATED IN "INSOBL". C C Non-Standard Intrinsic Funtion for f77 C But Standard Intrinsic Fuction for f90 IBITS IS USED. C C .. Arrays in Common .. INTEGER POLY(2:40),VINIT(2:40,8) C .. C .. Local Scalars .. INTEGER I C .. C .. Common blocks .. COMMON /SOBDAT/POLY,VINIT C .. C .. Save statement .. SAVE /SOBDAT/ C .. C .. Data statements .. C C DATA POLY/3,7,11,13,19,25,37,59,47,61,55,41,67,97,91,109,103,115, + 131,193,137,145,143,241,157,185,167,229,171,213,191,253,203, + 211,239,247,285,369,299/ DATA (VINIT(I,1),I=2,40)/39*1/ DATA (VINIT(I,2),I=3,40)/1,3,1,3,1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3, + 1,3,3,1,3,1,3,1,3,1,1,3,1,3,1,3,1,3/ DATA (VINIT(I,3),I=4,40)/7,5,1,3,3,7,5,5,7,7,1,3,3,7,5,1,1,5,3,3, + 1,7,5,1,3,3,7,5,1,1,5,7,7,5,1,3,3/ DATA (VINIT(I,4),I=6,40)/1,7,9,13,11,1,3,7,9,5,13,13,11,3,15,5,3, + 15,7,9,13,9,1,11,7,5,15,1,15,11,5,3,1,7,9/ DATA (VINIT(I,5),I=8,40)/9,3,27,15,29,21,23,19,11,25,7,13,17,1,25, + 29,3,31,11,5,23,27,19,21,5,1,17,13,7,15,9,31,9/ DATA (VINIT(I,6),I=14,40)/37,33,7,5,11,39,63,27,17,15,23,29,3,21, + 13,31,25,9,49,33,19,29,11,19,27,15,25/ DATA (VINIT(I,7),I=20,40)/13,33,115,41,79,17,29,119,75,73,105,7, + 59,65,21,3,113,61,89,45,107/ DATA (VINIT(I,8),I=38,40)/7,23,39/ C .. C END SUBROUTINE GNCRML(MAXS,LSM,SHIFT) C GENERATING LOWER TRIANULAR SCRMABLING MATRICES AND SHIFT VECTORS. C INPUTS : C FROM INSOBL : MAXS C FROM BLOCK DATA "SOBOL" : S, MAXCOL, C C OUTPUTS : C TO INSOBL : LSM, SHIFT C .. Scalar Arguments .. INTEGER MAXS C .. C .. Array Arguments .. INTEGER LSM(40,31),SHIFT(40) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Local Scalars .. INTEGER I,J,L,LL,P,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Common blocks .. COMMON /SOBOL/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /SOBOL/ C .. DO 30 P = 1,S SHIFT(P) = 0 L = 1 DO 20 I = MAXS,1,-1 LSM(P,I) = 0 STEMP = MOD((INT(UNI()*1000.0)),2) SHIFT(P) = SHIFT(P) + STEMP*L L = 2*L LL = 1 DO 10 J = MAXCOL,1,-1 IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.LT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF LSM(P,I) = LSM(P,I) + TEMP*LL LL = 2*LL 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GNCRMU(USM,USHIFT) C GENERATING UPPER TRIANGULAR SCRMABLING MATRICES AND C SHIFT VECTORS. C INPUTS : C FROM BLOCK DATA "SOBOL" : S, MAXCOL, C C OUTPUTS : C TO INSOBL : USM, USHIFT C .. Array Arguments .. INTEGER USHIFT(31),USM(31,31) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Local Scalars .. INTEGER I,J,STEMP,TEMP C .. C .. External Functions .. DOUBLE PRECISION UNI EXTERNAL UNI C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD C .. C .. Common blocks .. COMMON /SOBOL/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /SOBOL/ C .. DO 20 I = 1,MAXCOL STEMP = MOD((INT(UNI()*1000.0)),2) USHIFT(I) = STEMP DO 10 J = 1,MAXCOL IF (J.EQ.I) THEN TEMP = 1 ELSE IF (J.GT.I) THEN TEMP = MOD((INT(UNI()*1000.0)),2) ELSE TEMP = 0 END IF USM(I,J) = TEMP 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE GNSSOB(DIMEN,ATMOST,IFLAG,MAXS,OUTS) C C Subroutine for generating Sobol' sequence. C C User Define: C DIMEN : dimension C ATMOST : sequence length C MAXS : Maximum Digits of Scrambling Of Owen type Scrambling C IFLAG: User Choice of Sequences C IFLAG = 0 : No Scrambling C IFLAG = 1 : Owen type Scrambling C IFLAG = 2 : Faure-Tezuka type Scrambling C IFLAG = 3 : Owen + Faure-Tezuka type Scrambling C C .. Local Scalars .. INTEGER ATMOST,DIMEN,I,IFLAG,J,MAXS,TAUS C .. C .. Local Arrays .. DOUBLE PRECISION QUASI(40),OUTS(40,10000) LOGICAL FLAG(2) C .. C .. External Subroutines .. EXTERNAL GOSOBL,INSOBL CALL INSOBL(FLAG,DIMEN,ATMOST,TAUS,QUASI,MAXS,IFLAG) DO 10 J = 1,DIMEN OUTS(J,1) = QUASI(J) 10 CONTINUE DO 30 I = 2,ATMOST CALL GOSOBL(QUASI) DO 20 J = 1,DIMEN OUTS(J,I) = QUASI(J) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE GOSOBL(QUASI) C C THIS SUBROUTINE GENERATES A NEW C QUASIRANDOM VECTOR WITH EACH CALL C C IT ADAPTS THE IDEAS OF ANTONOV AND SALEEV, C USSR COMPUT. MATHS. MATH. PHYS. 19 (1980), C 252 - 256 C C THE USER MUST CALL "INSOBL" BEFORE CALLING C "GOSOBL". AFTER CALLING "INSOBL", TEST C FLAG(1) AND FLAG(2); IF EITHER IS FALSE, C DO NOT CALL "GOSOBL". "GOSOBL" CHECKS C THAT THE USER DOES NOT MAKE MORE CALLS C THAN HE SAID HE WOULD : SEE THE COMMENTS C TO "INSOBL". C C INPUTS: C FROM USER'S CALLING PROGRAM: C NONE C C FROM LABELLED COMMON /SOBOL/: C SV TABLE OF DIRECTION NUMBERS C S DIMENSION C MAXCOL LAST COLUMN OF V TO BE USED C COUNT SEQUENCE NUMBER OF THIS CALL C LASTQ NUMERATORS FOR LAST VECTOR GENERATED C RECIPD (1/DENOMINATOR) FOR THESE NUMERATORS C C C FIND THE POSITION OF THE RIGHT-HAND ZERO IN COUNT C C .. Array Arguments .. DOUBLE PRECISION QUASI(40) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),SV(40,31) C .. C .. Local Scalars .. INTEGER I,L C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. Intrinsic Functions .. INTRINSIC MOD C .. C .. Common blocks .. COMMON /SOBOL/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /SOBOL/ C .. L = 0 I = COUNT 10 L = L + 1 IF (MOD(I,2).EQ.1) THEN I = I/2 GO TO 10 END IF C C CHECK THAT THE USER IS NOT CHEATING ! C IF (L.GT.MAXCOL) STOP ' TOO MANY CALLS ON GOSOBL' C C CALCULATE THE NEW COMPONENTS OF QUASI, C FIRST THE NUMERATORS, THEN NORMALIZED C DO 20 I = 1,S LASTQ(I) = EXOR(LASTQ(I),SV(I,L)) C C IF A FULL-WORD EXCLUSIVE-OR, SAY .XOR., IS AVAILABLE C THEN REPLACE THE PRECEDING STATEMENT BY C C LASTQ(I) = LASTQ(I) .XOR. SV(I,L) C C TO GET A FASTER, EXTENDED FORTRAN PROGRAM C QUASI(I) = LASTQ(I)*RECIPD 20 CONTINUE C COUNT = COUNT + 1 C RETURN END SUBROUTINE INSOBL(FLAG,DIMEN,ATMOST,TAUS,QUASI,MAXS,IFLAG) C C THIS IS MODIFIED ROUTINE OF "INSOBL". C FIRST CHECK WHETHER THE USER-SUPPLIED C DIMENSION "DIMEN" OF THE QUASI-RANDOM C VECTORS IS STRICTLY BETWEEN 1 AND 41. C IF SO, FLAG(1) = .TRUE. C C NEXT CHECK "ATMOST", AN UPPER BOUND ON THE NUMBER C OF CALLS THE USER INTENDS TO MAKE ON "GOSOBL". IF C THIS IS POSITIVE AND LESS THAN 2**30, THEN FLAG(2) = .TRUE. C (WE ASSUME WE ARE WORKING ON A COMPUTER WITH C WORD LENGTH AT LEAST 31 BITS EXCLUDING SIGN.) C THE NUMBER OF COLUMNS OF THE ARRAY V WHICH C ARE INITIALIZED IS C MAXCOL = NUMBER OF BITS IN ATMOST. C IN "GOSOBL" WE CHECK THAT THIS IS NOT EXCEEDED. C C THE LEADING ELEMENTS OF EACH ROW OF V ARE C INITIALIZED USING "VINIT" FROM "BDSOBL". C EACH ROW CORRESPONDS TO A PRIMITIVE POLYNOMIAL C (AGAIN, SEE "BDSOBL"). IF THE POLYNOMIAL HAS C DEGREE M, ELEMENTS AFTER THE FIRST M ARE CALCULATED. C C THE NUMBERS IN V ARE ACTUALLY BINARY FRACTIONS. C LSM ARE LOWER TRIAUGULAR SCRAMBLING MATRICES. C USM ARE UPPER TRIAUGULAR SCRMABLING MATRIX. C SV ARE SCAMBLING GENERATING MATRICES AND THE NUMBERS C ARE BINARY FRACTIONS. C "RECIPD" HOLDS 1/(THE COMMON DENOMINATOR OF ALL C OF THEM). C C C "INSOBL" IMPLICITLY COMPUTES THE FIRST SHIFTED C VECTOR "LASTQ", AND RETURN IT TO THE CALLING C PROGRAM. SUBSEQUENT VECTORS COME FROM "GOSOBL". C "LASTQ" HOLDS NUMERATORS OF THE LAST VECTOR GENERATED. C C "TAUS" IS FOR DETERMINING "FAVORABLE" VALUES. AS C DISCUSSED IN BRATLEY/FOX, THESE HAVE THE FORM C N = 2**K WHERE K .GE. (TAUS+S-1) FOR INTEGRATION C AND K .GT. TAUS FOR GLOBAL OPTIMIZATION. C C INPUTS : C FROM USER'S PROGRAM : DIMEN, ATMOST, MAXS, IFLAG C FROM BLOCK DATA "BDSOBL" : POLY, VINIT C C OUTPUTS : C TO USER'S PROGRAM : FLAG, TAUS, LASTQ,QUASI C TO "GOSOBL" VIA /SOBOL/ : C SV, S, MAXCOL, COUNT, LASTQ, RECIPD C C C .. Scalar Arguments .. INTEGER ATMOST,DIMEN,IFLAG,MAXS,TAUS C .. C .. Array Arguments .. DOUBLE PRECISION QUASI(40) LOGICAL FLAG(2) C .. C .. Scalars in Common .. DOUBLE PRECISION RECIPD INTEGER COUNT,MAXCOL,S C .. C .. Arrays in Common .. INTEGER LASTQ(40),POLY(2:40),SV(40,31),VINIT(2:40,8) C .. C .. Local Scalars .. DOUBLE PRECISION LL INTEGER I,J,K,L,M,MAXX,NEWV,P,PP,TEMP1,TEMP2,TEMP3,TEMP4 C .. C .. Local Arrays .. INTEGER LSM(40,31),SHIFT(40),TAU(13),TV(40,31,31),USHIFT(31), + USM(31,31),V(40,31) LOGICAL INCLUD(8) C .. C .. External Functions .. INTEGER EXOR EXTERNAL EXOR C .. C .. External Subroutines .. EXTERNAL GNCRML,GNCRMU C .. C .. Intrinsic Functions .. INTRINSIC IBITS,MOD C .. C .. Common blocks .. COMMON /SOBDAT/POLY,VINIT COMMON /SOBOL/S,MAXCOL,SV,COUNT,LASTQ,RECIPD C .. C .. Save statement .. SAVE /SOBDAT/,/SOBOL/ C .. C .. Data statements .. DATA TAU/0,0,1,3,5,8,11,15,19,23,27,31,35/ C .. C C CHECK PARAMETERS C S = DIMEN FLAG(1) = (S.GE.1 .AND. S.LE.40) FLAG(2) = (ATMOST.GT.0 .AND. ATMOST.LT.2**30) IF (.NOT. (FLAG(1).AND.FLAG(2))) RETURN IF (S.LE.13) THEN TAUS = TAU(S) ELSE TAUS = -1 C RETURN A DUMMY VALUE TO THE CALLING PROGRAM END IF * C C FIND NUMBER OF BITS IN ATMOST C I = ATMOST MAXCOL = 0 10 MAXCOL = MAXCOL + 1 I = I/2 IF (I.GT.0) GO TO 10 C C INITIALIZE ROW 1 OF V C DO 20 I = 1,MAXCOL V(1,I) = 1 20 CONTINUE C C INITIALIZE REMAINING ROWS OF V C DO 80 I = 2,S C C THE BIT PATTERN OF POLYNOMIAL I GIVES ITS FORM C (SEE COMMENTS TO "BDSOBL") C FIND DEGREE OF POLYNOMIAL I FROM BINARY ENCODING C J = POLY(I) M = 0 30 J = J/2 IF (J.GT.0) THEN M = M + 1 GO TO 30 END IF C C WE EXPAND THIS BIT PATTERN TO SEPARATE COMPONENTS C OF THE LOGICAL ARRAY INCLUD. C J = POLY(I) DO 40 K = M,1,-1 INCLUD(K) = (MOD(J,2).EQ.1) J = J/2 40 CONTINUE C C THE LEADING ELEMENTS OF ROW I COME FROM VINIT C DO 50 J = 1,M V(I,J) = VINIT(I,J) 50 CONTINUE C C CALCULATE REMAINING ELEMENTS OF ROW I AS EXPLAINED C IN BRATLEY AND FOX, SECTION 2 C DO 70 J = M + 1,MAXCOL NEWV = V(I,J-M) L = 1 DO 60 K = 1,M L = 2*L IF (INCLUD(K)) NEWV = EXOR(NEWV,L*V(I,J-K)) C C IF A FULL-WORD EXCLUSIVE-OR, SAY .XOR., IS AVAILABLE, C THEN REPLACE THE PRECEDING STATEMENT BY C C IF (INCLUD(K)) NEWV = NEWV .XOR. (L * V(I, J-K)) C C TO GET A FASTER, EXTENDED FORTRAN PROGRAM C 60 CONTINUE V(I,J) = NEWV 70 CONTINUE C 80 CONTINUE C C MULTIPLY COLUMNS OF V BY APPROPRIATE POWER OF 2 C L = 1 DO 100 J = MAXCOL - 1,1,-1 L = 2*L DO 90 I = 1,S V(I,J) = V(I,J)*L 90 CONTINUE 100 CONTINUE C C COMPUTING GENERATOR MATRICES OF USER CHOICE C IF (IFLAG.EQ.0) THEN DO 120 I = 1,S DO 110 J = 1,MAXCOL SV(I,J) = V(I,J) 110 CONTINUE SHIFT(I) = 0 120 CONTINUE LL = 2.0** (MAXCOL) ELSE IF ((IFLAG.EQ.1) .OR. (IFLAG.EQ.3)) THEN CALL GNCRML(MAXS,LSM,SHIFT) DO 160 I = 1,S DO 150 J = 1,MAXCOL L = 1 TEMP2 = 0 DO 140 P = MAXS,1,-1 TEMP1 = 0 DO 130 K = 1,MAXCOL TEMP1 = TEMP1 + (IBITS(LSM(I,P),K-1,1)* + IBITS(V(I,J),K-1,1)) 130 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L L = 2*L 140 CONTINUE SV(I,J) = TEMP2 150 CONTINUE 160 CONTINUE LL = 2.0** (MAXS) END IF IF ((IFLAG.EQ.2) .OR. (IFLAG.EQ.3)) THEN CALL GNCRMU(USM,USHIFT) IF (IFLAG.EQ.2) THEN MAXX = MAXCOL ELSE MAXX = MAXS END IF DO 220 I = 1,S DO 180 J = 1,MAXCOL P = MAXX DO 170 K = 1,MAXX IF (IFLAG.EQ.2) THEN TV(I,P,J) = IBITS(V(I,J),K-1,1) ELSE TV(I,P,J) = IBITS(SV(I,J),K-1,1) END IF P = P - 1 170 CONTINUE 180 CONTINUE DO 210 PP = 1,MAXCOL TEMP2 = 0 TEMP4 = 0 L = 1 DO 200 J = MAXX,1,-1 TEMP1 = 0 TEMP3 = 0 DO 190 P = 1,MAXCOL TEMP1 = TEMP1 + TV(I,J,P)*USM(P,PP) IF (PP.EQ.1) THEN TEMP3 = TEMP3 + TV(I,J,P)*USHIFT(P) END IF 190 CONTINUE TEMP1 = MOD(TEMP1,2) TEMP2 = TEMP2 + TEMP1*L IF (PP.EQ.1) THEN TEMP3 = MOD(TEMP3,2) TEMP4 = TEMP4 + TEMP3*L END IF L = 2*L 200 CONTINUE SV(I,PP) = TEMP2 IF (PP.EQ.1) THEN IF (IFLAG.EQ.3) THEN SHIFT(I) = EXOR(TEMP4,SHIFT(I)) ELSE SHIFT(I) = TEMP4 END IF END IF 210 CONTINUE 220 CONTINUE LL = 2.0** (MAXX) END IF END IF C C RECIPD IS 1/(COMMON DENOMINATOR OF THE ELEMENTS IN V) C RECIPD = 1.0/LL C C SET UP FIRST VECTOR AND VALUES FOR "GOSOBL" C COUNT = 0 DO 230 I = 1,S LASTQ(I) = SHIFT(I) QUASI(I) = LASTQ(I)*RECIPD 230 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. cd .. # End of shell archive exit 0