Hi Folks Eric and I have just posted a new script to the REBOL script library. You can download it directly from: http://www.rebol.org/math/decimal.r The purpose of this package of functions is to improve REBOL'S support for floating-point mathematics. It provides many of the capabilities present in languages like C and FORTRAN, but which are currently unavailable in REBOL, such as native binary I/O. The package is the result of a collaboration between Eric Long, myself, and Gerald Goertzel. Comments, questions and suggestions are welcome. We are especially interested in whether these functions (which have only been tested with the Windows version of REBOL) produce the same results on other platforms and non-Intel CPU's. Please address any comments to the REBOL mail-list or send them directly to me at [EMAIL PROTECTED] The package is heavily documented with comments and examples. The functionality is encapsulated in just three objects: REAL, IEEE, and EMPIRICAL. These are the only names that are added to the global context. The REAL object--------------------- This object contains many functions for the manipulation and examination of REBOL decimal values. Among other things, these functions provide full support for native binary floating-point file IO, compatible with C, FORTRAN, etc. REBOL decimal and money values may also be saved transparently in a native binary hex-string and loaded without incurring roundoff error. These functions are necessary because the standard REBOL representation of decimal values may round off one or two significant digits, and the rounding of money values is an even greater source of error. Outline of functions in the REAL object--------------------- Main Interface Functions FORM Returns precise decimal representation (requires format.r). SHOW Returns a string showing the bits in the IEEE representation. TO-NATIVE Converts a decimal value into a native IEEE binary FROM-NATIVE Converts an 8-byte native binary into a decimal value. SAVE Saves a value with all decimal values converted. LOAD Loads a value and restores 8-byte binaries to decimals. WRITE Writes a block of decimal values to a file in native form. READ Reads a file containing a series of decimal values expressed in native form. Helper Functions SPLIT Returns a block with the three components of the IEEE double floating point representation of a decimal value. CONVERT Converts all decimal values to native binaries, leaves other values untouched. RESTORE Restores 8-byte binaries to decimal values. TO-BIN Returns a binary string representation or a numeric value. FROM-BIN Returns the integer represented by string of 1's and 0's. TO-MATRIX Converts a flat block into nested blocks of any depth. Most of the functions come in pairs: SAVE <--> LOAD WRITE <--> READ TO-NATIVE <--> FROM-NATIVE CONVERT <--> RESTORE TO-BIN <--> FROM/BIN These pairs have the property of reversing exactly (i.e., FROM-BIN TO-BIN X yields X). The IEEE and EMPIRICAL objects--------------- These provide several utility functions useful for more in-depth exploration of the details of floating-point numbers in REBOL and testing the REBOL comparison functions. The two objects provide similar functionality but differ in implementation. In the IEEE object all of the values are calculated according to the IEEE754 standard for double-precision floating-point. In the EMPIRICAL object, the values are calculated using an iterative empirical method (similar to root-finding). Most of the calculations in both are done in binary. Outline of functions in the IEEE and ENPIRICAL objects------------ REBTEST A function that tests the sensitivity of the REBOL EQUAL? comparison function. GET-NEXT A function that makes the minimum increment to a decimal value. GET-LAST A function that makes the maximum decrement to a decimal value. MAX-REAL The largest decimal value. MIN-NORM The smallest positive normalized decimal. (IEEE only) MIN-REAL The smallest positive decimal value. EPS The minimum increment to 1 (smallest x such that (x+1)>1). REBTEST illustrates some flaws in the implementation of the REBOL comparison functions (as of version 2.2), but note that it is easy to get exact comparisons by using ZERO? NEGATIVE? and POSITIVE? to test the results of subtraction. Using FORMAT Several of the functions can optionally use the full-precision formatting available in Eric's FORMAT function. These functions are: real/form real/show/decimal IEEE/rebtest empirical/rebtest You can download FORMAT at http://www.rebol.org/math/format.r If you wish to use it in conjunction with DECIMAL, you should load format.r before loading decimal.r. Note: None of the examples given in this post make use of FORMAT. There is way too much in this package to describe in this post, so I will just present a few examples below. If there is interest, I am willing to write a few short tutorial posts about floating-point in REBOL. Have fun exploring REBOL numbers:-) Larry EXAMPLES ================================ REAL/SPLIT calculates the three components of the IEEE representation of a double floating point value. These are the actual values used by the CPU in numerical calculations. >> set [sign exponent fraction] real/split pi == [0 1024 2.57063812465794E+15] >> (-1 ** sign) * (2 ** (exponent - 1023)) * (1 + (fraction * (2 ** -52))) == 3.14159265358979 Numbers smaller in absolute magnitude than 2 ** -1022 (denormals, which have an exponent component of zero) use a different formula: >> set [sign exponent fraction] real/split probe 2 ** -1030 8.69169475979376E-311 == [0 0 17592186044416] >> (-1 ** sign) * (2 ** -1074) * fraction == 8.69169475979376E-311 REAL/SHOW returns a bit-string which shows the internal IEEE754 double format. We have separated the 3 components with spaces for clarity. The format is 64 bits wide with the leading bit being a sign bit (0 pos and 1 neg). The next 11 bits are the base-2 biased exponent. It is biased by adding 1023 to the actual exponent so that internally all exponents are non-negative. The remaining 52 bits are the fractional part of the binary mantissa. There are two types of doubles: normalized and denormalized. The normalized values cover the range from 2 ** -1022 up to +Inf. The denormalized values are equally spaced between 0 and 2 ** -1022 and have an exponent of 0. The normalized values have an implicit binary point and leading bit of 1 (i.e., the mantissa 010... represents 1.010...). >> real/show ieee/min-real == {0 00000000000 0000000000000000000000000000000000000000000000000001} >> real/show 0 == {0 00000000000 0000000000000000000000000000000000000000000000000000} >> real/show negate ieee/min-real == {1 00000000000 0000000000000000000000000000000000000000000000000001} >> real/show ieee/get-next 1 == {0 01111111111 0000000000000000000000000000000000000000000000000001} >> real/show 1 == {0 01111111111 0000000000000000000000000000000000000000000000000000} >> real/show ieee/get-last 1 == {0 01111111110 1111111111111111111111111111111111111111111111111111} >> real/show ieee/max-real == {0 11111111110 1111111111111111111111111111111111111111111111111111} REAL/TO-NATIVE returns the binary representation of a double. On some machines(Intel and Motorola 68xxx), the 8 bytes are stored in RAM and on disk with the bytes reversed (i.e., the least significant byte is the first one encountered in memory). Reading in a file of doubles of unknown format may require trying both reversed and and unreversed formats. >> real/to-native pi == #{400921FB54442D18} >> real/to-native/rev pi == #{182D4454FB210940} ; configuration used in PC's >> real/from-native real/to-native pi == 3.14159265358979 >> pi - real/from-native real/to-native pi == 0 ; converts back to exact value REAL/WRITE and REAL/READ are meant for writing and reading decimal values in native binary format. REAL/WRITE accepts a single numerical value or a block of numerical values that may be nested to any depth, but the return value of REAL/READ is always a flat block. >> b1: [] for x 1 16 1 [append b1 square-root x] == [1 1.4142135623731 1.73205080756888 2 2.23606797749979 ... >> real/write %test.bin b1 16 decimal values written >> bb1: real/read %test.bin == [1 1.4142135623731 1.73205080756888 2 2.23606797749979 ... >> b1/7 - probe bb1/7 2.64575131106459 ; square root of 7 == 0 ; is restored exactly >> b2: real/to-matrix b1 4 ; convert to 4x4 matrix == [[1 1.4142135623731 1.73205080756888 2] [2.23606797749979 ... >> real/write %test.bin b2 16 decimal values written ; the same values are written, >> bb2: real/read %test.bin ; and read back into a flat block == [1 1.4142135623731 1.73205080756888 2 2.23606797749979 ... REAL/SAVE and REAL/LOAD are meant to be completely compatible with the standard functions SAVE and LOAD, except that there will be no roundoff error from saving and loading decimal or money values. >> amount: [100 * $1.504] == [100 * $1.50] >> save %test-x.dat amount >> do load %test-x.dat == $150.00 ; large round-off error >> real/save %test-y.dat amount >> do real/load %test-y.dat == $150.40 ; no error IEEE/REBTEST makes ten successive minimum increments and decrements to the target value, and tests whether EQUAL? returns TRUE when comparing the results to the original value. In all cases EQUAL? should return false. This example is done without having FORMAT loaded. Let's use 2 ** 53 as the target. All positive integers from 1 to 2 ** 53 are exactly represented as binary doubles. All doubles greater than 2 ** 52 are integers. The doubles from 2 ** 52 thru 2 ** 53 are consecutive integers. The first column shows the result of REBOL's test for equality of the target to (target + increment). The 2nd column shows the REBOL display value and the 3rd column shows the increment. >> ieee/rebtest 2 ** 53 ; the exact value is 9007199254740992 target 9.00719925474099E+15 ; REBOL display rounds to 15 sig. digits positive increments t=t+i target+increment increment true 9.00719925474099E+15 2 ; positive increments are of size 2 true 9.007199254741E+15 4 ; REBOL display value changes true 9.007199254741E+15 6 ; REBOL equality test is wrong true 9.007199254741E+15 8 true 9.007199254741E+15 10 true 9.007199254741E+15 12 true 9.00719925474101E+15 14 ; REBOL display changes again true 9.00719925474101E+15 16 ; but equality test is still wrong false 9.00719925474101E+15 18 ; equality test correct once we add 18 false 9.00719925474101E+15 20 negative increments t=t+i target+increment increment false 9.00719925474099E+15 -1 ; negative increments are of size 1 false 9.00719925474099E+15 -2 false 9.00719925474099E+15 -3 false 9.00719925474099E+15 -4 false 9.00719925474099E+15 -5 false 9.00719925474099E+15 -6 ; REBOL equality test OK for all inc. false 9.00719925474099E+15 -7 false 9.00719925474098E+15 -8 ; REBOL display changes false 9.00719925474098E+15 -9 false 9.00719925474098E+15 -10 Other interesting targets are: 0, 1, ieee/max-real, ieee/min-real, 2 ** 52, and ieee/min-norm. >> empirical/rebtest 2 ** 53 Increments: Smallest NOT EQUAL? Factor Positive: 2 18 9 ; EQUAL? inc 9 times larger than true Negative: -1 -1 1
