housel@en.ecn.purdue.edu (Peter S. Housel) (07/29/89)
#! /bin/sh
# This is a shell archive. Remove anything before this line, then feed it
# into a shell via "sh file" or similar. To overwrite existing files,
# type "sh file -c".
# The tool that generated this appeared in the comp.sources.unix newsgroup;
# send mail to comp-sources-unix@uunet.uu.net if you want that tool.
# If this archive is complete, you will see the following message at the end:
# "End of archive 1 (of 3)."
# Contents: README MANIFEST NOTES cfu.x mlf8.x vfprintf.c
# Wrapped by housel@en.ecn.purdue.edu on Fri Jul 28 13:35:03 1989
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'README' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(1463 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
XThis package contains a floating-point arithmetic package for Minix-PC
Xversion 1.3, and a few ``mostly ANSI-compatible'' additions to the C
Xlibrary. Installation instructions are in the file INSTALL, and notes
Xon implementation are in the file NOTES.
X
XThis is more-or-less a medium quality implementation of floating point
Xarithmetic for the Minix C compiler. No 80x87 math coprocessor is
Xrequired (or even useful) to use this package. The binary format
Xshould be identical to that of most IEEE 754 implementations, and many
Xoperations give bit-for-bit identical results. It is by no means,
Xhowever, fully standard-conforming. (There is no support for
Xinfinities, NaN's, or denormalized numbers, for example.)
X
XWhenever compiling a program which uses floating-point constants,
Xand/or when linking any program which uses floating-point arithmetic,
Xthe "-f" option of the C compiler should be used. (This causes the
Xfloating point constant postprocessor to be run, and adds the floating
Xpoint library to the list of libraries linked with the program.)
X
XAlthough the goal is to eventually provide all of the math functions
Xmentioned in the math.h file (and in section B4 of K&R2), not all of
Xthese are available yet. (Due to portability problems, the BSD 4.3
Xfreely-redistributable portions of libm are not yet usable.)
X
XAny bug reports, suggestions, questions, or improved code would be
Xappreciated.
X
X-Peter S. Housel- housel@ecn.purdue.edu ...!pur-ee!housel
END_OF_FILE
if test 1463 -ne `wc -c <'README'`; then
echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'MANIFEST' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'MANIFEST'\"
else
echo shar: Extracting \"'MANIFEST'\" \(1756 characters\)
sed "s/^X//" >'MANIFEST' <<'END_OF_FILE'
X File Name Archive # Description
X----------------------------------------------------------
X README 1
X MANIFEST 1 This shipping list
X CRCLIST 3
X INSTALL 3
X Makefile 2
X NOTES 1
X PATCHDATES 3
X _mult.c 3
X _poly.c 3
X addsub.x 2
X atof.c 3
X cc.c.cdif 2
X ceil.c 3
X cff.x 2
X cfi.x 2
X cfu.x 1
X cif.x 2
X cmf8.x 3
X cuf.x 2
X dvf8.x 2
X exp.c 2
X fabs.c 3
X fat.x 3
X float.h 2
X floor.c 3
X fpp.c 2
X fprintf.c 3
X frexp.x 3
X ldexp.x 3
X log.c 3
X log10.c 3
X math.h 3
X mlf8.x 1
X modf.x 2
X ngf8.x 3
X norm4.x 2
X norm8.x 2
X pow.c 3
X printf.c 3
X ret8.x 3
X return.x 3
X sprintf.c 3
X sqrt.c 3
X stdarg.h 3
X stddef.h 3
X strtod.c 2
X strtod_aux.x 2
X trp.x 2
X vfprintf.c 1
X vprintf.c 3
X vsprintf.c 3
X zrf4.x 3
X zrf8.x 3
END_OF_FILE
if test 1756 -ne `wc -c <'MANIFEST'`; then
echo shar: \"'MANIFEST'\" unpacked with wrong size!
fi
# end of 'MANIFEST'
fi
if test -f 'NOTES' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'NOTES'\"
else
echo shar: Extracting \"'NOTES'\" \(23050 characters\)
sed "s/^X//" >'NOTES' <<'END_OF_FILE'
XHere are some notes on the implementation of the PC-Minix floating
Xpoint routines, and a (rather extensive) list of bibliographical
Xreferences on floating point arithmetic.
X
X1. IMPLEMENTATION
X
X1.1 A Brief Overview of Floating-Point Math
X
XSeveral of the references have better explanations of what floating
Xpoint numbers are and how arithmetic works on them. The Sterbenz book
Xgives a fairly rigorous mathematical treatment; the first of the
XPlaugher references provides an excellent in-depth explanation from the
Xuser's point of view. Briefly, my attempt:
X
XMost of the computations done by computers are done on integers. These
Xare simply represented - each bit of a binary quantity is used to
Xrepresent a power of two. The least significant bit represents the
Xnumber 1; the next, the number 2; and so on.
X
XFloating-point numbers are often referred to as "reals". This is a
Xmisnomer; the number of distinct values that a 32-bit floating-point
Xnumber can take on is the same as, or even less than, the number of
Xdistinct 32-bit integers. However, floats can represent a much larger
Xrange of values. As an example, 32-bit binary integers can range from
X-2^31 to 2^31. A common format for 32-bit floats can have magnitudes from
X2^-127 to 2^128.
X
XThis is accomplished by splitting up the bits of the binary quantity
Xinto three fields: a sign, a biased exponent, and a significand
X(commonly called the mantissa). The significand represents a number
Xbetween zero and one, and is usually normalized so that it ranges from
X1/2 to one. The number represented is then:
X (sign) mantissa * 2^exponent
X
XThe exponent is stored in the bit field with a bias added. The bias is
Xgreater than or equal to the minimum actual exponent, so that the
Xbiased exponent is always an unsigned quantity, and no exponent sign
Xbit is required.
X
XZero is usually represented by a zero biased exponent and a zero
Xmantissa. This allows a floating point zero to be represented by an
X"all-zeroes" binary number.
X
X1.2 IEEE Standard Number Representation
X
XThe IEEE (Institute of Electrical and Electronics Engineers) Standards
XCommittee has defined a standard behaviour for floating point
Xarithmetic implementations. This standard, IEEE 754-1985, defines the
Xformats, operations, conversions, and exceptions supported by a
Xconforming implementation. Four formats are defined: single precision,
Xsingle extended, double precision, and double extended. Field widths
Xand interpretations are defined for the non-extended formats as follows:
X
X Single Double
XMantissa bits 24 53
XExponent max +127 +1023
XExponent min -126 -1022
XExponent width +127 +1023
XExponent width in bits 8 11
XFormat width in bits 32 64
X
X(The IEEE standard numbers the bias and exponent differently from a lot
Xof other floating point literature, including the ANSI float.h and the
Xsource code of this package.)
X
XIn addition, representations are defined for positive and negative
Xinfinity. There are also quantities called NaN's, for Not a Number.
XThese result from illegal operations, and can either cause an exception
Xor propagate as operations are performed on them.
X
XThe standard defines several arithmetic operations: add, subtract,
Xmultiply, divide, remainder, square root, and format conversion. The
Xbehaviour of each these operations when given a number, infinity, or
XNaN is detailed. Rounding modes and exceptions are also defined.
XNumbers can gracefully underflow as "denormalized" numbers with a zero
Xbiased-exponent.
X
XMany of these features are nontrivial or slow to implement in software.
XFor this reason, this implementation does not strictly conform to the
Xstandard. The number of bits for single and double precision, the
Xbitfield layout, and the exponent bias are the same. However,
Xinfinities, NaN's, and denormalized numbers are not recognized. There
Xis only one rounding mode, round-to-nearest with round-to-even in the
Xtie case.
X
X(A newer standard, IEEE 854-1987, is similar to 754-1985 but is more
Xgeneral. For example, numbers may be represented either in binary or
Xdecimal, and the formats are not specified exactly.)
X
X1.3 Floating-Point Math in C
X
XThe C language defines three floating point sizes: float, double, and
Xlong double. (long double is new with the draft ANSI C standard, and
Xis not available in the Minix C compilers.) The standard states that
Xprecision(float) <= precision(double) <= precision(long double), but
Xthat they may have the same representation.
X
XThe primary floating point datatype in C is double. In K&R C
X(including the Minix compiler), all floating point quantities are
Xextended to double before doing arithmetic on them. This is part of
Xthe "usual arithmetic conversions." (ANSI Standard C will relax this,
Xso that expressions containing strictly floats may be done using
Xfloat-arithmetic, provided the result is the same "as if" doubles had
Xbeen used.) Floating point arguments to functions are also widened to
Xdouble, except in Standard C when an in-scope prototype says that the
Xfunction takes a float argument. All of the standard library functions
Xreturn doubles.
X
X1.4 The ACK routines
X
XThe Minix C (and Pascal) compilers are based on the Amsterdam Compiler
XKit, or ACK. The compiler generates 8086/8088 assembly code files from
XC program text using several passes (each found in /lib or /usr/lib) as
Xfollows:
X
X file.c <file.i> file.k file.m
X--------> cpp --------> cem --------> opt --------> cg --------> file.s
X
X(The "-v" option of the "cc" command allows the invocation of each of
Xthese passes to be seen.)
X
XThe "cpp" program is a version of the usual C macro preprocessor, which
Xdoes "#-command" processing and macro substitution on the C source
Xcode.
X
XThe ACK compiler "frontends" take as input a program in the source
Xlanguage (in the case of "cem", C) and generate EM machine code. The
XEM ("Encoding Machine") machine (which does not exist in hardware, but
Xwas defined as a generic machine for the ACK system to generate code
Xfor) is a byte-addressable, stack-oriented processor. It supports
Xoperations on integers (signed and unsigned), pointers, sets, and
Xfloating-point numbers. Each data type is available in one or more
Xsizes, which are identified by their lengths in bytes. The
Xarchitecture was designed based on research on instructions used by
Xblock-structured high-level languages.
X
XThe output of "cem" is in EM compact assembly format. This is EM
Xassembly-language represented in a compact byte-code representation.
XThe definition of this format, as well as the format of the ASCII
Xassembly language and the object code format, is described in a
Xdocument included with the ACK system. An excerpt of this document was
Xposted to comp.os.minix as articles <1650@ast.cs.vu.nl> and
X<1651@ast.cs.vu.nl>.
X
XThe "opt" program takes the "cem" output and performs some
Xmachine-independent peephole optimizations on the EM code. There are
Xas many as 400 EM instruction sequences which the optimizer can
Xtransform into more efficient ones. This take the pressure of
Xgenerating optimal code off of the frontend and reduces code
Xduplication between frontends.
X
XThe ACK system also includes a global optimizer, which is capable of
Xmany more powerful optimizations. This program is, however, too big to
Xrun under Minix and is not included or available.
X
XThe code generator, "cg", reads the output of "opt" and translates the
XEM assembly into code for the target machine. For Minix-PC, this is
X8086 assembly language in the Minix compressed-assembly format. Most
XEM instructions translate into a small number of 8086 instructions.
XHowever, there are several which would translate to a larger number,
Xperhaps six or more. These are broken out into subroutines.
X
XThese auxiliary subroutines are found "near the end" of the C library,
Xand do such things as implement switch-statement jumptables, convert
Xbetween various sizes of integers, do 32-bit arithmetic, and so on.
X
XOne of the files in this set of auxiliary routines is called "fakfp.s,"
Xfor "fake floating point." The "cg" program is capable of generating
Xcalls to a set of auxiliary routines which implement floating-point
Xarithmetic, but no real routines are included. Instead, fakfp.s
Xincludes a stub-routine which generates an EM "Illegal EM instruction"
Xtrap and aborts the program.
X
XThis package includes many of the omitted routines. There are stubs in
Xfakfp.s for the following:
X
X.mlf .mlf4 .mlf8 multiply
X.dvf .dvf4 .dvf8 divide
X.ngf .ngf4 .ngf8 negate
X.adf .adf4 .adf8 add
X.sbf .sbf4 .sbf8 subtract
X.cmf .cmf4 .cmf8 compare
X.zrf .zrf4 .zrf8 zero
X.fif .fif4 .fif8 split into exponent and fraction part
X.fef .fef4 .fef8 multiply and split integer and fraction part
X.cif convert integer to float
X.cfi convert float to integer
X.cuf convert unsigned to float
X.cfu convert float to unsigned
X.cff convert float to float (float <==> double)
X
XBecause of the "usual conversions" rules for expressions in pre-ANSI C,
Xall arithmetic is done in "double" floating point, even when all of the
Xelements are of type "float". This makes it unnecessary to include
Xmost of the auxiliary routines for doing arithmetic on float-sized
Xquantities. (The exceptions to this are the conversion routines, which
Xconvert between all different sizes of floating point and integer numbers.)
XThe "C" compiler does not generate FEF or FIF instructions either, so
Xthese are also not implemented.
X
X1.5 Implementation notes on the ACK routines
X
XAll of the routines named above expect their arguments on the stack,
Xand return their results there. This, and the use of the normalization
Xsubroutines, determine most of the structure of the routines.
X
XMost of the floating point ACK-auxiliary routines call .norm8 or .norm4,
Xwhich take a denormalized mantissa, an exponent, and a sign, and pack
Xthem into a standard-format normalized number. (Note that these are
Xnot actually EM instructions.)
X
XWhen calling .norm8 (.norm4), there should be a 64-bit (32-bit)
Xmantissa on the stack. The radix point is assumed to be in front of
Xbit 52 (bit 23). The bx register should contain a properly biased
Xexponent. The high bit of dh is 1 for negative numbers and 0 for
Xpositive numbers. The dl register should contain eight bits of information
Xto be used for rounding.
X
XMost of the work of normalization involves shifting the mantissa up or
Xdown until there is a 1-bit in the position where the implied 1-bit
Xwould be, and packing the exponent and sign into the number. When bits
Xare lost (by shifting off of the lower end), it is necessary to round
Xthe number to minimize the error introduced. (Simply truncating those
Xbits which are lost, known as "chopping", has a high average error and
Xa high error bias.) Round-to-nearest is used, and tie cases round to
Xeven (lowest bit is jammed to zero). When 1-bits are shifted out of the
Xeight-bit rounding register, its least significant bit "sticks" to 1,
Xto allow those bits to participate in the rounding process. (Only three
Xbits, usually one each of "guard," "round," and "sticky," are necessary
Xfor this type of rounding. However, coding for a three-bit register
Xwould be larger and slower when there are already instructions for
Xmanipulating eight-bit registers.) Overflows cause EM traps, and
Xunderflows return a zero.
X
XThe .zrf4 and .zrf8 routines simply push a zero of the appropriate
Xsize. The .ngf8 routine flips the sign bit of its argument.
X
XMost of the rest of the routines have a couple of code segments in
Xcommon which set the si and/or di registers to point to the
Xarguments or result, and which unpack a number into its
Xsign/exponent/mantissa components.
X
XThe addition and subtraction routines, .adf8 and .sbf8, are the most
Xdifficult. The algorithm used is described in pseudocode in Rick
XGrehan's "Floating Point without a Coprocessor" (see the references).
XBasically the scheme is to line up the radix points of the two numbers
Xand add or subtract them depending on their relative magnitudes. Most
Xof the problems are caused by the scarcity and small size of the
Xregisters in the 8086.
X
XMultiplication (mlf8.x) consists, basically, of adding the exponents
X(being careful with the bias), exclusive-or-ing the signs, and
Xmultiplying the mantissas. In order to take advantage of the 8086 MUL
Xinstruction, Algorithm M from section 4.3.1 of Knuth's Seminumerical
XAlgorithms (again, see the references) is used. The loop is unrolled
Xand some special cases are commented out or otherwise avoided. (No
Xdoubt a faster multiply would have been possible using the algorithms
Xin section 4.3.3, but that would have required much more study than I
Xwas willing to do, as well as requiring anyone else who maintained the
Xcode to understand it. The method used is reasonably fast and
Xreasonably straightforward.) A 106-bit mantissa is generated; the most
Xsignificant of these are used for the answer and the rest are packed
Xinto a rounding byte.
X
XDivision (dvf8.x) is like multiplication except that the exponents are
Xsubtracted and the mantissas are divided. Again, see the Grehan
Xarticle for the pseudocode. Division by zero causes an EM trap.
X
XThe convert-float-to-float instructions (cff.x) convert either an
Xeight-byte float into a four-byte float, or vice versa. This basically
Xinvolves unpacking the number, realigning the mantissas, and adjusting
Xthe bias. The normalization routines are used to re-pack the result.
X
XConverting a floating point number into an integer (cfi.x) or an
Xunsigned integer (cfu.x) involves realigning the mantissa until
Xthe radix point is at a known point, and taking everything before
Xthat point as the integer.
X
XThe inverse operations, in cif.x and cfi.x, are implemented by placing
Xthe integers into a mantissa field, placing the radix point at the
Xend, and normalizing.
X
XThe floating-point comparison routine, in cmf8.x, uses a quick
Xalgorithm outlined by Plaugher in the "Safe Math" article. It takes
Xadvantage of the way floating point numbers are packed into memory, in
Xthat larger floating point numbers, when taken as integers, compare as
Xgreater than smaller ones. The main point of difficulty is that
Xfloating point numbers use sign-magnitude representation, so they must
Xbe converted to pseudo-2's-complement before the comparison. The
Xreturn value is a 16-bit integer on the stack. If it is less than
Xzero, the first argument was less than the second; if zero, they
Xcompared equal; if greater than zero, the first was greater than the
Xsecond.
X
XWhen an EM instruction is supposed to indicate failure, the auxiliary
Xroutine will usually call the EM trap routines. A "trp.x" file is
Xalready included with the library, but the one presented here provides
Xmessages for more trap codes. All of the codes are defined in the EM
Xdefinition documents.
X
X(In general, in all of these routines, watch out for the bits whizzing
Xthis way and that. Be especially careful of the ones that fall off of
Xthe ends.)
X
X1.6 Assorted other library routines and their implementation
X
XThe above routines are all called implicitly by C code. Most C
Xlibraries also provide several explicitly-callable subroutines that do
Xfloating-point math.
X
XThe ldexp() and frexp() routines are implemented in assembly code.
Xfrexp() is basically a frontend for the number-unpacking routines,
Xand ldexp() is a frontend for the normalization routine's repacking
Xfunction. The ldexp() function is careful to keep the exponent within
Xthe legal range, and returns HUGE_VAL (defined in <math.h>) of the proper
Xsign when it is not.
X
XThe modf() function is also written in assembly. It unpacks an
Xarbitrary floating-point number into two floating-point point components,
Xan integer and a fractional part. Numbers which are large enough or small
Xenough to be entirely one or the other are special cases. Otherwise,
Xthe number is shifted until the radix point lines up with the end of
Xthe integer part, the result is split into two parts and each half is
Xrenormalized.
X
XThe Standard C routines used to convert human-readable forms of
Xfloating point numbers into their internal representations are atof()
Xand strtod(). They are mostly the same, except that strtod() can store
Xa pointer to the end of the number being converted. Therefore, atof()
Xis defined in terms of strtod().
X
XPart of strtod() is in C, although the entire subroutine was translated
Xfrom VAX/11 assembly (the freely-redistributable 4.3BSD-Tahoe version
Xof atof(), which is not included here). The C portion uses _mul10add()
X(written in assembly, and included in strtod_aux.x) to accumulate
Xdigits in a 64-bit accumulator until it overflows. Then the optional
Xexponent is parsed, the base-10 exponent is computed, and _adjust()
X(also in strtod_aux.x) is called to multiply the 64-bit integer by
X10^decexp, and pack the result into a double-format floating point
Xnumber. All sorts of clever tricks are used in _adjust() to maximize
Xspeed and precision; credit for these goes to the author of atof.s.
X
XStandard C defines vfprintf() as a for producing formatted output from
Xa variable-length vector of arguments. The version included here is
Xa very slight modification of 4.3BSD-Tahoe's _doprnt() code (which
Xis freely redistributable.) vfprintf.c should be compiled with
Xthe "-DFLOATLIB" flag if the floating-point conversion options
Xare to be included.
X
XAll of the other members of the "printf() family" are defined in terms
Xof vfprintf(). To my knowledge, all of these (along with the other
Xmath.h functions provided) conform to the relevant portions of the
Xdraft proposed ANSI C standard.
X
XPreliminary versions of several functions traditionally provided in
X"libm" are also included. A few of them are not especially accurate,
Xperhaps to only about 40 bits. They were written in haste and
Xignorance, and will be replaced with provably better routines at the
Xearliest opportunity.
X
Xceil() and floor() truncate numbers to integral floating-point
Xnumbers in different ways. They use modf() to split up the
Xinteger and fraction parts of their arguments.
X
X_poly() and _mult() are internal routines, and are not meant to be
Xcalled by user-level programs. The names are inspired by the Plaugher
X"Safe Math" article. _poly() evaluates a polynomial function with
Xgiven coefficients at a given value. No overflow checking is done, so
Xthe value should be "reduced". _mult() multiplies two numbers and
Xreturns HUGE_VAL (from <math.h>) with errno properly set if the result
Xis out of range, or the result otherwise.
X
Xlog() returns the natural logarithm of its argument. The series expansion
Xused was found in A Survey of Numerical Mathematics, and the constants
Xare from the BSD 4.3-Tahoe libm source. log10() is trivially defined
Xusing log().
X
Xexp() uses the argument reduction described in the previous reference.
XThe constants, except for the polynomial coefficient vector, come from
Xthe libm source.
X
Xpow() uses log() and exp(), with a bit of code to handle special cases.
X
Xsqrt() was inspired quite directly by the Plaugher article "Do-it-yourself
XMath Functions". The linear-approximation constants are a least-squares
Xestimator for the square root over [.5,1).
X
X1.7 The EM-file postprocessor
X
XFor portability and machine-independence, the "cem" compiler frontend
Xoutputs floating-point constants in human-readable ASCII notation
Xexactly as they were included in the original source program. It is
Xintended that the code generator handle the machine-dependent
Xconversion to binary representation. However, the "cg" program
Xdistributed with Minix-PC is not able to do this conversion, and simply
Xtreats all float constants as 0.0. (The 1.2-1.3 generator prints a
Xwarning, "Warning: dummy floating point constant(s)".) Someone
Xwith the source to "cg" could add this capability, but this option
Xis not available to everyone (myself included).
X
XTo solve this problem, it is necessary to introduce another pass into
Xthe C compilation process. The "fpp" postprocessor takes the
XEM-assembly output of the "opt" program, and copies it to the output.
XHowever, floating point constants are intercepted and converted to
Xbinary format using atof(). The result is split up and output as
Xeither two or four 16-bit integer constants.
X
XThe result is that the "cg" program never sees floating point constant
Xarguments. The conversion and copy is very fast, but since it is not
Xalways necessary, the "fpp" pass is only run when the "-f" option is
Xgiven to the C compiler driver.
X
XThe EM format documentation mentioned above is, unfortunately, required
Xto understand the logic of the fpp.c program. (I have also written a
Xprogram to convert compact EM-assembly into a sort of human-readable
XEM-assembly (though not identical to the format used in the ACK because
XI have not seen enough to get it right). If you would like a copy of this
Xprogram please write and ask.)
X
X2. REFERENCES
X
X1. ANSI/IEEE, IEEE Standard for Binary Floating-Point Arithmetic, ANSI/IEEE
XStd. 754-1985. New York, N.Y.: The Institute of Electrical and Electronics
XEngineers, Inc., 1985.
X
X2. ANSI/IEEE, IEEE Standard for Radix-Independent Floating-Point Arithmetic,
XANSI/IEEE Std. 854-1987. New York, N.Y.: The Institute of Eletrical and
XElectronics Engineers, Inc., 1987.
X
X3. Cody, William J., and William Waite. Software Manual for the
XElementary Functions, Englewood Cliffs, NJ: Prentice-Hall, 1980.
X
X4. Grehan, Rick. "Some Assembly Required: Floating-Point Without a
XCoprocessor, Part 1," Byte 13(9): pp. 313-319, Sep. 1988.
X
X5. Grehan, Rick. "Some Assembly Required: Floating-Point Without a
XCoprocessor, Part 2," Byte 13(10): pp. 293-298, Oct. 1988.
X
X6. Hart, John F., E.W. Cheney, et al. Computer Approximations, New York:
XRobert E. Krieger Publishing Company, 1978.
X
X7. Kernighan, Brian W., Dennis M. Richie. The C Programming Language,
XSecond Edition. Englewood Cliffs, NJ: Prentice-Hall, 1988.
X
X7. Knuth, Donald E. The Art of Computer Programming, Second Ed., Vol. 2,
XChapter 4. Reading, Mass.: Addison-Wesley, 1981.
X
X8. Motorola, Inc. MC68881/MC68882 Floating-Point Coprocessor User's Manual,
XEnglewood Cliffs, NJ: Prentice-Hall, 1987.
X
X9. Ochs, Tom. "Theory and Practice," Computer Language, 6(3): pp. 67-81,
XMarch 1989.
X
X10. Plaugher, P.J. "Computer Arithmetic," Programming on Purpose, Computer
XLanguage, 5(2): pp. 17-23, Feb. 1988.
X
X11. Plaugher, P.J. "Properties of Floating-point Arithmetic,"
XProgramming on Purpose, Computer Language, 5(3): pp. 17-22, Mar. 1988.
X
X12. Plaugher, P.J. "Safe Math," Programming on Purpose, Computer
XLanguage, 5(5): pp. 17-21, May 1988.
X
X13. Plaugher, P.J. "Do-it-yourself Math Functions," Programming on Purpose,
XComputer Language, 5(6): pp. 17-22, Jun. 1988.
X
X14. Sterbenz, Pat. Floating Point Computation, Englewood Cliffs, NJ:
XPrentice-Hall, 1972.
X
X15. Tanenbaum, Andrew, et al. "A Practical Tool Kit for Making Portable
XCompilers," Communications of the ACM 26(9): pp. 654-660, September 1983.
X
X16. Wilson, Pete. "Floating-Point Survival Kit," Byte 13(3): pp. 217-226,
XMar. 1988.
X
X17. Young, David M., Robert Todd Gregory. A Survey of Mumerical Mathematics,
XVolume I, Reading, Mass: Addison-Wesley, 1972.
END_OF_FILE
if test 23050 -ne `wc -c <'NOTES'`; then
echo shar: \"'NOTES'\" unpacked with wrong size!
fi
# end of 'NOTES'
fi
if test -f 'cfu.x' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'cfu.x'\"
else
echo shar: Extracting \"'cfu.x'\" \(3044 characters\)
sed "s/^X//" >'cfu.x' <<'END_OF_FILE'
X.define .cfu
X|
X| floating point => unsigned conversion routines
X| author: Peter S. Housel 6/12/89,6/14/89
X|
X
XBIAS4 = 0x7F - 1
XBIAS8 = 0x3FF - 1
X
X .text
X .globl .cfu
X|
X| on entry dx=source size, cx=dest. size
X|
X.cfu:
X push si | save register variable
X mov si,sp
X
X cmp dx,#4
X jne cfu8
Xcfu4:
X mov bx,4+2(si) | extract exp
X movb dh,bh | extract sign
X shl bx,*1 | shift up one, then down 8
X movb bl,bh
X xorb bh,bh
X
X mov ax,4+2(si) | remove exponent from mantissa
X and ax,#0x7F
X test bx,bx | check for zero exponent - no leading "1"
X jz 0f | for denormalized numbers
X or ax,#0x80 | restore implied leading "1"
X0: mov 4+2(si),ax
X
X cmp bx,#BIAS4 | check exponent
X jl zer4 | strictly fractional, no integer part?
X cmp bx,#BIAS4+32 | is it too big to fit in a 32-bit integer?
X jg toobig
X
X1: cmp bx,#BIAS4+24 | shifted all the way down yet?
X jge 2f
X shr 4+2(si),*1 | shift down to align radix point;
X rcr 4+0(si),*1 | extra bits fall off the end (no rounding)
X inc bx
X jmp 1b
X
X2: cmp bx,#BIAS4+24 | do we have to shift up?
X jle 3f
X shl 4+0(si),*1 | shift up to align radix point
X rcl 4+2(si),*1
X dec bx
X jmp 2b
Xzer4:
X xor ax,ax | make the whole thing zero
X mov 4+0(si),ax
X mov 4+2(si),ax
X
X3: jmp cfu8b | amazingly, we can share the rest of the code
X
Xcfu8:
X cmp dx,#8
X jne ill | illegal EM op
X
X mov bx,4+6(si) | extract exp
X movb dh,bh | extract sign
X shr bx,*1
X shr bx,*1
X shr bx,*1
X shr bx,*1
X and bx,#2047 | kill sign bit
X
X mov ax,4+6(si) | remove exponent from mantissa
X and ax,#0x0F
X test bx,bx | check for zero exponent - no leading "1"
X jz 0f | for denormalized numbers
X or ax,#0x10 | restore implied leading "1"
X0: mov 4+6(si),ax
X
X cmp bx,#BIAS8 | check exponent
X jl zer8 | strictly fractional, no integer part?
X cmp bx,#BIAS8+32 | is it too big to fit in a 32-bit integer?
X jg toobig
X
X1: cmp bx,#BIAS8+53 | shifted all the way down yet?
X jge cfu8b
X shr 4+6(si),*1 | shift down to align radix point;
X rcr 4+4(si),*1 | extra bits fall off the end (no rounding)
X rcr 4+2(si),*1
X rcr 4+0(si),*1
X inc bx
X jmp 1b
Xzer8:
X xor ax,ax | make the whole thing zero
X mov 4+0(si),ax
X mov 4+2(si),ax
X mov 4+4(si),ax
X mov 4+6(si),ax
X
Xcfu8b: cmp cx,#2
X jne cfu8_4
Xcfu8_2: | 8-byte float to 2-byte integer
X mov ax,4+2(si)
X or ax,ax
X jnz toobig | too big to fit into a 16-bit integer?
X mov ax,4+0(si)
X testb dh,dh | is it negative?
X js toobig
X pop si | restore si
X pop bx | save program counter
X xorb dh,dh | dl still contains source size
X add sp,dx | get rid of source value
X push ax | put return value on stack
X push bx | put back return address
X ret
Xcfu8_4: | 8-byte float to 4-byte integer
X cmp cx,#4
X jne ill
X mov ax,4+0(si) | put integer into registers
X mov cx,4+2(si)
X testb dh,dh | is it negative?
X js toobig
X pop si | restore si
X pop bx | save program counter
X xorb dh,dh | dl still contains source size
X add sp,dx | get rid of source value
X push cx | put return value on stack
X push ax
X push bx | put back return address
X ret
Xill:
X pop bx
X jmp .trpilin
X
XECONV = 10
X .globl .fat
Xtoobig:
X mov ax,#ECONV
X push ax
X jmp .fat
END_OF_FILE
if test 3044 -ne `wc -c <'cfu.x'`; then
echo shar: \"'cfu.x'\" unpacked with wrong size!
fi
# end of 'cfu.x'
fi
if test -f 'mlf8.x' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'mlf8.x'\"
else
echo shar: Extracting \"'mlf8.x'\" \(5794 characters\)
sed "s/^X//" >'mlf8.x' <<'END_OF_FILE'
X.define .mlf8
X|
X| floating point multiply routine
X| author: Peter S. Housel 4/1/89
X|
X
XU = 14 | offset of multiplier
XV = 6 | offset of multiplicand
XBIAS8 = 0x3FF - 1
X
X .text
X .globl .mlf8
X.mlf8:
X push si
X push di
X mov si,sp | point to arguments
X
X mov bx,6+U(si) | extract u.exp
X movb dh,bh | extract u.sign
X shr bx,*1
X shr bx,*1
X shr bx,*1
X shr bx,*1
X and bx,#2047 | kill sign bit
X
X mov cx,6+V(si) | extract v.exp
X movb dl,ch | extract v.sign
X shr cx,*1
X shr cx,*1
X shr cx,*1
X shr cx,*1
X and cx,#2047 | kill sign bit
X
X mov ax,6+U(si) | remove exponent from u.mantissa
X and ax,#0x0F
X test bx,bx | check for zero exponent - no leading "1"
X jz 0f
X or ax,#0x10 | restore implied leading "1"
X0: mov 6+U(si),ax
X or ax,4+U(si)
X or ax,2+U(si)
X or ax,0+U(si)
X jz retz | multiplying by zero
X
X mov ax,6+V(si) | remove exponent from v.mantissa
X and ax,#0x0F
X test cx,cx | check for zero exponent - no leading "1"
X jz 0f
X or ax,#0x10 | restore implied leading "1"
X0: mov 6+V(si),ax
X or ax,4+V(si)
X or ax,2+V(si)
X or ax,0+V(si)
X jz retz2 | multiplying by zero
X
X xorb dh,dl | xor sign bits for resultant sign
X add bx,cx | add exponents,
X sub bx,#BIAS8+5 | remove excess bias, account for
X | repositioning
X
X xor ax,ax | initialize 128-bit product to zero
X push ax
X push ax
X push ax
X push ax
X push ax
X push ax
X push ax
X push ax
X mov di,sp | di => result
X push dx | save sign
X push bx | save exponent
X
X| see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
X
Xbig0: mov ax,0+U(si) | bigit from multiplier
X test ax,ax
X jz big1 | don't multiply by zero
X mov cx,ax | hold in cx to reduce memory reads
X mul 0+V(si)
X mov 0(di),ax | store into result
X mov bx,dx | save "carry" bigit
X
X mov ax,cx | restore multiplier bigit
X mul 2+V(si)
X| add ax,2(di) | add in bigit so far
X| adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 2(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 4+V(si)
X| add ax,4(di) | add in bigit so far
X| adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 4(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 6+V(si)
X| add ax,6(di) | add in bigit so far
X| adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 6(di),ax | store into result
X mov 8(di),dx | store carry
X|--
Xbig1: mov ax,2+U(si) | bigit from multiplier
X test ax,ax
X jz big2 | don't multiply by zero
X mov cx,ax | hold in cx to reduce memory reads
X mul 0+V(si)
X add ax,2(di) | add in bigit so far
X adc dx,#0
X mov 2(di),ax | store into result
X mov bx,dx | save "carry" bigit
X
X mov ax,cx | restore multiplier bigit
X mul 2+V(si)
X add ax,4(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 4(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 4+V(si)
X add ax,6(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 6(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 6+V(si)
X add ax,8(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 8(di),ax | store into result
X mov 10(di),dx | store carry
X|--
Xbig2: mov ax,4+U(si) | bigit from multiplier
X test ax,ax
X jz big3 | don't multiply by zero
X mov cx,ax | hold in cx to reduce memory reads
X mul 0+V(si)
X add ax,4(di) | add in bigit so far
X adc dx,#0
X mov 4(di),ax | store into result
X mov bx,dx | save "carry" bigit
X
X mov ax,cx | restore multiplier bigit
X mul 2+V(si)
X add ax,6(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 6(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 4+V(si)
X add ax,8(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 8(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 6+V(si)
X add ax,10(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 10(di),ax | store into result
X mov 12(di),dx | store carry
X|--
Xbig3: mov ax,6+U(si) | bigit from multiplier
X test ax,ax
X jz mdone | don't multiply by zero
X mov cx,ax | hold in cx to reduce memory reads
X mul 0+V(si)
X add ax,6(di) | add in bigit so far
X adc dx,#0
X mov 6(di),ax | store into result
X mov bx,dx | save "carry" bigit
X
X mov ax,cx | restore multiplier bigit
X mul 2+V(si)
X add ax,8(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 8(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 4+V(si)
X add ax,10(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 10(di),ax | store into result
X mov bx,dx
X
X mov ax,cx | restore multiplier bigit
X mul 6+V(si)
X add ax,12(di) | add in bigit so far
X adc dx,#0
X add ax,bx | add in carry bigit
X adc dx,#0
X mov 12(di),ax | store into result
X mov 14(di),dx | store carry - SHOULD BE ZERO
X|----
Xmdone: pop bx | restore exponent
X pop dx | restore sign
X
X2:
X mov ax,12(di) | multiply (shift) until
X and ax,#0xfff0 | one in "implied" position
X jnz 3f
X shl 0(di),*1
X rcl 2(di),*1
X rcl 4(di),*1
X rcl 6(di),*1
X rcl 8(di),*1
X rcl 10(di),*1
X rcl 12(di),*1
X dec bx | decrement exponent
X jmp 2b
X3:
X movb dl,5(di) | get rounding bits
X mov ax,0(di)
X orb al,4(di)
X or ax,2(di)
X jz 4f
X orb dl,*1 | set "sticky bit" if any low-order set
X4:
X add si,#U | si => u
X add di,#6 | di => most sig words of result
X xchg si,di | reverse for move into u's location
X mov
X mov
X mov
X mov
X add sp,#16 | remove result
X9: pop di | restore register variables
X pop si
X pop cx | save return address
X add sp,#8 | remove v
X push cx | restore return address
X jmp .norm8
X
Xretz2:
X xor ax,ax
X mov 0+U(si),ax
X mov 2+U(si),ax
X mov 4+U(si),ax
X mov 6+U(si),ax
X
Xretz: xor bx,bx | exponent = 0
X xor dx,dx | sign = 0, round bits = 0
X jmp 9b
END_OF_FILE
if test 5794 -ne `wc -c <'mlf8.x'`; then
echo shar: \"'mlf8.x'\" unpacked with wrong size!
fi
# end of 'mlf8.x'
fi
if test -f 'vfprintf.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'vfprintf.c'\"
else
echo shar: Extracting \"'vfprintf.c'\" \(15492 characters\)
sed "s/^X//" >'vfprintf.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1988 Regents of the University of California.
X * All rights reserved.
X *
X * Redistribution and use in source and binary forms are permitted
X * provided that the above copyright notice and this paragraph are
X * duplicated in all such forms and that any documentation,
X * advertising materials, and other materials related to such
X * distribution and use acknowledge that the software was developed
X * by the University of California, Berkeley. The name of the
X * University may not be used to endorse or promote products derived
X * from this software without specific prior written permission.
X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
X * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
X */
X
X#ifdef LIBC_SCCS
Xstatic char sccsid[] = "@(#)doprnt.c 5.35 (Berkeley) 6/27/88";
X#endif /* LIBC_SCCS */
X
X#include <sys/types.h>
X#include <stdarg.h>
X#include <stdio.h>
X#include <ctype.h>
X
Xtypedef unsigned char u_char;
Xtypedef unsigned long u_long;
X
X/* 11-bit exponent (VAX G floating point) is 308 decimal digits */
X#define MAXEXP 308
X/* 128 bit fraction takes up 39 decimal digits; max reasonable precision */
X#define MAXFRACT 39
X
X#define DEFPREC 6
X
X#define BUF (MAXEXP+MAXFRACT+1) /* + decimal point */
X
X#define PUTC(ch) (void) putc(ch, fp)
X
X#define ARG() \
X _ulong = flags&LONGINT ? va_arg(argp, long) : \
X flags&SHORTINT ? va_arg(argp, short) : va_arg(argp, int);
X
X#define todigit(c) ((c) - '0')
X#define tochar(n) ((n) + '0')
X
X#define LONGINT 0x01 /* long integer */
X#define LONGDBL 0x02 /* long double; unimplemented */
X#define SHORTINT 0x04 /* short integer */
X#define ALT 0x08 /* alternate form */
X#define LADJUST 0x10 /* left adjustment */
X#define ZEROPAD 0x20 /* zero (as opposed to blank) pad */
X#define HEXPREFIX 0x40 /* add 0x or 0X prefix */
X
Xint vfprintf(fp, fmt0, argp)
X register FILE *fp;
X u_char *fmt0;
X va_list argp;
X{
X register u_char *fmt; /* format string */
X register int ch; /* character from fmt */
X register int cnt; /* return value accumulator */
X register int n; /* random handy integer */
X register char *t; /* buffer pointer */
X double _double; /* double precision arguments %[eEfgG] */
X u_long _ulong; /* integer arguments %[diouxX] */
X int base; /* base for [diouxX] conversion */
X int dprec; /* decimal precision in [diouxX] */
X int fieldsz; /* field size expanded by sign, etc */
X int flags; /* flags as above */
X int fpprec; /* `extra' floating precision in [eEfgG] */
X int prec; /* precision from format (%.3d), or -1 */
X int realsz; /* field size expanded by decimal precision */
X int size; /* size of converted field or string */
X int width; /* width from format (%8d), or 0 */
X char sign; /* sign prefix (' ', '+', '-', or \0) */
X char softsign; /* temporary negative sign for floats */
X char *digs; /* digits for [diouxX] conversion */
X char buf[BUF]; /* space for %c, %[diouxX], %[eEfgG] */
X
X fmt = fmt0;
X digs = "0123456789abcdef";
X for (cnt = 0;; ++fmt) {
X if (!(ch = *fmt))
X return (cnt);
X if (ch != '%') {
X PUTC(ch);
X continue;
X }
X flags = 0; dprec = 0; fpprec = 0; width = 0;
X prec = -1;
X sign = '\0';
X
Xrflag: switch (*++fmt) {
X case ' ':
X /*
X * ``If the space and + flags both appear, the space
X * flag will be ignored.''
X * -- ANSI X3J11
X */
X if (!sign)
X sign = ' ';
X goto rflag;
X case '#':
X flags |= ALT;
X goto rflag;
X case '*':
X /*
X * ``A negative field width argument is taken as a
X * - flag followed by a positive field width.''
X * -- ANSI X3J11
X * They don't exclude field widths read from args.
X */
X if ((width = va_arg(argp, int)) >= 0)
X goto rflag;
X width = -width;
X /* FALLTHROUGH */
X case '-':
X flags |= LADJUST;
X goto rflag;
X case '+':
X sign = '+';
X goto rflag;
X case '.':
X if (*++fmt == '*')
X n = va_arg(argp, int);
X else {
X n = 0;
X while (isascii(*fmt) && isdigit(*fmt))
X n = 10 * n + todigit(*fmt++);
X --fmt;
X }
X prec = n < 0 ? -1 : n;
X goto rflag;
X case '0':
X /*
X * ``Note that 0 is taken as a flag, not as the
X * beginning of a field width.''
X * -- ANSI X3J11
X */
X flags |= ZEROPAD;
X goto rflag;
X case '1': case '2': case '3': case '4':
X case '5': case '6': case '7': case '8': case '9':
X n = 0;
X do {
X n = 10 * n + todigit(*fmt);
X } while (isascii(*++fmt) && isdigit(*fmt));
X width = n;
X --fmt;
X goto rflag;
X case 'L':
X flags |= LONGDBL;
X goto rflag;
X case 'h':
X flags |= SHORTINT;
X goto rflag;
X case 'l':
X flags |= LONGINT;
X goto rflag;
X case 'c':
X *(t = buf) = va_arg(argp, int);
X size = 1;
X sign = '\0';
X goto pforw;
X case 'D':
X flags |= LONGINT;
X /*FALLTHROUGH*/
X case 'd':
X case 'i':
X ARG();
X if ((long)_ulong < 0) {
X _ulong = -_ulong;
X sign = '-';
X }
X base = 10;
X goto number;
X#ifdef FLOATLIB
X case 'e':
X case 'E':
X case 'f':
X case 'g':
X case 'G':
X _double = va_arg(argp, double);
X /*
X * don't do unrealistic precision; just pad it with
X * zeroes later, so buffer size stays rational.
X */
X if (prec > MAXFRACT) {
X if (*fmt != 'g' && *fmt != 'G' || (flags&ALT))
X fpprec = prec - MAXFRACT;
X prec = MAXFRACT;
X }
X else if (prec == -1)
X prec = DEFPREC;
X /*
X * softsign avoids negative 0 if _double is < 0 and
X * no significant digits will be shown
X */
X if (_double < 0) {
X softsign = '-';
X _double = -_double;
X }
X else
X softsign = 0;
X /*
X * cvt may have to round up past the "start" of the
X * buffer, i.e. ``intf("%.2f", (double)9.999);'';
X * if the first char isn't NULL, it did.
X */
X *buf = NULL;
X size = cvt(_double, prec, flags, &softsign, *fmt, buf,
X buf + sizeof(buf));
X if (softsign)
X sign = '-';
X t = *buf ? buf : buf + 1;
X goto pforw;
X#endif FLOATLIB
X case 'n':
X if (flags & LONGINT)
X *va_arg(argp, long *) = cnt;
X else if (flags & SHORTINT)
X *va_arg(argp, short *) = cnt;
X else
X *va_arg(argp, int *) = cnt;
X break;
X case 'O':
X flags |= LONGINT;
X /*FALLTHROUGH*/
X case 'o':
X ARG();
X base = 8;
X goto nosign;
X case 'p':
X /*
X * ``The argument shall be a pointer to void. The
X * value of the pointer is converted to a sequence
X * of printable characters, in an implementation-
X * defined manner.''
X * -- ANSI X3J11
X */
X /* NOSTRICT */
X _ulong = (u_long)va_arg(argp, void *);
X base = 16;
X goto nosign;
X case 's':
X if (!(t = va_arg(argp, char *)))
X t = "(null)";
X if (prec >= 0) {
X /*
X * can't use strlen; can only look for the
X * NUL in the first `prec' characters, and
X * strlen() will go further.
X */
X char *p, *memchr();
X
X if (p = memchr(t, 0, prec)) {
X size = p - t;
X if (size > prec)
X size = prec;
X } else
X size = prec;
X } else
X size = strlen(t);
X sign = '\0';
X goto pforw;
X case 'U':
X flags |= LONGINT;
X /*FALLTHROUGH*/
X case 'u':
X ARG();
X base = 10;
X goto nosign;
X case 'X':
X digs = "0123456789ABCDEF";
X /* FALLTHROUGH */
X case 'x':
X ARG();
X base = 16;
X /* leading 0x/X only if non-zero */
X if (flags & ALT && _ulong != 0)
X flags |= HEXPREFIX;
X
X /* unsigned conversions */
Xnosign: sign = '\0';
X /*
X * ``... diouXx conversions ... if a precision is
X * specified, the 0 flag will be ignored.''
X * -- ANSI X3J11
X */
Xnumber: if ((dprec = prec) >= 0)
X flags &= ~ZEROPAD;
X
X /*
X * ``The result of converting a zero value with an
X * explicit precision of zero is no characters.''
X * -- ANSI X3J11
X */
X t = buf + BUF;
X if (_ulong != 0 || prec != 0) {
X do {
X *--t = digs[_ulong % base];
X _ulong /= base;
X } while (_ulong);
X digs = "0123456789abcdef";
X if (flags & ALT && base == 8 && *t != '0')
X *--t = '0'; /* octal leading 0 */
X }
X size = buf + BUF - t;
X
Xpforw:
X /*
X * All reasonable formats wind up here. At this point,
X * `t' points to a string which (if not flags&LADJUST)
X * should be padded out to `width' places. If
X * flags&ZEROPAD, it should first be prefixed by any
X * sign or other prefix; otherwise, it should be blank
X * padded before the prefix is emitted. After any
X * left-hand padding and prefixing, emit zeroes
X * required by a decimal [diouxX] precision, then print
X * the string proper, then emit zeroes required by any
X * leftover floating precision; finally, if LADJUST,
X * pad with blanks.
X */
X
X /*
X * compute actual size, so we know how much to pad
X * fieldsz excludes decimal prec; realsz includes it
X */
X fieldsz = size + fpprec;
X if (sign)
X fieldsz++;
X if (flags & HEXPREFIX)
X fieldsz += 2;
X realsz = dprec > fieldsz ? dprec : fieldsz;
X
X /* right-adjusting blank padding */
X if ((flags & (LADJUST|ZEROPAD)) == 0 && width)
X for (n = realsz; n < width; n++)
X PUTC(' ');
X /* prefix */
X if (sign)
X PUTC(sign);
X if (flags & HEXPREFIX) {
X PUTC('0');
X PUTC((char)*fmt);
X }
X /* right-adjusting zero padding */
X if ((flags & (LADJUST|ZEROPAD)) == ZEROPAD)
X for (n = realsz; n < width; n++)
X PUTC('0');
X /* leading zeroes from decimal precision */
X for (n = fieldsz; n < dprec; n++)
X PUTC('0');
X
X /* the string or number proper */
X for (n = size; --n >= 0; )
X PUTC(*t++);
X /* trailing f.p. zeroes */
X while (--fpprec >= 0)
X PUTC('0');
X /* left-adjusting padding (always blank) */
X if (flags & LADJUST)
X for (n = realsz; n < width; n++)
X PUTC(' ');
X /* finally, adjust cnt */
X cnt += width > realsz ? width : realsz;
X break;
X case '\0': /* "%?" prints ?, unless ? is NULL */
X return (cnt);
X default:
X PUTC((char)*fmt);
X cnt++;
X }
X }
X /* NOTREACHED */
X}
X
X#ifdef FLOATLIB
Xstatic
Xcvt(number, prec, flags, signp, fmtch, startp, endp)
X double number;
X register int prec;
X int flags;
X u_char fmtch;
X char *signp, *startp, *endp;
X{
X register char *p, *t;
X register double fract;
X int dotrim, expcnt, gformat;
X double integer, tmp, modf();
X char *exponent(), *round();
X
X dotrim = expcnt = gformat = 0;
X fract = modf(number, &integer);
X
X /* get an extra slot for rounding. */
X t = ++startp;
X
X /*
X * get integer portion of number; put into the end of the buffer; the
X * .01 is added for modf(356.0 / 10, &integer) returning .59999999...
X */
X for (p = endp - 1; integer; ++expcnt) {
X tmp = modf(integer / 10, &integer);
X *p-- = tochar((int)((tmp + .01) * 10));
X }
X switch(fmtch) {
X case 'f':
X /* reverse integer into beginning of buffer */
X if (expcnt)
X for (; ++p < endp; *t++ = *p);
X else
X *t++ = '0';
X /*
X * if precision required or alternate flag set, add in a
X * decimal point.
X */
X if (prec || flags&ALT)
X *t++ = '.';
X /* if requires more precision and some fraction left */
X if (fract) {
X if (prec)
X do {
X fract = modf(fract * 10, &tmp);
X *t++ = tochar((int)tmp);
X } while (--prec && fract);
X if (fract)
X startp = round(fract, (int *)NULL, startp,
X t - 1, (char)0, signp);
X }
X for (; prec--; *t++ = '0');
X break;
X case 'e':
X case 'E':
Xeformat: if (expcnt) {
X *t++ = *++p;
X if (prec || flags&ALT)
X *t++ = '.';
X /* if requires more precision and some integer left */
X for (; prec && ++p < endp; --prec)
X *t++ = *p;
X /*
X * if done precision and more of the integer component,
X * round using it; adjust fract so we don't re-round
X * later.
X */
X if (!prec && ++p < endp) {
X fract = 0;
X startp = round((double)0, &expcnt, startp,
X t - 1, *p, signp);
X }
X /* adjust expcnt for digit in front of decimal */
X --expcnt;
X }
X /* until first fractional digit, decrement exponent */
X else if (fract) {
X /* adjust expcnt for digit in front of decimal */
X for (expcnt = -1;; --expcnt) {
X fract = modf(fract * 10, &tmp);
X if (tmp)
X break;
X }
X *t++ = tochar((int)tmp);
X if (prec || flags&ALT)
X *t++ = '.';
X }
X else {
X *t++ = '0';
X if (prec || flags&ALT)
X *t++ = '.';
X }
X /* if requires more precision and some fraction left */
X if (fract) {
X if (prec)
X do {
X fract = modf(fract * 10, &tmp);
X *t++ = tochar((int)tmp);
X } while (--prec && fract);
X if (fract)
X startp = round(fract, &expcnt, startp,
X t - 1, (char)0, signp);
X }
X /* if requires more precision */
X for (; prec--; *t++ = '0');
X
X /* unless alternate flag, trim any g/G format trailing 0's */
X if (gformat && !(flags&ALT)) {
X while (t > startp && *--t == '0');
X if (*t == '.')
X --t;
X ++t;
X }
X t = exponent(t, expcnt, fmtch);
X break;
X case 'g':
X case 'G':
X /* a precision of 0 is treated as a precision of 1. */
X if (!prec)
X ++prec;
X /*
X * ``The style used depends on the value converted; style e
X * will be used only if the exponent resulting from the
X * conversion is less than -4 or greater than the precision.''
X * -- ANSI X3J11
X */
X if (expcnt > prec || !expcnt && fract && fract < .0001) {
X /*
X * g/G format counts "significant digits, not digits of
X * precision; for the e/E format, this just causes an
X * off-by-one problem, i.e. g/G considers the digit
X * before the decimal point significant and e/E doesn't
X * count it as precision.
X */
X --prec;
X fmtch -= 2; /* G->E, g->e */
X gformat = 1;
X goto eformat;
X }
X /*
X * reverse integer into beginning of buffer,
X * note, decrement precision
X */
X if (expcnt)
X for (; ++p < endp; *t++ = *p, --prec);
X else
X *t++ = '0';
X /*
X * if precision required or alternate flag set, add in a
X * decimal point. If no digits yet, add in leading 0.
X */
X if (prec || flags&ALT) {
X dotrim = 1;
X *t++ = '.';
X }
X else
X dotrim = 0;
X /* if requires more precision and some fraction left */
X if (fract) {
X if (prec) {
X do {
X fract = modf(fract * 10, &tmp);
X *t++ = tochar((int)tmp);
X } while(!tmp);
X while (--prec && fract) {
X fract = modf(fract * 10, &tmp);
X *t++ = tochar((int)tmp);
X }
X }
X if (fract)
X startp = round(fract, (int *)NULL, startp,
X t - 1, (char)0, signp);
X }
X /* alternate format, adds 0's for precision, else trim 0's */
X if (flags&ALT)
X for (; prec--; *t++ = '0');
X else if (dotrim) {
X while (t > startp && *--t == '0');
X if (*t != '.')
X ++t;
X }
X }
X return(t - startp);
X}
X
Xstatic char *
Xround(fract, exp, start, end, ch, signp)
X double fract;
X int *exp;
X register char *start, *end;
X char ch, *signp;
X{
X double tmp;
X
X if (fract)
X (void)modf(fract * 10, &tmp);
X else
X tmp = todigit(ch);
X if (tmp > 4)
X for (;; --end) {
X if (*end == '.')
X --end;
X if (++*end <= '9')
X break;
X *end = '0';
X if (end == start) {
X if (exp) { /* e/E; increment exponent */
X *end = '1';
X ++*exp;
X }
X else { /* f; add extra digit */
X *--end = '1';
X --start;
X }
X break;
X }
X }
X /* ``"%.3f", (double)-0.0004'' gives you a negative 0. */
X else if (*signp == '-')
X for (;; --end) {
X if (*end == '.')
X --end;
X if (*end != '0')
X break;
X if (end == start)
X *signp = 0;
X }
X return(start);
X}
X
Xstatic char *
Xexponent(p, exp, fmtch)
X register char *p;
X register int exp;
X u_char fmtch;
X{
X register char *t;
X char expbuf[MAXEXP];
X
X *p++ = fmtch;
X if (exp < 0) {
X exp = -exp;
X *p++ = '-';
X }
X else
X *p++ = '+';
X t = expbuf + MAXEXP;
X if (exp > 9) {
X do {
X *--t = tochar(exp % 10);
X } while ((exp /= 10) > 9);
X *--t = tochar(exp);
X for (; t < expbuf + MAXEXP; *p++ = *t++);
X }
X else {
X *p++ = '0';
X *p++ = tochar(exp);
X }
X return(p);
X}
X#endif FLOATLIB
END_OF_FILE
if test 15492 -ne `wc -c <'vfprintf.c'`; then
echo shar: \"'vfprintf.c'\" unpacked with wrong size!
fi
# end of 'vfprintf.c'
fi
echo shar: End of archive 1 \(of 3\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 2 3 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 3 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0