	SUBROUTINE MULMUL (PT1,PT2,RETCD,ENTRY)
C COPYRIGHT (C) 1983 GLENN EVERHART
C PERMISSION IS GIVEN TO ANYONE TO USE, DISTRIBUTE, OR COPY THIS
C PROGRAM FREELY BUT NOT TO SELL IT COMMERICALLY.
	INCLUDE 'VKLUGPRM.FTN'
C	PARAMETER RRW = 32
C	PARAMETER RCL = 32
C RRW=MAX REAL ROWS
C RCL=MAX REAL COLS
C RRW MUST BE 1 LARGER TO HANDLE 1ST 27 VARIABLES IN AVBLS
C VBLS AND TYPE DIMENSIONED RRW,RCL
C **************************************************
C *                                                *
C *  SUBROUTINE  MULMUL (PT1,PT2,RETCD,ENTRY)      *
C *                                                *
C **************************************************
C
C
C  MULTIPLE PRECISION MULTIPLY ROUTINES
C
C  	ARGUMENT LIST IS (PT1,PT2,RETCD,ENTRY) WHERE OPERATION
C  	OF * IS PERFORMED AS FOLLOWS:
C
C		STACK1(,PT1) GETS VALUE OF STACK1(,PT1)*STACK2(,PT2)
C
C  NOTE:  STACK2 IS NOT CLEANED UP BY THE OPERATION
C
C  RETCODES	VALUE OF RETCD	MEANING
C
C			1	NORMAL
C			2	ERROR (OVERFLOW)
C
C
C
C ENTRY SPECIFIES BASE: 8, 10, OR 16
C
C
C MODIFY CODES: M3, M4, M10
C
C
C
C
C
C
C MULMUL CALLS ERRMSG TO PRINT ERROR MESSAGES.
C
C
C
C
C MULMUL IS CALLED BY CALBIN
C
C
C
C
C
C  VARIABLE    USE
C
C    BASE    BASE UNDER WHICH OPERATION IS PERFORMED.
C    CARRY   CARRY INTO NEXT POWER OF BASE.
C    ENTRY   SPECIFIES BASE IN ARGUMENT OF SUBTRACTION.
C    I,J     TEMPORARY VALUES.
C    PT1     STACK 1 POINTER TO OPERAND 1.
C    PT2     STACK 2 POINTER TO OPERAND 2.
C    PSUM    VECTOR THAT SUMS PARTIAL PRODUCTS.
C    RETCD   RETURN CODE: 1=O.K.,  2=ERROR.
C    TEMP    HOLDS INTEGER*4 TEMPORARY VALUES.
C    ZL1     POINTS TO HIGH ORDER NON-ZERO DIGIT OF OPERAND 1.
C    ZL2     POINTS TO HIGH ORDER NON-ZERO DIGIT OF OPERAND 2.
C
C
C
C
C
C ************************
C *++++++++++++++++++++++*
C *+                    +*
C *+     WARNING        +*
C *+                    +*
C *++++++++++++++++++++++*
C ************************
C
C    IF THE NUMBER OF MULTIPLE PRECISION DIGITS IS INCREASED TO N
C OR INTEGER*4 NOT AVAILABLE ETC., BE CERTAIN THAT 2*N*(BASE-1)**2
C CAN BE HELD BY EACH ELEMENT OF PSUB, TEMP, AND CARRY. IF NOT, THEN
C THE REDUCTION TO CANONICAL FORM IN PSUM MUST BE DONE AFTER EACH
C PARTIAL PRODUCT IS ADDED IN.
C
C
C
C
C MODIFIED 4-DEC-1979 P.B.
C CHANGED LINES 460 AND 510 TO USE TEMP TO FORCE EVALUATION OF PRODUCT
C AS INTEGER*4. THIS FIXED BUG THAT PREVENTED PROPER EVALUATION OF
C 000000000000000F*0F
C
C
C
C
C
C
C	SUBROUTINE MULMUL (PT1,PT2,RETCD,ENTRY)
C
	INTEGER*4 PSUM(19)
	INTEGER*4 BASE,TEMP,ZL1,ZL2,CARRY
