This file is stored on FLT.COM
This is the comment file used by SVFLT to save all floating point routines.
STORED ON: FLT.VLI
AUTHOR: PAUL MENNEN
DATE: 29-Nov-83
The programs in this file contain routines to do floating point computations with any desired precision.
When FLT.VLI is input, the FLTDEF function is executed which causes the following message to be asked:
Enter decimal precision required (digits):
The response to this question only affects the operation of DECBIN and BINDEC and any function which
calls them. All other functions still can be carried out to any precision desired.
A floating point number is stored as a list of n elements.
The first element is a 2 element list composed of the exponent and the sign in that order.
The remaining n-1 elements of the list contain the fraction bits stored as an unsigned number 16 bits per element
with the least significant word (lsw) first and the most significant word (msw) last.
The exponent is a signed two's compliment number signifying a power of 2 assuming the binary point
is to the left of the most significant fraction bit.
The sign element is a 16 bit word whose value is either a 0 or a 1 where 0 represents a positive
number and 1 represents a negative number.
Floating point zero is treated as a special case and is not even represented by a list.
Zero is represented simply as the single atom 0.
***** functions defined in the utility file - FLT.VLI *****
(EXPO FPN) 1EXPR extracts the exponent from a Floating Point Number (FPN)
(SIGN FPN) 1EXPR extracts the sign from a floating point number
(FRACTION FPN) 1EXPR extracts the fraction from a floating point number
(MAKEFLT EX SGN FRAC) 3EXPR this is the floating point constructor function.
returns the floating point number with exponent EX, sign SGN, and fraction FRAC
(NORM FPN) 1EXPR normalize the floating point number in FPN (FPN is destroyed)
(ADD FPN1 FPN2) 2EXPR return the sum of the 2 normalized floating point numbers FPN1 and FPN2
where FPN1 and FPN2 must have the same length fraction
(ADDD FPN1 FPN2) 2EXPR returns the sum of 2 normalized floating point numbers FPN1 and FPN2
where FPN1 and FPN2 may have different length fractions
(MUL FPN1 FPN2 Z) 3EXPR return the product of the floating point numbers in FPN1 and FPN2.
if Z is not given, the length of the resulting fraction will be rounded back so that it is
equal to the length of the fraction of FPN1. if Z is given, then 16*Z bits will be removed
from the least significant end of the fraction before it is rounded.
(e.g. if Z=0 the full precision is returned)
(NEG FPN) 1EXPR return the floating point number -FPN
(SUB FPN1 FPN2) 2EXPR return the floating point number FPN1-FPN2
where FPN1 and FPN2 must have the same precision
(SUBB FPN1 FPN2) 2EXPR return the floating point number FPN1-FPN2
where FPN1 and FPN2 may have different precisions
(DIV FPN1 FPN2) 2EXPR return the floating point number FPN1/FPN2 (different precisions are allowed)
(RECIP FPN) 1EXPR return the floating point number 1/FPN
(FLOAT N LENGTH) 2EXPR convert integer N into a floating point number whose fraction has length LENGTH
Note: if you are using VLISP.C (or other FIS version) then the standard 1SUBR float function
will be overwritten when FLT.VLI is loaded.
(FIX1 FPN) 1EXPR convert the floating point number FPN into a single precision integer (16 bits)
Note: if you are using VLISP.C (or other FIS version) then the standard 1SUBR fix function
will be overwritten when FLT.VLI is loaded.
(BINDEC FPN) 1EXPR converts the floating point number FPN into an ascii string that is a decimal representation of FPN
(DECBIN DCS Z) 2EXPR converts the Decimal Character String (DCS) into a binary floating point number.
if Z is given the length of the fraction of the resulting FPN will be 16*Z bits, if not given,
the length will be chosen baised on the number of digits in the DCS. DSC may be in almost any
reasonable form that you might think of. examples of allowed forms for DCS are:
"f123521876" "f.00135" "f+123." "f-654.1234"
"f2,333,444.2E-92" "f5E7" "f1004 E 31" "f5.E7"
"f-.3E-1" "f0"
note: the leading "f" is optional if you are using VLISP.A or VLISP.B (or other non-FIS version).
If VLISP.C is used, the leading "f" is required so that VLISP can distinguish between these
variable precision floating point numbers and the standard FIS 32 bit floating point numbers.
(SUM DCS1 DCS2 Z) 3EXPR returns a decimal representation of the sum of the two floating point numbers whose
decimal representations are DCS1 and DCS2. if Z is specified, the arithmetic will
be carried out using FPN's with fractions of length 16*Z bits. if Z is unspecified the lengths
of the FPN's will be chosen to be appropriate for the number of digits specifed in DCS1 and DCS2.
(DIF DCS1 DCS2 Z) 3EXPR returns a decimal representation of the difference of the two floating point numbers whose
decimal representations are DCS1 and DCS2.
(PROD DCS1 DCS2 Z) 3EXPR returns a decimal representation of the product of the two floating point numbers whose
decimal representations are DCS1 and DCS2.
(QUOT DCS1 DCS2 Z) 3EXPR returns a decimal representation of the quotient of the two floating point numbers whose
decimal representations are DCS1 and DCS2
(FLTDEF) 0EXPR allows the user to change the decimal precision. the user is prompted for the required input.
since this function must be executed before any of the other functions in FLT.VLI, FLTDEF is
executed automatically when the FLT.VLI file is loaded.
FLTDEF
defines selector functions and sets decimal precision
.SYNONYM EXPO extracts the exponent from a floating point number
.SYNONYM SIGN extracts the sign from a floating point number
.SYNONYM FRACTION extracts the fraction from a floating point number
.NDIGITS accept decimal precision from user
.E+1 define E+1 = FLOAT(10)
.E+8 define E+8 = FLOAT(100000000)
.EW In binary EW = 1100110011001100, which is repeating pattern
found in the fraction of the floating point number .1
.NFRAC NFRAC = 2 + (NDIGITS+1) / 16 * LOG(2)
.E-1W E-1W = .1 with sufficient precision to yield NDIGITS
of precision with DECBIN and BINDEC
.LAMBDA The lambda function repeats the EW pattern N times
.E-8W .E-8W = .01
.E-8W .E-8W = .0001
.E-8W .E-8W = .00000001 (used by DECBIN and BINDEC)
.RPLACB remove the 16 L.S. bits from the fraction
.RPLACB remove the 16 L.S. bits from the fraction
MAKEFLT
this is the floating point constructor function.
returns the floating point number with exponent EX, sign SGN, and fraction FRAC
REMOVL
return the list X with its last element removed. X is destroyed as a result
SHDOWN
downshift the fraction list A by N bits and OR A1 into the msw, A is destroyed
.A its easier to do this with a reversed order fraction
SHDN
This does the work for SHDOWN. A is a fraction in reverse order (i.e. msw first)
which is downshifted by -N bits. The integer A1 is or'ed into the msw after the shifting.
.ZEROP if we are shifting by 0 bits, return A
.-16 if more than 15 shifts are required, we can save time by
shifting by a whole word
.N16 N16 is the upshift count for aligning the portion of the word to be
sent down to the next lower significant word
.LAMBDA repeat this function for each element of A. first
downshift this element and or in the bits (A1) from the previous word. The result
.A of the LOGOR function is assembled in a list
by Mapcar and is returned as the result of SHDN
SHUP
upshift the fraction A by N bits
.ZEROP if shifting by 0 bits, return A unmodified
.16 if more than 15 shifts are required, we can save time by
shifting by a whole word
.N16 N16 is the downshift count for aligning the portion of the word to be
sent up to the next higher significant word
.LAMBDA repeat this function for each element of A. first
upshift this element and or in the bits (A1) from the previous word. The result of the
.A LOGOR function is assembled in a list
by Mapcar and is returned as the result of SHUP
NORM
normalize the floating point number in A (A is destroyed)
.ZEROP is it 0 (zero has special representation)
.0 yes, you cant normalize zero, just return it.
.NSHIF NSHIF is set to the number of leading zeroes in the fraction
.ZEROP is NSHIF zero?
.A if yes, return original number
.MAKEFLT if not, make a new floating point number using the original exponent reduced by NSHIF,
the original sign, and a downshifted version of the original fraction
LEAD0
returns the number of leading zeroes in fraction list A
.AL set AL to the number of leading 0's in the msw of A
.AL if all 16 bits were zero, return 16 plus the number of leading 0's in
the fraction with the msw removed. if AL was < 16 then AL is the answer
LEAD
returns the number of leading 0's in integer X by use of a binary search
.LT detect the case where the answer is 0 first, since it is the most likely
.1 is the answer a 1,2,3, or 4?
.5 is the answer a 5,6,7, or 8?
.9 is the answer a 9,10,11, or 12?
.13 is the answer a 13,14,15, or 16?
ADD
return the sum of the 2 normalized floating point numbers X and Y
where X and Y must have the same length fraction
.EQ if X is zero, X+Y = Y
.EQ if Y is zero, X+Y = X
.DEXPO subtract the exponents
.ZEROP do X and Y have the same exponent?
.IF yes, come here and
check for a zero result (i.e. X= -Y)
.EQ the result will be 0 if the signs are opposite
.EQUAL and the fractions are identical
.0 if so, return zero
.ADDA add 2 numbers with the same exponent (but non-zero result)
.GT is the exponent of X larger than the exponent of Y?
.SHDOWN if yes, we must downshift Y
.SHDOWN if no, we must downshift X
ADDA
add 2 fractions of the same precision and common exponent (EX).
The 2 fractions are F1 and F2 and their signs are S1 and S2 respectively
(i.e. sign magnitude notation)
.FR1 FR1 will be fraction 1 in twos compliment
.FR2 FR2 will be fraction 2 in twos compliment
.S12 set S12 to the sum of the two sign words plus the last carry out.
thus S12 is the most significant 2 bits of the twos compliment result
.LOGXOR do the two fractions have the same sign?
.COND if yes, come here.
do this if the result was positive with no overflow
.EQ do this if the result was positive with an overflow
.EQ do this if the result was negative with an overflow
.T do this if the result was negative with no overflow
.IF do this if the fraction had different signs (no overflows are possible here).
is the result positive?
.0 yes
.1 no, its negative
NEGFR
returns the twos compliment of fraction F by complimenting all bits and adding a 1
.CRY here is the one which is to be added in
.CRY compliment all bits, add in the 1, and propagate the carry
ADDD
returns the sum of 2 normalized floating point numbers X and Y
where X and Y may have different length fractions
.EQ if X is 0, X+Y = Y
.EQ if Y is 0, X+Y = X
.LDIF LDIF is the number of words longer that fraction X is than fraction Y
.ZEROP if X and Y are the same length we can add them with ADD
.GT come here if fraction X is the longest
.FR fraction Y will be extended into FR until it is the same length as fraction X
.0 append enough zeroes to the beginning of fraction Y
.ADD now we can use ADD
.T come here if fraction Y is the longest
.FR fraction X will be extended into FR until it is the same length as fraction Y
.0 append enough zeroes to the beginning of fraction X
.ADD now we can use ADD
MUL
return the product of the floating point numbers in X and Y.
if Z is not given, the length of the resulting fraction will be rounded back so that it is
equal to the length of the fraction of X. if Z is given, then Z elements will be removed
from the least significant end of the fraction before it is rounded.
(e.g. if Z=0 the full precision is returned)
.OR are either multiplier or multiplicand zero?
.0 if yes, the answer is zero
.EX to multiply, we add exponents and exclusive-or the signs
.OR now check for overflow or underflow
.GE we underflowed if both exponents were negative but their sum was positive
.LT we overflowed if both exponents were positive but their sum was negative
.1 L is a list used to collect partial products
.LAMBDA do this function for each element of fraction Y.
we now have the D1'th element of fraction Y
stored in D
.LAMBDA do this function for each element of fraction X
we now have the (E1-D1+1) th element of fraction X stored in E
.MULDE multiply the 2 unsigned numbers D and E
.ATTACH attach LOW to the list of partial products
.ATTACH attach HIGH to the list of partial products one position higher
.CRY clear the carry to initialize
.MXY now add up the partial products into list MXY
.LAMBDA for each list inside L, add up all the elements along with the carry from the
.L previous add. save the new carry in CRY
.ROUND now round the answer MXY back to desired length
.NORM and return the answer normalized
TXT1
this function merely provides some text which helps explain MUL
In case you could not see how MUL works, this examle may help.
Suppose (FRACTION X) = (f e d) and (FRACTION Y) = (c b a)
where a,b,c,d,e,f are all sixteen bit integers.
I will use the notation cdL and cdH to denote the Low and High 16 bits of the 32 bit product (c times d)
The multiplication is done in parts similar to hand decimal multiplation:
d e f
x a b c
___________________________________
cfH cfL
ceH ceL
cdH cdL
bfH bfL
beH beL
bdH bdL
afH afL
aeH aeL
adH adL
____________________________________________
xy6 xy5 xy4 xy3 xy2 xy1
the columns are added (and carrys propagated) to generate the six word result
the partial product list (L) for this problem would look like this:
L = ((0 cfL)
(cfH ceL bfL)
(ceH cdL bfH beL afL)
(cdH beH bdL afH aeL)
(bdH aeH adL)
(adH))
then the columns are added up to produce the list MXY:
MXY = (xy1 xy2 xy3 xy4 xy5 xy6)
ROUND
the length of MXY is now the sum of (length (fraction x)) and (length (fraction y))
If Z is given, this routine shortens this fraction (MXY) by removing Z of its least sig. elements.
If Z is not given, then MXY is shortened so that it is the same precision as (fraction X)
ROUND then increments the shortened result by one if necessary to implement rounding
.RND do this in case Z is zero
.Z if Z is not given, set Z to the length of fraction Y
.UNTIL strip the unwanted words from MXY, when done RND will contain the last word deleted
.LT is the most sig. bit of RND a one?
.SETQ if yes, do this, otherwize we are done
.MXY add a one to MXY and then propagate the carrys thru
ATTACH
add E to the Nth list in L as its new last element (L is a list of lists).
if L does not have N lists in it, then add a new list at the end of L
which contains the single element E.
NEG
return the floating point number -X
SUB
return the floating point number X-Y
where X and Y must have the same precision
SUBB
return the floating point number X-Y
where X and Y may have different precisions
DIV
return the floating point number X/Y (different precisions are allowed)
RECIP
return the floating point number 1/X
a successive approximation technique is used based on the fact that the
solution to the equation: R = 2R - X*R*R is: R = 1/X
.ZEROP is X zero?
.PRINT cant divide by 0, print message
.ROLD set ROLD to the first guess
.RNEW set RNEW to the result of the first iteration
.RCOUNT this iteration counter is only for
measurement use and may be removed
.UNTIL repeat this loop until the reciprical has converged to 8 bits beyond the precision of X.
.EQUAL are all fraction bits the same, not counting the least significant word?
.LE do the least significant words agree in the lower 8 bits?
.SETQ if one of the above conditions failed, continue the loop,
move RNEW into ROLD
.RCOUNT increment the iteration counter
.RNEW compute the next iteration
.RNEW return RNEW as the answer
GUESS
returns the value used as the initial guess of the reciprical of X.
A fraction is formed with the same sign and fraction as X but with a complimented exponent.
This is a lousy guess. We should be able to do better!
FLOAT
convert integer N into a floating point number whose fraction has length LENGTH
.0 zero is represented as 0
.GE is it a positive number?
.NORM yes
.MAKEFLT build the floating pt. number to be normalized
.16 exponent = 16
.0 sign = positive
.PROGN add zeros to the front of N to get the fraction
.NORM no, its a negative number
.MAKEFLT build the floating pt. number to be normalized
.16 exponent = 16
.1 sign = negitive
.PROGN add zeros to the front of -N to get the fraction
FIX
convert the floating point number X into a single precision integer (16 bits)
.0 if it is zero, return 0
.SETQ compute shift count
.GT will it fit in an integer?
.LT is it less than 1 in magnitude?
.CAR if its positive just downshift the msw of fraction by EX16 positions
.CAR if its negative. do the same thing, but return the negation.
BINDEC
return a decimal representation of the floating point number S as a character string having NDIG
digits. if NDIG is not specified, the number of decimal digits will be chosen to be
appropriate for the precision in the input variable S. In either case, the number of decimal
digits computed will never be greater than NDIGITS as specified by the FLTDEF function. The
first character in the decimal representation is "f" so vlisp wont confuse it with its
internal floating numbers.
.ZEROP Is S zero?
.0 Yes, simply return zero
.SLEN SLEN is the length of the fraction of S
.SGN1 SGN1 saves the sign of S for later use
.-1 initialize EXP10, it will be exponent of answer
.DONE we are not done yet
.DEC initialize char string to contain answer
.NDIG if NDIG was unspecified choose NDIG=(16*LOG(2)*SLEN)-3
.NDIG dont allow NDIG to be larger than NDIGITS
.GETSF get the fraction bits from S into SF to be converted to decimal
.UNTIL do this once for each decimal digit desired
.ADDFR multiply SF by 10 (SF = 2*SF + 8*SF)
.NEWL the integer part of result is the next decimal digit, save in DEC
.RPLACA zero out the integer part of SF and go loop
.ADDFR compute one more digit for rounding purposes
.GE is the last digit 5 or greater?
.SETQ if yes, come here to round the result up
.ADD1 increment the last digit computed
.WHILE if that digit just changed to a 10, change it again to a zero and
go to the next digit and
.0 increment it
.NEXTL continue this process until a digit other than
.ADD1 9 is incremented or the end of the string is reached
.SETQ add the E to indicate power of 10 and reverse string
.COND get rid of the first digit if it a zero
.CONS insert the decimal point between the first and second digits
.SGN1 put a minus sign in front if S was negative
.NEWL use f as 1st character, so implode wont internalize result
as a fis floating pt number
.NCONC attach exponent and implode the result to give the desired string
GETSF
this routine extracts the binary bits from S (into SF) that we want
to convert to decimal. S is repetativly multiplied by 10, 1E8, .1, or 1E-8
as appropriate to get S withen the desired range (i.e. an exponent between 0 and -3)
.1 is exponent > or = to 1?
.28 yes, is exponent > or = to 28?
.MUL yes, multiply by 1E-8 and adjust EXP10 to keep track
.MUL no, multiply by .1
.INCR and adjust EXP10 to keep track
.-4 no, is exponent < or = -4?
.-30 yes, is exponent < or = to -30 ?
.MUL yes, multiply S by 1E8 and adjust EXP10 to keep track
.MUL no, multiply S by 10
.DECR and adjust EXP10 to keep track
.T exponent is in desired range, set DONE
.SF shift the fraction of S right 0,1,2, or 3 positions (given by exponent of S)
then attach an extra word of zero's to the most signifiacant end.
put fraction into SF (the extra word is the integer part of SF)
ADJUST
used by DECBIN to adjust the result to account for the decimal exponent
.GT is the exponent positive?;
if yes, multiply SF by 10 EXP10 times. to
.WHILE make this faster, do it in groups of 8 for as long as possible,
(i.e. multiply by 1.0E8) and
.UNTIL then finish up by multiplying by 10 the correct number of times
.LT is the exponent negative?
.WHILE If yes multiply by .1 -EXP10 times using groups
of 8 as above
DECBIN
return a floating point number with N words of precision in its fraction that
corresponds to the decimal representation of the character string DEC. If N is not
specified, it will be chosen to be appropriate baised on the number of digits in DEC.
In either case, the number of fraction words computed is never larger than NFRAC as
specified by the FLTDEF function
.IF is DEC zero (a special case)
.0 yes, return zero
.0 initialize multiple length
integer SF, change DEC from character string to a list
.WHILE do this loop once for each input character
.DEC1 DEC1 is the next input character
.NUMBP is the next character a decimal digit? If yes, do the following.
.PERIOD adjust EXP10 if we are after the decimal point
.INCR NUMDIG counts number of digits encountered
.0 add a 0 to the most sig. end of SF
.ADDFR SF = (SF*2 + SF*8) = SF * 10
.SF Now add the decimal digit to the multiple length integer
.UNTIL and propagate the carry's thru
.REMOVL remove the most sig word form SF if its 0, thus SF grows
only as fast as needed
.69 The character isnt a digit, is it the letter E?;
yes, interpret all characters following the E
.SETQ as a decimal exponent. set DEC to nil indicating that we
are done scanning the input
.PERIOD If the character is a decimal point, remember
the fact by PERIOD = True
.EQ If the character is a minus sign, set end result to negative
.IFN if N, the desired length of fraction is unspecified,
choose it on the basis of NUMDIG
.NFRAC dont allow N to be larger than NFRAC
.EXSF EXSF = # of bits in the integer SF
.LT if SF is shorter than desired fraction, append zeroes
at least sig end until proper length
.GT if SF is longer than desired, delete least sig words
until SF is the proper length
.SF turn the multiple precision integer SF into a floating
point number
.ADJUST adjust result to take into account the decimal exponent
.SF return SF as the result
SVFLT
save all the above routines in the file FLT.VLL