C-------------- SAMPLE PROGRAM TO ENCIPHER 72-BIT PASSWORDS ACROSS 00000010 C MACHINES WITH 3 DIFFERENT WORD SIZES. OUTPUT IS TABLE 2A OR 2B. 00000020 C TO RUN THIS PROGRAM COMMENTS BEGINNING C32, C36, OR C60 MUST BE 00000030 C ACTIVATED BY REPLACING THESE FIRST THREE CHARACTERS WITH BLANKS. 00000040 C TABLE 2A WILL BE PRODUCED ON IBM 360/370 OR XEROX SIGMA SERIES 00000050 C 32-BIT MACHINES BY ACTIVATING COMMENTS BEGINNING WITH C32. 00000060 C TABLE 2B WILL BE PRODUCED ON UNIVAC 1100, HONEYWELL 600/6000, AND 00000070 C DEC SYSTEM 10/20 36-BIT MACHINES WHEN ACTIVATING COMMENTS BEGINNING 00000080 C WITH C36. TABLE 2B WILL BE PRODUCED ON THE CDC 6600/7600 60-BIT 00000090 C MACHINES WHEN ACTIVATING COMMENTS BEGINNING WITH C60. 00000100 C 00000110 C INPUT TO THIS PROGRAM ARE THE FOLLOWING 10 CARDS BEGINNING WITH C= 00000120 C= HEX(000000000000000000) OCTAL(000000000000000000000000) DECIMAL(0) 00000130 C= HEX(000000000000000001) OCTAL(000000000000000000000001) DECIMAL(1) 00000140 C= HEX(000000000000000002) OCTAL(000000000000000000000002) DECIMAL(2) 00000150 C= HEX(000000000000000003) OCTAL(000000000000000000000003) DECIMAL(3) 00000160 C= HEX(FFFFFFFFFFFFFFFFA3) OCTAL(777777777777777777777643) DECIMAL(P) 00000170 C= HEX(FFFFFFFFFFFFFFFFA4) OCTAL(777777777777777777777644) DECIMAL(P+1) 00000180 C= HEX(FFFFFFFFFFFFFFFFA5) OCTAL(777777777777777777777645) DECIMAL(P+2) 00000190 C= HEX(00000000000000005C) OCTAL(000000000000000000000134) DECIMAL(92) 00000200 C= HEX(FFFFFFFFFFFFFFFFFF) OCTAL(777777777777777777777777) DECIMAL(P+92)00000210 C= HEX(555555555555555555) OCTAL(252525252525252525252525) BINARY(01...)00000220 C 00000230 C-----------------------------------------------------------------------00000240 INTEGER A,NN,N(6) 00000250 C32 INTEGER X(3),FX(3),P(3),AC(3,6),WORK(43),NP,IW 00000260 C36 INTEGER X(2),FX(2),P(2),AC(2,6),WORK(31),NP,IW 00000270 C60 INTEGER X(2),FX(2),P(2),AC(2,6),WORK(31),NP,IW 00000280 C-----------A IS CHOSEN TO MAKE P=2**72-A A PRIME INTEGER IN THE 00000290 C INTERVAL (2**(NP*M)-2**(M-1), 2**(NP*M)). NP*M=72 HERE. 00000300 DATA A/93/, NN/6/ 00000310 DATA N(1)/1048569/,N(2)/524285/,N(3)/3/,N(4)/2/,N(5)/1/,N(6)/0/ 00000320 C.32.....IBM, XEROX 32-BIT MACHINES, HEXADECIMAL CONSTANTS(Z). 00000330 C32 DATA NP/3/, IW/43/ 00000340 C32 DATA AC(1,1),AC(2,1),AC(3,1)/Z00000000,Z00000000,Z00000001/, 00000350 C32 * AC(1,2),AC(2,2),AC(3,2)/Z00FFFFFF,Z00FFFFFF,Z00FFFFA1/, 00000360 C32 * AC(1,3),AC(2,3),AC(3,3)/Z00FFFFFF,Z00FFFFFF,Z00FFFF83/, 00000370 C32 * AC(1,4),AC(2,4),AC(3,4)/Z00FFFFFF,Z00FFFFFF,Z00FFFF7B/, 00000380 C32 * AC(1,5),AC(2,5),AC(3,5)/Z00FFFFFF,Z00FFFFFF,Z00FFFF9E/, 00000390 C32 * AC(1,6),AC(2,6),AC(3,6)/Z00FFFFFF,Z00FFFFFF,Z00FFFF02/ 00000400 C 00000410 C.36......UNIVAC, HONEYWELL, DEC 36-BIT MACHINES, OCTAL CONSTANTS(O). 00000420 C36 DATA NP/2/, IW/31/ 00000430 C36 DATA AC(1,1),AC(2,1)/O000000000000,O000000000001/, 00000440 C36 * AC(1,2),AC(2,2)/O777777777777,O777777777641/, 00000450 C36 * AC(1,3),AC(2,3)/O777777777777,O777777777603/, 00000460 C36 * AC(1,4),AC(2,4)/O777777777777,O777777777573/, 00000470 C36 * AC(1,5),AC(2,5)/O777777777777,O777777777636/, 00000480 C36 * AC(1,6),AC(2,6)/O777777777777,O777777777402/ 00000490 C 00000500 C.60......CDC 6600/7600 60-BIT +-, 48-BIT */, OCTAL CONSTANTS(B). 00000510 C60 DATA NP/2/, IW/31/ 00000520 C60 DATA AC(1,1),AC(2,1)/00000000000000000000B,00000000000000000001B/,00000530 C60 * AC(1,2),AC(2,2)/00000000777777777777B,00000000777777777641B/,00000540 C60 * AC(1,3),AC(2,3)/00000000777777777777B,00000000777777777603B/,00000550 C60 * AC(1,4),AC(2,4)/00000000777777777777B,00000000777777777573B/,00000560 C60 * AC(1,5),AC(2,5)/00000000777777777777B,00000000777777777636B/,00000570 C60 * AC(1,6),AC(2,6)/00000000777777777777B,00000000777777777402B/ 00000580 C 00000590 WRITE(6,1) 00000600 C32 1 FORMAT(45H1HEXADECIMAL INPUT(X) HEXADECIMAL OUTPUT(FX)/1X) 00000610 C36 1 FORMAT(15H1OCTAL INPUT(X),14X,16HOCTAL OUTPUT(FX)/1X) 00000620 C60 1 FORMAT(15H1OCTAL INPUT(X),14X,16HOCTAL OUTPUT(FX)/1X) 00000630 C---------INPUT AND ENCIPHER 10 BIT STRINGS. 00000640 DO 4 I=1,10 00000650 READ(5,2) X 00000660 C32 2 FORMAT(7X,3Z6) 00000670 C36 2 FORMAT(33X,2O12) 00000680 C60 2 FORMAT(33X,2O12) 00000690 CALL PURDY(X,FX,A,P,N,AC,WORK,NP,NN,IW,NP,IERR) 00000700 IF(IERR.NE.0) GO TO 5 00000710 WRITE(6,3) X,FX 00000720 C32 3 FORMAT(1X,3Z6,4X,3Z6) 00000730 C36 3 FORMAT(1X,2O12,4X,2O12) 00000740 C60 3 FORMAT(1X,2O12,4X,2O12) 00000750 4 CONTINUE 00000760 STOP 00000770 5 WRITE(6,6) IERR 00000780 6 FORMAT(46H0****ERROR RETURN FROM SUBROUTINE PURDY. IERR=,I2) 00000790 STOP 00000800 END 00000810 BLOCK DATA 00000820 COMMON/MACHIN/ ISIZE,M,MAXPOS,MINNEG,MASK 00000830 C32 DATA ISIZE/32/,M/24/,MAXPOS/Z007FFFFF/,MINNEG/Z00800000/, 00000840 C32 * MASK/Z01000000/ 00000850 C36 DATA ISIZE/36/,M/36/,MAXPOS/O377777777777/,MINNEG/O400000000000/, 00000860 C36 * MASK/0/ 00000870 C60 DATA ISIZE/60/,M/36/,MAXPOS/00000000377777777777B/, 00000880 C60 * MINNEG/00000000400000000000B/,MASK/00000001000000000000B/ 00000890 END 00000900 SUBROUTINE PURDY(X,FX,A,PRIME,N,ACOEFF,WORK,NP,NN,NW,IDIMA,IERR) 00000010 C=======ROUTINE TO EVALUATE PURDY'S IRREVERSIBLE ENCIPHERING FUNCTION... C FX=F(X)=P(X) MOD PRIME WHERE C P(X)= SUM(ACOEFF(I)*X**N(I)), I=1,2,3,...,NN. C THAT IS, THIS ROUTINE EVALUATES SPARSE POLYNOMIALS OF VERY LARGE C DEGREE MODULO A PRIME. SUCH FAMILIES OF ENCIPHERING FUNCTIONS AND C RELATED PROPERTIES ARE DISCUSSED IN... C A HIGH SECURITY LOG-IN PROCEDURE, BY GEORGE B. PURDY, AUGUST 1974, C CACM 17(8), PAGES 442-444. C C C PARAMETERS... C THE LABELED COMMON/MACHIN/ DESCRIBES PARAMETERS RELATING TO THE C HOST MACHINE'S WORD SIZE AND MUST BE INITIALIZED TO THE FOLLOWING C VALUES BEFORE SUBROUTINE PURDY IS CALLED... C C ISIZE IS GIVEN AS THE TOTAL NUMBER OF BITS IN AN INTEGER COMPUTER C WORD. C C M IS GIVEN AS AN EVEN INTEGER WHOSE VALUE IS THE NUMBER OF C SIGNIFICANT (LOW-ORDER) BITS IN EACH COMPUTER WORD TO BE C USED TO REPRESENT A MULTIPRECISION INTEGER. +,-,/,* MUST BE C VALID OVER THE NORMAL RANGES FOR INTEGER ARITHMETIC. C THAT IS, INTEGER ARITHMETIC MUST BE VALID OVER C -2**(M-1)-1 .LE. R .LE. 2**(M-1)-1 FOR 1'S COMPLEMENT, C AND SIGNED MAGNITUDE MACHINES AND VALID OVER THE RANGE C -2**(M-1) .LE. R .LE. 2**(M-1)-1 FOR 2'S COMPLEMENT MACHINES C WHERE R IS ANY COMPUTED RESULT. FOR SIGNED MAGNITUDE C MACHINES M MUST BE EVEN AND NO GREATER THAN ISIZE-2. C OTHERWISE M MUST BE EVEN AND LESS THAN OR EQUAL TO ISIZE. C C MAXPOS IS GIVEN AS THE BIT STRING VALUE OF 2**(M-1)-1, THAT IS THE C BIT STRING 0111...111 . C C MINNEG IS GIVEN AS THE BIT STRING VALUE OF 2**(M-1), THAT IS THE C BIT STRING 1000...000. C C MASK IS GIVEN AS THE BIT STRING VALUE OF 2**M MOD 2**ISIZE. THE C HOST COMPUTER MUST SUPPORT INTEGER ADDITION AND SUBTRACTION C WITH THE VALUE OF MASK AND ALL INTEGERS IN THE RANGE C DESCRIBED FOR M ABOVE. C C DEFINITION: A NP WORD MULTIPRECISION INTEGER IS AN UNSIGNED (NP*M)- C BIT NUMBER STORED IN NP SUCCESSIVE INTEGER COMPUTER C WORDS, THE HIGHEST ORDER PART BEGINNING WITH THE FIRST C WORD, THE LOWEST ORDER ENDING IN THE LAST. THAT IS, C EACH ELEMENT OF THE ARRAY REPRESENTING A MULTIPRECISION C INTEGER IS TREATED AS A BASE 2**M DIGIT WITH ISIZE-M C LEADING ZERO BITS. NORMALLY ISIZE=M OR M IS CHOSEN SUCH C THAT ISIZE-M IS AN EXACT MULTIPLE OF THE NUMBER OF BITS C USED TO REPRESENT A CHARACTER. FOLLOWING PURDY'S C NOTATION, IN FURTHER DISCUSSIONS HERE Q, THE PASSWORD C LENGTH, WILL BE ASSUMED TO BE Q=NP*M. C C C X IS GIVEN AS AN UNSIGNED NP-WORD MULTIPRECISION INTEGER C REPRESENTING A STRING OF Q BITS TO BE ENCIPHERED. C X NORMALLY REPRESENTS A CHARACTER STRING (PASSWORD). C C FX IS RETURNED AS AN UNSIGNED NP-WORD MULTIPRECISION INTEGER C WHOSE Q-BIT VALUE IS THE VALUE OF THE ENCIPHERING C FUNCTION, F(X), ABOVE. C C A IS GIVEN AS A POSITIVE INTEGER .LE. MAXPOS WHOSE VALUE IS THE C PRIME OFFSET CONSTANT DESCRIBED IN THE ABOVE REFERENCE. THAT C IS, A MUST BE CHOSEN TO FORCE THE VARIABLE PRIME TO BE PRIME C WHERE PRIME IS DEFINED AS PRIME=2**(M*NP) - A. SOME VALUES C ARE TABULATED IN KNUTH, ART OF COMPUTER PROGRAMMING, VOL. 2, C TABLE 4.5.4-1. NO CHECK IS MADE BY THE PROGRAM TO INSURE C THIS PRIMALITY CONDITION AS IT ONLY AFFECTS THE KNOWN C DEGENERACY OF THE ENCIPHERING FUNCTION, F. C C PRIME IS RETURNED AS AN UNSIGNED NP-WORD MULTIPRECISION INTEGER C WHOSE Q-BIT VALUE IS THE VALUE OF 2**(NP*M)-A. C C N IS GIVEN AS AN INTEGER ARRAY OF LENGTH NN WHOSE ELEMENTS ARE C THE EXPONENTS OF X. THE ELEMENTS OF N MUST BE NON-NEGATIVE. C TYPICALLY, AS SUGGESTED BY PURDY, THE FIRST ONE OR TWO C ELEMENTS OF N ARE VERY LARGE NON-PRIME INTEGERS, SAY C GREATER THAN 10**5. C C ACOEFF IS GIVEN AS AN IDIMA BY NN ARRAY REPRESENTING THE UNSIGNED C MULTIPRECISION INTEGER COEFFICIENTS OF THE ENCIPHERING C POLYNOMIAL F(X) ABOVE. THESE COEFFICIENTS SHOULD BE C CHOSEN TO HAVE (LARGE) VALUES, LESS THAN 2**Q-A. THE C PRECISION OF THE ELEMENTS OF ACOEFF IS NP AND THEREFORE C IDIMA USUALLY WOULD BE CHOSEN SUCH THAT IDIMA=NP. C C WORK IS GIVEN AS AN WORKSPACE ARRAY OF NW=12*NP+7 INTEGER WORDS. C C NP IS GIVEN AS AN INTEGER WHOSE VALUE IS THE PRECISION OF C MULTIPRECISION INTEGERS. THAT IS, NP IS THE NUMBER OF C INTEGER COMPUTER WORDS NECESSARY TO STORE THE UNSIGNED C MULTIPRECISION INTEGER 2**Q-1. NP MUST BE GREATER THAN 1. C C NN IS GIVEN AS AN INTEGER WHOSE VALUE IS THE NUMBER OF TERMS IN C P(X) ABOVE. NN MUST BE GREATER THAN 0. C C NW IS GIVEN AS AN INTEGER WHOSE VALUE IS GREATER THAN OR EQUAL C TO 12*NP+7. NW IS THE DIMENSION OF THE WORKSPACE ARRAY, WORK. C C IDIMA IS GIVEN AS AN INTEGER WHOSE VALUE IS EQUAL TO THE C FIRST DIMENSION CONSTANT OF THE ARGUMENT IN THE CALLING C PROGRAM CORRESPONDING TO THE PARAMETER, ACOEFF. NORMALLY C IDIMA=NP. C C IERR IS RETURNED AS AN ERROR INDICATOR AS FOLLOWS... C IERR=1 IF NP .LE. 1 OR NN IS NOT POSITIVE, IF NW IS LESS THAN 12*NP+7 C OR IF ANY OF THE NN ELEMENTS OF N ARE NEGATIVE. C IERR=2 IF A IS NOT IN ITS RANGE DESCRIBED ABOVE. C IERR=3 IF M IS UNEVEN OR GREATER THAN ISIZE OR IF THE HOST MACHINE C PERFORMS SIGNED MAGNITUDE ARITHMETIC AND M IS UNEVEN OR C GREATER THAN ISIZE-2. C IERR=4 IF THE LEADING ISIZE-M BITS OF ANY ELEMENT OF X ARE NOT ZERO. C IERR=0 OTHERWISE. FX AND PRIME ARE RETURNED UNCHANGED IF IERR C IS RETURNED NON-ZERO. C C NOTES: THIS ALGORITHM IS PRESENTED HERE AS A GENERAL EVALUATION C SCHEME AND THE FORTRAN CODE IS INTENDED TO SERVE AS A MODEL C FOR CONSTRUCTING A SIMILAR SYSTEM UTILITY FUNCTION. C TO PROGRAM AN EFFICIENT RE-ENTERABLE ASSEMBLER VERSION OF C THIS ALGORITHM TO BE USED AS PART OF A SYSTEM PASSWORD C AUTHORIZATION SCHEME THE FOLLOWING IS RECOMMENDED... C C ---- CHOOSE M=ISIZE (EVEN FOR SIGNED MAGNITUDE MACHINES) AND C USE EQUIVALENT MACHINE INSTRUCTIONS IN PLACE OF SUBROUTINES C LADD, LSUB, AND KOMP. ELIMINATE SUBROUTINES SPLIT, JOIN, C ISIGNM AND ALL REFERENCES TO COMMON BLOCK /MACHIN/. C IN SUBROUTINE MPMLT CUT THE PRECISION FROM 2*NP TO NP BY USING C THE EQUIVALENT OF A LOGICAL MULTIPLY MACHINE INSTRUCTION. C C ---- CHOOSE A FAMILY OF ENCIPHERING POLYNOMIALS AND FACTOR IT C BEFORE EVALUATION. FOR EXAMPLE, IF N(1)=2*N(2)-1, N(3)=3, C N(4)=2, N(5)=1, N(6)=0, NN=6, THEN P(X) FACTORS TO... C X**(N(2)-1)*(X*(X**(N(2)-1)+ACOEFF(1,1))) + C X*(ACOEFF(1,4)+X*(ACOEFF(1,3)+ACOEFF(1,2)*X)). C THIS REDUCES THE NUMBER OF CALLS TO EXPP THUS SHARPLY REDUCING C COMPUTING TIME. C C ---- WHEN PRECISION IS CONSTANT, WORKSPACE AND RELATED PARAMETERS C BECOME CONSTANTS AND THE PROGRAM IS CONSIDERABLY SIMPLIFIED. C ONCE THIS IS DONE PASS ALL WORKSPACE VARIABLES IN COMMON AND C REMOVE ALL DIMENSION PARAMETERS AND WORKSPACE PARAMETERS FROM C THE ARGUMENT/PARAMETER LISTS. THIS SUBSTANTIALLY REDUCES C LINKAGE TIME. C C WRITTEN BY H. D. KNOBLE, PENN STATE UNIVERSITY COMPUTATION CENTER, C JUNE 1977. C======================================================================= INTEGER X(NP),FX(NP),Q,A,PRIME(NP),N(NN),ACOEFF(IDIMA,NN),WORK(NW) COMMON/MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK C-----------CHECK PARAMETERS FOR DOMAIN ERRORS. IERR=0 IF(NP.LE.1.OR.NN.LE.0) IERR=1 IF(IDIMA.LT.NP.OR.NW.LT.12*NP+7)IERR=1 IF(A.LE.0.OR.A.GT.MAXPOS) IERR=2 DO 1 I=1,NN IF(N(I).LT.0) IERR=1 1 CONTINUE C------------CHECK HOST MACHINE ASSUMPTIONS ASSUMING MINNEG IS CORRECT. IF((MOD(M,2).NE.0).OR.(M.GT.ISIZE).OR. * (M.GT.ISIZE-2.AND.MINNEG.EQ.0)) IERR=3 C------------CHECK IF THE LEADING (ISIZE-M) BITS OF THE WORDS OF X C ARE NON-ZERO. IF(ISIZE.EQ.M) GO TO 3 DO 2 I=1,NP IF(X(I).GE.MASK) IERR=4 2 CONTINUE 3 IF(IERR.NE.0) RETURN C-----------FOR Q=NP*M, COMPUTE P=2**Q-A. NP1=NP+1 NP2=NP*2 DO 4 I=1,NP2 4 WORK(I)=0 WORK(NP2)=A CALL MPSUB(WORK,WORK(NP1),PRIME,NP) C----------INITIALIZE N2=NP2+1 IWE=10*NP+7 IWM=8*NP+7 IWA=3*NP+2 C----------EVALUATE F(X) ONE TERM AT A TIME. DO 5 I=1,NN CALL EXPP(X,N(I),WORK(NP1),PRIME,WORK(N2),NP,IWE) CALL MULTP(ACOEFF(1,I),WORK(NP1),WORK(NP1),PRIME,WORK(N2),NP,IWM) 5 CALL ADDP(WORK,WORK(NP1),WORK,PRIME,WORK(N2),NP,IWA) C--------------SET RESULT = F(X). FX AND X MAY BE SAME LOCATIONS. DO 6 I=1,NP 6 FX(I)=WORK(I) RETURN END SUBROUTINE EXPP(X,K,Y,P,WORK,NP,NW) 00001920 C========COMPUTE Y=X**K MOD P FOR P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS,X,Y,P OF PRECISION NP. WORK IS A C WORKSPACE ARRAY OF NW=NP*10+7 WORDS. C SEE KNUTH ALGORITHM 4.6.3-A. C INTEGER X(NP),Y(NP),P(NP),WORK(NW) NPP1=NP+1 IWM=NP*8+7 N=K C---------Y=1, Z=X. DO 1 I=1,NP Y(I)=0 1 WORK(I)=X(I) Y(NP)=1 C---------X**0=1 FOR X .GE. 0 . IF(N.LT.1) RETURN 2 L=N N=N/2 IF(MOD(L,2).EQ.0) GO TO 5 CALL MULTP(Y,WORK,Y,P,WORK(NPP1),NP,IWM) IF(N.EQ.0) RETURN 5 CALL MULTP(WORK,WORK,WORK,P,WORK(NPP1),NP,IWM) GO TO 2 END SUBROUTINE MULTP(R,S,RS,P,WORK,NP,NW) 00002170 C=========COMPUTE RS=R*S MOD P FOR P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS R,S,RS,P OF PRECISION NP. WORK C IS A WORKSPACE ARRAY OF NW=NP*8+7 WORDS. INTEGER R(NP),S(NP),RS(NP),P(NP),WORK(NW) NP2=2*NP IW1=NP2+1 IW2=IW1+NP2 IW3=IW2+NP2 CALL MPMLT(R,S,WORK,WORK(IW1),WORK(IW2),WORK(IW3),NP,NP, * NP2,NP2,NP2,2*NP2) CALL MOD2Q(WORK,RS,P,WORK(IW1),NP2,NP,NP*6+7) RETURN END SUBROUTINE ADDP(R,S,RPS,P,WORK,NP,NW) 00002310 C=========COMPUTE RPS=R+S MOD P FOR P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS R,S,RPS,P OF PRECISION NP. C WORK IS A WORKSPACE ARRAY OF NW=NP*3+2 WORDS. INTEGER R(NP),S(NP),RPS(NP),P(NP),WORK(NW) NPP1=NP+1 CALL MPADD(R,S,WORK,NP,NPP1) CALL MODQ1(WORK,RPS,P,WORK(NPP1+1),NPP1,NP,2*NP+1) RETURN END SUBROUTINE MOD2Q(W,R,P,WORK,NP2,NP,NW) 00002410 C=========COMPUTE R=W MOD P FOR 2Q-BIT W, P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS W OF PRECISION NP2=NP*2, C AND R,P OF PRECISION NP. WORK IS A WORKSPACE ARRAY OF C NW=NP*6+7 WORDS. INTEGER W(NP2),R(NP),P(NP),WORK(NW),U,A,JP COMMON /MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK NPP1=NP+1 NPM2=NPP1*2 U=1 JP=U+NPP1 A=JP IWU=JP+NPP1 IW1=IWU+1 IW2=IW1+2 IW3=IW2+NP2 C-----------FORM U=A*W1. WORK(A)= MASK-P(NP) CALL MPMLT(WORK(A),W,WORK(U),WORK(IW1),WORK(IW2),WORK(IW3), * 1,NP,NPP1,2,NP2,NPM2) C---------------USE KNUTH'S THEOREM 4.3.1.B AND LEADING DIGIT THEOREM C TO SOLVE FOR J SUCH THAT J*P.LE. U .LT.(J+1)*P. R(1)=WORK(U)+1 CALL MPMLT(R(1),P,WORK(JP),WORK(IW1),WORK(IW2),WORK(IW3), * 1,NP,NPP1,2,NP2,NPM2) IF(KOMP(WORK(JP),WORK(U),NPP1).LE.0) GO TO 1 CALL MPMLT(WORK(U),P,WORK(JP),WORK(IW1),WORK(IW2),WORK(IW3), * 1,NP,NPP1,2,NP2,NPM2) C----------------COMPUTE U MOD 2**Q - JP MOD 2**Q. 1 CALL MPSUB(WORK(U+1),WORK(JP+1),R,NP) CALL ADDP(W(NPP1),R,R,P,WORK(IW1),NP,NPM2+NP) RETURN END SUBROUTINE MODQ1(S,Y,P,WORK,NP1,NP,NW) 00002740 C===========COMPUTE Y=S MOD P FOR (Q+1)-BIT S AND P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS S OF PRECISION NP1 AND C Y,P OF PRECISION NP. WORK IS A WORKSPACE OF NW=2*NP+1 C WORDS. INTEGER S(NP1),Y(NP),P(NP),WORK(NW) COMMON/MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK CALL MODQ(S(2),Y,P,NP) C-------------S(1), THE HIGH-ORDER PART OF S, IS EITHER ZERO OR ONE. IF(S(1).EQ.0) RETURN C-------------RECOVER A=2**Q-P DO 1 I=1,NP 1 WORK(I)=0 WORK(NP)=MASK-P(NP) C-------------COMPUTE SUM=MODQ(S MOD 2**Q)+A*S(1). CALL MPADD(Y,WORK,WORK(NP1),NP,NP1) CALL MODQ(WORK(NP1+1),Y,P,NP) RETURN END SUBROUTINE MODQ(X,Y,P,NP) 00002930 C===========COMPUTE X MOD P FOR Q-BIT X AND P=2**Q-A FOR UNSIGNED C MULTIPRECISION INTEGERS X,Y,P OF PRECISION NP. INTEGER X(NP),Y(NP),P(NP) IF(KOMP(X,P,NP).GE.0) GO TO 2 C------------Y=X. DO 1 I=1,NP 1 Y(I)=X(I) RETURN C---------Y=X-P 2 CALL MPSUB(X,P,Y,NP) RETURN END SUBROUTINE MPMLT(IX,IY,IZ,U,V,W,NX,MY,NPM,N2,M2,NPM2) 00003060 C===========UNSIGNED MULTIPRECISION INTEGER MULTIPLICATION C IZ=IX*IY. NX,MX,NPM ARE THE PRECISIONS OF THE UNSIGNED C INTEGERS IX,IY,IZ RESPECTIVELY. U,V,W ARE C WORKSPACES OF N2=NX*2, M2=MY*2, AND MPN2=(NX+MY)*2 WORDS C RESPECTIVELY. NO INTEGER OVERFLOWS ARE GENERATED. C C SEE KNUTH'S ALGORITHM 4.3.1 M. C INTEGER IX(NX),IY(MY),IZ(NPM),U(N2),V(M2),W(NPM2) INTEGER T(2),K,SHIFT1,V0 COMMON /MACHIN/ ISIZE,M,MAXPOS,MINNEG,MASK SHIFT1=2**(M/2-1) C----------DOUBLE THE PRECISION SINCE FORTRAN CANNOT ACCESS THE C PRODUCT OF TWO M-BIT WORDS. DO 10 I=1,NX I2=I*2 10 CALL SPLIT(IX(I),U(I2-1),U(I2)) DO 11 I=1,MY I2=I*2 11 CALL SPLIT(IY(I),V(I2-1),V(I2)) C----------BEGIN KNUTH'S ALGORITHM. DO 1 I=1,NPM2 1 W(I)=0 J=M2 C----------THE PROBABILITY OF V(J)=0 IS SMALL. SKIP STEP 2. 3 I=N2 K=0 4 IPJ=I+J C----------COMPUTE T=U(I)*V(J)+W(I+J)+K WITHOUT INTEGER OVERFLOW. V0=MOD(V(J),SHIFT1) CALL LADD(U(I)*V0, (V(J)-V0)*U(I), IZ) CALL LADD(IZ(2), W(IPJ)+K, T) C-----------COMPUTE K = FLOOR(T/2**(M/2)), W(IPJ)=T MOD 2**(M/2), C THE HIGH AND LOW-ORDER HALFS OF T(2). CALL SPLIT(T(2),K,W(IPJ)) 5 I=I-1 IF(I.GT.0) GO TO 4 W(J)=K 6 J=J-1 IF(J.GT.0) GO TO 3 C----------CONVERT THE RADIX 2**(M/2) DIGITS TO RADIX 2**M DIGITS. DO 12 I=1,NPM I2=I*2 12 IZ(I)=JOIN(W(I2-1),W(I2)) RETURN END SUBROUTINE MPADD(U,V,W,N,NP1) 00003530 C=============UNSIGNED MULTIPRECISION INTEGER ADDITION, C W=U+V. N IS THE PRECISION OF UNSIGNED NUMBERS, IX, IY. C W IS THE UNSIGNED NP1-WORD PRODUCT. C C SEE KNUTH'S ALGORITHM 4.3.1-A. INTEGER U(N),V(N),W(NP1),SUM1(2),SUM2(2) 1 J=N K=0 2 CALL LADD(U(J),V(J),SUM1) CALL LADD(SUM1(2),K,SUM2) K=SUM1(1)+SUM2(1) W(J+1)=SUM2(2) 3 J=J-1 IF(J.GT.0) GO TO 2 W(1)=K RETURN END SUBROUTINE MPSUB(U,V,W,N) 00003710 C===========UNSIGNED MULTIPRECISION INTEGER SUBTRACTION, C W=U-V IGNORING POSSIBLE BORROW IF U IS LESS THAN V. C N IS THE PRECISION OF UNSIGNED NUMBERS U,V,W. C C SEE KNUTH'S ALGORITHM 4.3.1-S. INTEGER U(N),V(N),W(N),DIF1(2),DIF2(2) 1 J=N K=0 2 CALL LSUB(U(J),V(J),DIF1) CALL LSUB(DIF1(2),K,DIF2) K=DIF1(1)+DIF2(1) W(J)=DIF2(2) 3 J=J-1 IF(J.GT.0) GO TO 2 RETURN END FUNCTION KOMP(U,V,N) 00003880 CXXXXXXXXX FUNCTION TO COMPARE UNSIGNED MULTIPRECISION C INTEGERS U,V OF PRECISION N. KOMP=-1 IF U IS LESS THAN V, C +1 IF U IS GREATER THAN V, AND 0 IF U=V. C C SEE KNUTH PROBLEM 4.3.1-11 INTEGER U(N),V(N) DO 1 J=1,N C----------IF HIGH-ORDER BITS DIFFER U .NE. V. IF(ISIGN(1,U(J)).NE.ISIGN(1,V(J))) GO TO 2 C----------ELSE IF U IS LESS THAN V, THEN KOMP=-1 IF(U(J).LT.V(J)) GO TO 3 C----------ELSE IF U IS GREATER THAN V THEN KOMP=+1. IF(U(J).GT.V(J)) GO TO 4 1 CONTINUE C----------U=V. KOMP=0 RETURN C----------THE OPERAND WITH HIGH-ORDER BIT ON IS GREATER. 2 IF(U(J).LT.V(J))GO TO 4 3 KOMP=-1 RETURN 4 KOMP=+1 RETURN END FUNCTION JOIN(HALF1,HALF2) 00004130 CXXXXXXXXXX ROUTINE TO CONVERT TWO RADIX 2**(M/2) DIGITS HALF1, C HALF2 TO A RADIX 2**M DIGIT, JOIN. THAT IS THIS C ROUTINE FORMS A M-BIT NUMBER WITH THE HIGH-ORDER ISIZE-M C BITS EQUAL TO ZERO, THE M/2 HIGH-ORDER BITS BEING THE M/2 C LOW-ORDER BITS OF HALF1, AND THE M/2 LOW-ORDER BITS BEING C THE M/2 LOW-ORDER BITS OF HALF2. C THIS IS DONE WITHOUT CAUSING AN INTEGER OVERFLOW AND C ASSUMING M IS EVEN, AND HIGH-ORDER (ISIZE-M)+M/2 BITS OF C EACH HALF ARE ZERO. INTEGER HALF1,HALF2,ISUM(2),SHIFT1 COMMON/MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK SHIFT1=2**((M/2)-1) C------------FORM JOIN=HALF2+(2**(M/2))*HALF1 IX=HALF1*SHIFT1 CALL LADD(IX,IX,ISUM) JOIN=ISUM(2)+HALF2 RETURN END SUBROUTINE SPLIT(IX,HALF1,HALF2) 00004320 CXXXXXXXXXXXX ROUTINE TO CONVERT A RADIX 2**M DIGIT, IX TO TWO C RADIX 2**(M/2) DIGITS HALF1,HALF2. THAT IS THIS C ROUTINE SPLITS THE M-BIT NUMBER IX INTO HIGH-ORDER HALF C HALF1 AND LOW-ORDER HALF2 WITHOUT CAUSING AN INTEGER C OVERFLOW. THE HIGH-ORDER ISIZE-M BITS OF IX ARE ASSUMED C TO BE ZERO. INTEGER X,HALF1,HALF2,SHIFT,ISUM(2) COMMON/MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK SHIFT=2**(M/2) X=IX IF(X.GE.0.AND.X.LE.MAXPOS) GO TO 2 C------------REMOVE SIGN BIT BEFORE DIVISION. CALL LADD(MINNEG,X,ISUM) HALF1=ISUM(2)/SHIFT+SHIFT/2 C------------FORM X MOD 2**(M/2) HALF2=MOD(ISUM(2),SHIFT) RETURN 2 HALF1=X/SHIFT HALF2=MOD(X,SHIFT) RETURN END SUBROUTINE LSUB(IU,IV,IW) 00004540 CXXXXXXXXX SUBROUTINE TO PERFORM A LOGICAL SUBTRACT OF 2 M-BIT C NUMBERS IU, IV. IW(2) IS THEIR LOGICAL DIFFERENCE, C IW(1) IS +1 IF THERE WAS A BORROW, ZERO OTHERWISE. C SEE SUBROUTINE LADD. INTEGER IU,IV,IX,IY,IW(2),ISUM(2) COMMON /MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK IX=IU IY=IV IW(1)=0 C----------SORT OUT SPECIAL CASES FOR LADD. IF(IX.EQ.IY) GO TO 3 IF(IY.EQ.0) GO TO 1 IF(IY.EQ.MINNEG) GO TO 2 CALL LADD(IX,MASK-IY,ISUM) IF(ISUM(1).EQ.0) IW(1)=1 IW(2)=ISUM(2) RETURN 1 IW(2)=IX RETURN 2 CALL LADD(IX,IY,ISUM) IF(ISUM(1).EQ.0) IW(1)=1 IW(2)=ISUM(2) RETURN 3 IW(2)=0 RETURN END SUBROUTINE LADD(IU,IV,IW) 00004810 CXXXXXXXXXXXXX SUBROUTINE TO PERFORM A LOGICAL ADD OF THE LOW-ORDER C M BITS OF TWO COMPUTER WORDS, IU, IV. IW(2) IS THEIR C LOGICAL SUM, IW(1) IS THE CARRY BIT VALUE. C NO INTEGER OVERFLOWS OCCUR EVEN IF M=ISIZE. C EXAMPLE - 0111 LADD 1011 YIELDS IW(1)=1, IW(2)=0010, C FOR M=4, MAXPOS=0111, MINNEG=1000, MASK=0 OR 10000. C INTEGER IX,IY,U,IV,IW(2),DIF,NEG,POS,ISUM,BIG,SMALL,ISX,ISY,ISD COMMON/MACHIN/ISIZE,M,MAXPOS,MINNEG,MASK IX=IU IY=IV IW(1)=0 C-------------ELIMINATE ZERO OPERANDS. IF(IX.EQ.0.OR.IY.EQ.0) GO TO 2 C-------------DETERMINE IF IX,IY IS ++, --, OR +- ISX=ISIGNM(IX) ISY=ISIGNM(IY) IF(ISX.NE.ISY) GO TO 5 IF(ISX.LT.0) GO TO 3 C-------------++ (BOTH OPERANDS HAVE HIGH-ORDER BITS OFF) 1 DIF=MIN0(IX,IY)-(MAXPOS-MAX0(IX,IY)) ISD=ISIGNM(DIF) IF(ISD.LE.0) IW(2)=IX+IY IF(ISD.GT.0) IW(2)=(MINNEG+DIF)-1 RETURN C------------AT LEAST ONE OPERAND IS ZERO. 2 IW(2)=IX+IY RETURN C------------ -- (BOTH OPERANDS HAVE HIGH-ORDER BITS ON) 3 IW(1)=1 BIG=MAX0(IX,IY) SMALL=MIN0(IX,IY) DIF=0 IF(SMALL.NE.MINNEG) DIF=MINNEG-SMALL DIF=BIG-DIF IF(ISIGNM(DIF).GE.0) GO TO 4 IW(2)=((MAXPOS+DIF)+1)-MASK RETURN 4 IW(2)=IX+IY RETURN C------------- +- (OPERANDS HAVE OPPOSITE HIGH-ORDER BIT VALUES). 5 IW(2)=IX+IY NEG=IX POS=IY IF(ISIGNM(IX).LT.0) GO TO 6 NEG=IY POS=IX C------------- IF ++(--(MINNEG,NEG),POS) IMPLIES A ++ INTEGER OVERFLOW C THEN THERE WILL BE A CARRY FOR +-. 6 ISUM=((MAXPOS+NEG)+1)-MASK IF(MIN0(ISUM,POS).LE.(MAXPOS-MAX0(ISUM,POS))) RETURN IW(2)=IW(2)-MASK IW(1)=1 RETURN END FUNCTION ISIGNM(IX) 00005370 CXXXXXXXX FUNCTION TO RETURN 0 IF IX=0, 1 IF 0