C
C
	InTEgeR*4 ST1TYP(40),ST2TYP(40)
	InTEgeR*4 RETCD,ENTRY
	InTEgeR*4 ST1PT,ST2PT,ST1LIM,ST2LIM
	InTEgeR*4 PT1,PT2
	InTEgeR*4 I
C
	LOGICAL*1 STACK1(20,40),STACK2(20,40)
C
	COMMON /STACK/STACK1,STACK2,ST1PT,ST2PT,ST1TYP,ST2TYP,
     ;          ST1LIM,ST2LIM
C
C
C
C
C
	BASE=ENTRY
C
C
C  DETERMINE SIGN
	RETCD=1
	TEMP=STACK1(20,PT1)+STACK2(20,PT2)
	STACK1(20,PT1)=TEMP-TEMP/2*2
C
C
C  ZERO PARTIAL PRODUCT SUM VECTOR
	DO 200 I=1,19
200	PSUM(I)=0
C
C
C  FIND WHERE FIRST NON-ZEROES ARE
	DO 210 I=1,19
	ZL1=20-I
	IF (STACK1(ZL1,PT1).NE.0) GOTO 220
210	CONTINUE
	ZL1=0
220	DO 250 I=1,19
	ZL2=20-I
	IF (STACK2(ZL2,PT2).NE.0) GOTO 300
250	CONTINUE
	ZL2=0
C
C
C
300	IF (ZL1.NE.0.AND.ZL2.NE.0) GOTO 400
C
C
C ****************************************************
C ****** ONE OF THE FACTORS IS 0 SO ANSWER IS 0 ******
C ****************************************************
	DO 310 I=1,20
310	STACK1(I,PT1)=0
	RETURN
C
C
C  MAKE SURE THINGS AREN'T TOO BIG ALREADY
400	IF (ZL1+ZL2.LE.20) GOTO 450
C
C
C
C
C **** ERROR RETURN ****
C
C RESULT IS GREATER THAN 99 DIGITS
410	RETCD=2
	CALL ERRMSG (22)
	RETURN
C
C
C  DETERMINE THE MULTIPLIER (FACTOR WITH FEWEST DIGITS TO SPEED OPERATION)
450	IF (ZL1.GT.ZL2) GOTO 500
C
C
C ******************************************
C ******  PERFORM THE MULTIPLICATION  ******
C ******************************************
	DO 460 I=1,ZL1
	DO 460 J=1,ZL2
C FORCE PRODUCT TO USE INTEGER*4 IN CASE 0F*0F=E1 AND NOT ABLE
C TO HOLD IN A SIGNED BYTE
	TEMP=STACK1(I,PT1)
460	PSUM(I+J-1)=PSUM(I+J-1)+TEMP*STACK2(J,PT2)
	GOTO 600
500	DO 510 I=1,ZL2
	DO 510 J=1,ZL1
C FORCE PRODUCT TO USE INTEGER*4 IN CASE TRYING TO CALCULATE THINGS LIKE
C 0F*0F=E1 WHICH DOES NOT FIT INTO A SIGNED BYTE
	TEMP=STACK1(J,PT1)
510	PSUM(I+J-1)=PSUM(I+J-1)+TEMP*STACK2(I,PT2)
C
C
C
C ************************************************************
C ****** REDUCE ANSWER TO STANDARD CANONICAL FORM WHERE ******
C ****** POWERS OF THE BASE ARE TIMES A DIGIT LESS      ******
C ****** THAN THE BASE.                                 ******
C ************************************************************
600	CARRY=0
	DO 650 I=1,19
	TEMP=PSUM(I)+CARRY
	CARRY=TEMP/BASE
650	STACK1(I,PT1)=TEMP-CARRY*BASE
	IF (CARRY.EQ.0) RETURN
	GOTO 410
	END
