|
Listing 4. MULTIPLE PRECISION MATH WORDS
This math package is derived from Dr. S. Y. Tang's work on multiple precision math distributed with the F-PC system, Version 2.25. It is adapted to the cmForth for RTX2000, and is verified on a single board computer based on the Forth Kit from Indelko, Sweden. The major difference, may be improvement, is that the components are represented by full 16 bit integers, instead of four decimal digits as used by Dr. Tang. Using decimal digits has the advantage that the internal representation can be converted to external (human readable) representation very conveniently. However, math operations on decimal components are often cumbersome and time consuming. As the primary goal of this implementation is to achieve th highest possible performance, binary components are more suitable. An additional advantage is the number size is slightly increased. Using binary component for internal representation requires that numbers can be entered in decimal digital formats and that they are converted to decimal format for displaying. The I/O processes are never time critical and special format conversion words are defined to convert number between the internal and external formats. A multiple precision number consists of a sign integer, a component integer, and followed by n components. It has a length of n+2 16 bit words. The components are arranged with the least significant component at the lowest memory. For positive numbers, the sign is 1, and for negative numbers the sign is -1. Zero is represented by an 1 in the sign integer and a zero in the component integer, without any component. The components represent the magnitude of a number. Address of a multiple precision number points to its component integer. The sign integer is at address-2, and the least significant component starts at address+2. The absolute value of a number is thus similar to a 'length string', in which the first 16 bit integer contains the length (number of components) of the multiple number. Since RTX2000 has an integral multiplier and very efficient math step instructions, this package tries to take the best advantage of these resources. The performance is more than 100 times faster than the similar implementation on a regular personal computers.
( number arrays, 07may89cht)
: CELLS ( component -- byte ) 2* ; : CELL! ( n addr -- ) 2 + ! ; : CELL@ ( addr -- n ) 2 + @ ; : CELLS! ( n addr component -- ) 2* + ! ; : CELLS@ ( addr component -- n ) 2* + @ ; : CELLS+! ( n addr component -- ) 2* + +! ;
400 CONSTANT MAXSIZE 400 bytes per number, about 1000 digits.
: M# ( -- ) Define name numbers. CREATE HERE MAXSIZE DUP ALLOT Allocate 400 bytes in dictionary. ERASE Clear to zero. DOES R> Retrieve parameter field address. 2 + ; Return address of component integer.
M# AA M# BB M# CC M# DD M# EE General purpose numbers M# ANS M# AA2 M# BB2 Accumulators M# PROD M# REM M# QUOT M# TMP Multiply and divide registers M# GUESS M# SQRT Square root registers M# NN Input number register
( field operators)
: SIGN@ ( aa -- sign) Return the sign integer of a number. 2 - @ ; : SIGN! ( sign aa -- ) Change the sign integer of a number. 2 - ! ; : COMP@ ( aa -- comp) Return the component integer. @ ; : COMP! ( n aa -- ) Store the component integer. ! ;
: S>M ( n aa -- ) Initiate a number to an integer n. 1 OVER SIGN! Make it positive. 1 OVER COMP! Need only 1 component. 1 CELLS! ; Store the component.
M# M1 1 M1 S>M M# M2 2 M2 S>M Most useful numbers, 1 and 2.
: M! ( aa bb -- ) Copy number aa to bb. 2 - >R Destination address DUP @ Number of components SWAP 2 - SWAP Source address R> SWAP 2 + MOVE ; Copy number
: :=0 ( aa -- ) Initiate a number to 0. 1 OVER SIGN! Sign of 0 must be positive. MAXSIZE 2 - ERASE ; Everything else are zero.
: MABS ( aa -- ) Make aa positive. 1 SWAP SIGN! ;
( get number, 15may89cht) Temporary number I/O words. Display the components as integers. Will be replace by MGET and M. when tools are available to define them.
: GET# ( aa n -- ) Get the following n integers and using them to build a number. DUP 0< Negate number has n<0. IF -1 ELSE 1 THEN >R OVER R> SWAP SIGN! Put in the correct sign. ABS 2DUP SWAP COMP! abs(n) is the component integer. 1 - FOR 32 WORD NUMBER Get next integer from input stream, OVER I 1 + CELLS! convert and store component. NEXT DROP ;
: .# ( aa -- ) Primitive number output. Display sign and then the components. DUP SIGN@ -1 = IF 45 EMIT ELSE SPACE THEN Print optional - sign. DUP COMP@ ?DUP IF 1 - Print components. FOR DUP I 1 + CELLS@ U. NEXT ELSE 48 EMIT Print a 0 if the number is zero. THEN DROP ;
( utility, 21may89cht)
: M0= ( aa -- f) @ 0= ;
: M0< ( aa -- f) DUP @ IF SIGN@ 0< Test sign. ELSE DROP 0 Test 0. THEN ;
: M0> ( aa -- f) DUP @ IF SIGN@ 0 > Test sign. ELSE DROP 0 Test 0. THEN ;
: M0<= ( aa -- f) DUP SIGN@ 0< Test sign. SWAP @ 0= OR ; Test 0.
: M0>= ( aa -- f) DUP SIGN@ 0 > SWAP @ 0= OR ;
: MNEGATE ( aa -- ) Negate a number. DUP @ IF DUP SIGN@ NEGATE Change sign, if number is not 0. SWAP SIGN! ELSE 1 SWAP SIGN! Make sure 0 has a positive sign. THEN ;
: #COMP ( aa -- aa ) Eliminate leading 0 components in a number by adjusting component integer. BEGIN DUP COMP@ 2DUP CELLS@ 0= Test the most significant component. SWAP 0 > AND WHILE -1 OVER +! Eliminate it if it is 0. REPEAT Continue till MSC is not zero. DUP @ 0= If result is 0, IF 1 OVER SIGN! THEN make sure it is positive. ;
( madd and msub, 08may89cht) Primitive add and subtract operators.
: MADD ( aa bb n -- ) Add n components in aa and bb. Sum is stored in ANS. aa and bb are assumed to be positive. DUP >R ANS ! ANS will have n components. ANS 2 + 4 G! Save destination address in MD. 2 + 6 G! Address of source 2. Saved in SR. 2 + Address of source 1. Left on stack. R> 1 - Retrieve n, which must be > 0. 0 + Clear carry. FOR 2 @+ SWAP Get a component from aa. 6 G@ 2 @+ Get a component from bb. 6 G! Save address to bb back in SR. +c Add tow components with carry. 4 G@ 2 !+ Store sum to ANS. 4 G! Save address to ANS back in MD. NEXT 0 0 +c Last carry? IF 1 4 G@ ! Yes, add another component to ANS. 1 ANS +! Extend component integer of ANS. THEN DROP ; Done.
: MCOMPLEMENT ( aa -- ) Complement aa. The resulting number and aa have a sum of n zero components and a final carry of 1. 2 @+ SWAP Fetch the component integer. ?DUP If it is not 0, complement all the IF 1 - components. FOR 0 @+ SWAP Get a component. 0SWAP-c Subtract it from zero with borrow. SWAP 2 !+ Put complement back. NEXT Repeat for all n components. THEN DROP ;
: MSUB ( aa bb n -- ) Subtract n components in bb from aa. Signs are ignored. DUP >R ANS ! Result in ANS will have n components. ANS 2 + 4 G! Save ANS pointer in MD. 2 + 6 G! Save bb pointer in SR. 2 + Leave aa pointer on stack. R> 1 - Run through n components. FOR Borrow is cleared during - operation. 2 @+ SWAP Fetch a component in aa. 6 G@ 2 @+ 6 G! Fetch a component in bb. -c Subtract with borrow. 4 G@ 2 !+ 4 G! Save results in ANS. NEXT 0 0 -c If last borrow is set, aa<bb. IF ANS MCOMPLEMENT Complement results. ANS MNEGATE Negate it. THEN DROP aa>bb. No problem. ANS #COMP DROP ; Adjust component integer if necessary.
( new math, 08may89cht)
: MEXTEND ( aa n -- ) Clear n components above aa, to be sure that no garbage is left there. DUP 0 > Skip if n is 0. IF OVER DUP @ Get the component integer. 1 + 2* + Address above MS component. OVER 2* ERASE Clear n components there above. SWAP +! Update the component integer. ELSE 2DROP THEN ;
: MALIGN ( aa bb -- aa bb n f) Make aa and bb the same length n. Return true if signs of aa and bb are the same. OVER @ OVER @ 2DUP - aa and bb have the same length? ?DUP IF DUP 0< Not same length. aa<bb? IF NEGATE Absolute difference. SWAP >R >R DROP Save length of bb and length difference. OVER Copy aa for MEXTEND. ELSE SWAP-DROP aa>bb. SWAP >R >R Save length of aa and length difference. DUP Copy bb for MEXTEND THEN R> MEXTEND Extend the shorter number to alignment. ELSE DROP >R Same length. No extension needed. THEN OVER SIGN@ Get sign of aa. DUP ANS SIGN! Make it the sign of ANS. OVER SIGN@ = Compare it with the sign of bb. R> SWAP ; Leave comparison flag as f.
: M+ ( aa bb -- ) Adding aa and bb. Sum is is ANS. DUP M0= If bb=0, ANS=aa. IF DROP ANS M! EXIT THEN MALIGN Align aa and bb. IF MADD Add if signs are the same. ELSE MSUB Subtract if signs are difference. THEN ;
: M- ( aa bb -- ) Subtract bb from aa. Difference is in ANS. DUP M0= If bb=0, ANS=aa. IF DROP ANS M! EXIT THEN MALIGN Align aa and bb. IF MSUB Same sign, subtract. ELSE MADD Different sign, add. THEN ;
( simple math, 15may89cht)
: M1+ ( aa -- ) Increment aa by 1. DUP M1 M+ Add M! to aa. ANS SWAP M! ; Copy sum back to aa.
: M1- ( aa -- ) Decrement aa by 1. DUP M1 M- Subtract M1 from aa. ANS SWAP M! ; Copy difference back to aa.
: M= ( aa bb -- flag) M- ANS M0= ; Subtract then compare with 0. : M> ( aa bb -- flag) M- ANS M0> ; : M< ( aa bb -- flag) M- ANS M0< ; : M>= ( aa bb -- flag) M- ANS M0>= ; : M<= ( aa bb -- flag) M- ANS M0<= ;
( multiply by integer, 17may89cht)
: SM* ( aa n --, product in AA2) Multiply aa by an integer n. OVER @ 0= Test for 0 cases. IF 2DROP AA2 :=0 EXIT THEN DUP 0= IF 2DROP AA2 :=0 EXIT THEN 6 G! Save n in MD. DUP SIGN@ AA2 SIGN! Copy sign. AA2 2 + 4 G! Save pointer to product. 0 SWAP DUP @ DUP AA2 ! Init component number of product. 1 - Loop count. SWAP 2 + SWAP Source component pointer. FOR 2 @+ >R Get a component from source. 6 G@ MULU MLR@ MHR@ Multiply with n. >R +c Add carry from last component. 4 G@ 2 !+ 4 G! Store to product. R> R> Restore carry and source pointer. NEXT DROP 0 +c Carry from last component? ?DUP IF Yes, 4 G@ ! Store MS component. 1 AA2 +! Extend component count. THEN ;
: MACCUMULATE ( pp -- ) Add AA2 to pp in PROD. Used by M*. AA2 @ ?DUP If AA2=0, do nothing. IF 1 - AA2 2 + 4 G! Similar to MADD except using pp. FOR 0 @+ >R 4 G@ 2 @+ 4 G! +c R> 2 !+ NEXT 0 0 +c ?DUP IF 1 SWAP ! ELSE DROP THEN ELSE DROP THEN ;
( m*, 17may89cht)
: M* ( aa bb --, product is PROD) Multiply aa and bb. Use SM* steps. PROD :=0 Init product. DUP M0= Test for 0 cases. IF 2DROP EXIT THEN OVER M0= IF 2DROP EXIT THEN OVER @ OVER @ + Product component count. PROD ! OVER SIGN@ OVER SIGN@ Sign of product. 0< IF NEGATE THEN PROD SIGN! PROD OVER @ 1 - Component count in bb. FOR >R 2 + Pick a component in bb. 2DUP @ SM* Multiply aa with it. R> 2 + DUP MACCUMULATE Accumulate partial product. NEXT Repeat through bb. 2DROP DROP PROD #COMP DROP ; Adjust component count of product.
HEX COMPILER
3 SHIFT 2*c Shift opcodes. 5 SHIFT c2/ : @- E540 SHORT ; Decrement fetch. : !- E5C0 SHORT ; Decrement store. DECIMAL FORTH
( shifts, 18may89cht) Shifts are much cheaper than division.
: M2* ( aa -- ) Left shift aa by one bit. DUP 2 @+ SWAP Get component. ?DUP Do shift only if aa is not 0. IF 1 - 0 2*c DROP Clear carry first. FOR 0 @+ Fetch a component. SWAP 2*c Shift left with carry. SWAP 2 !+ Store component. NEXT 0 2*c Overflow from MS component IF 1 SWAP ! Yes, extend aa. 1 SWAP +! ELSE 2DROP THEN ELSE 2DROP THEN ;
: M2/ ( aa -- ) Right shift aa by one bit. DUP @ Very similar to M2*. 2* OVER + OVER @ ?DUP IF 1 - 0 c2/ DROP FOR 0 @- SWAP c2/ Left shift with carry. SWAP 2 !- NEXT DROP #COMP DROP ELSE 2DROP THEN ;
: MSHL ( aa n -- ) Shift aa left by n bits. ?DUP Repeat M2* n times. IF 1 - FOR DUP M2* NEXT THEN DROP ;
: MSHR ( aa n -- ) Shift aa right by n bits. ?DUP Repeat M2/ n times. IF 1 - FOR DUP M2/ NEXT THEN DROP ;
( divide by integer, 18may89cht) SM* will be used to convert number string to binary, and SM/ will be used to convert binary to number String. Important for number I/O.
: SM/ ( aa n -- rem ) Divide aa by integer n. Quotient in aa. SWAP DUP @ Test for 0. DUP 0= IF SWAP-DROP SWAP-DROP EXIT THEN DUP >R Save component count. 2* + Address of MS component. SWAP 0 SWAP ( aa' 0 n -- ) Get everything in order. R> 1 - Run through aa from MS to LS component. FOR >R >R 0 @- SWAP Get a component. R> I M/MOD 32 bit mixed divide. >R SWAP 2 !- Save quotient. Leave rem on stack. R> R> NEXT DROP SWAP #COMP DROP ; Adjust component count.
: SM+ ( aa n -- ) Add n to aa. OVER 2 @+ Test for 0. SWAP ?DUP IF 1 - 0 + Clear carry. FOR 0 @+ >R Add n to first component. +c Carry over other components. R> 2 !+ 0 SWAP NEXT SWAP 0 +c In rarely occasion, carry over IF 1 SWAP ! the MS component. 1 SWAP +! ELSE 2DROP THEN ELSE ! If aa=0, make aa=n. 1 SWAP ! THEN ;
( m. , 18may89cht) Number must be converted to ASCII string for output to a terminal.
VARIABLE MPLD Output string pointer. HEX F000 CONSTANT MPAD Output string buffer. DECIMAL
: MHOLD ( n -- ) Add n to the output string. MPAD @ MPLD @ - 5 MOD 0= Insert a space after four digits, IF -1 MPLD +! to improve readability. 32 MPLD @ C! THEN -1 MPLD +! Add n to output string. MPLD @ C! ; Decrement string pointer.
: M. ( aa -- ) Display nn in four digit groups. DUP SIGN@ >R Save sign. AA2 M! Copy aa to AA2 for conversion. AA2 Replace aa by AA2. MPAD MPLD ! Init output string pointer. BEGIN DUP BASE @ SM/ Convert one digit. DIGIT MHOLD Add that digit to output string. DUP M0= Repeat until AA2 is exhausted. UNTIL DROP R> 0< Add sign. IF 45 MHOLD THEN MPLD @ Print output string. MPAD OVER - TYPE ;
: MGET ( aa -- ) Convert following string to a number and store it in aa. DUP :=0 Init aa and AA2. AA2 :=0 1 WORD Get the following number string. 1 C@+ OVER 0= Exit if string is null. IF 2DROP EXIT THEN DUP C@ 45 = First character a - sign? IF 1 + >R If so, adjust pointers. 1 - R> -1 and put -1 for sign. ELSE 1 No -, put 1 for sign. THEN >R ( sign) Save the sign. SWAP ?DUP Scan the input string. IF 1 - FOR 1 C@+ >R Get the next character. DUP 32 - Ignore spaces and control characters. IF OVER DUP BASE @ SM* Shift existing number to the left. AA2 SWAP M! -DIGIT Convert next digit. OVER SWAP SM+ Accumulate new digit to the sum. ELSE DROP THEN R> NEXT Repeat to end of string. THEN DROP R> SWAP SIGN! ; Put in the correct sign.
( divide primitives, 19may89cht) Division is pure hard work.
VARIABLE #SHIFT Bits to left justify the divisor. VARIABLE #OFFSET Offset to MS component of divisor. VARIABLE #QUOT Temporary place to store 16 bit quotient.
: (DIVIDE) ( -- ) Dividend in REM, divisor in TMP. A division step, processing one component in the remainder. REM #OFFSET @ + Get 2 MS components in remainder. 2 @+ @ ( dividend) TMP #OFFSET @ + @ Get MS component is divisor. M/MOD DROP Trial divide. DUP #QUOT ! Save trial quotient in #QUOT. TMP SWAP SM* Multiply divisor with trial quotient. REM AA2 M- Subtract results from remainder. ANS REM M! Save difference. ANS M0< If difference is negative, adjust quotient and remainder. IF -1 #QUOT +! As the divisor is left justified, REM TMP M+ no more than two additions ANS REM M! are needed to bring the remainder ANS M0< back to positive. This was IF -1 #QUOT +! suggested by Dr. Robert L Smith REM TMP M+ and proven by Dr. A. Y. Tang. ANS REM M! This is the most efficient method THEN to carry out multiply precision THEN division. TMP @ REM @ - Extend remainder so that the MS REM SWAP MEXTEND ; components are cleared.
: MADVANCE ( n aa -- ) Shift aa left 16 bits and add n as the LS component in aa. Replenish the remainder from the dividend. DUP @ ?DUP If aa is not zero, IF >R shift it left 16 bits by DUP 2 + coping components to the left. OVER 4 + R> 2* CMOVE> THEN 1 OVER +! Append n to LS component in aa. 2 + ! ;
: MLEFTJUSTIFY ( aa -- ) Find bits needed to left justify the MS component in aa. The bits is stored in #SHIFT. 0 #SHIFT ! Init #SHIFT. DUP @ CELLS@ Get MS component of aa. BEGIN DUP 0< 0= Shift it until negative. WHILE 2* 1 #SHIFT +! Increment #SHIFT if not yet justified. REPEAT DROP ;
( m/ , 19may89cht) Multiple precision division.
: M/ ( aa bb -- ) Uses TMP, AA2, ANS. Quotient is in QUOT and remainder in REM. aa and bb are preserved. OVER MABS Force aa and bb to be positive. DUP MABS DUP M0= Abort if bb=0. IF ABORT" div by 0" THEN TMP M! Save bb in TMP, the divisor. TMP DUP MLEFTJUSTIFY Find bits to left justify divisor. OVER #SHIFT @ MSHL Left shift dividend. DUP #SHIFT @ MSHL Left justify divisor. @ 2* #OFFSET ! Get offset to MS component in divisor. |