{ ==================================================================== sfp.f Software floating point package SwiftForth Copyright (c) Dr. Heinrich Hohl, Munich, Germany Originally written using LMI PC/FORTH 3.2 (Fall 1992) Rewritten for SwiftForth 3.5.6 (December 2014) -------------------------------------------------------------------- The floating point numbers used in this package are represented by a composite single length number. 24 bits are used for the significand, and the top 8 bits for the exponent. Significand: 24 bits (800000 ... FFFFFF 000000 000001 ... 7FFFFF) Decimal range: -8,388,608 ... +8,388,607 Exponent: 8 bits (80 ... FF 00 01 ... 7F) Decimal range -128 ... +127 The significand must be normalized in order to obtain a valid floating point number. It contains an implied decimal point five digits from the right. Valid significands are: Negative: -9.99999 ... -1.00000 Zero: 0.00000 Positive: +1.00000 ... +9.99999 This means that floating point numbers have a resolution of 6 decimal digits in total. Zero is always given by 0-sig 0-exp. Floating point numbers are handled like single length numbers with respect to stack operations, storage etc. Floating point conversion is appended to the standard number conversion routine. NUMBER is extended, NUMBER? remains unchanged. Decimal base is mandatory for input and output of fp numbers. Any base may be used for input and output of integers. However, this may cause confusion if numbers contain an "E". Decimal base: 4.38E2 --> fp number (4.38 x 10^2) Hex base: 4.38E2 --> d number (438E2 0) ==================================================================== } FORTH DEFINITIONS DECIMAL PACKAGE PRIVATE-SFP { -------------------------------------------------------------------- Basic data structures -------------------------------------------------------------------- } PUBLIC AKA CELLS FLOATS AKA CELL+ FLOAT+ \ size of a floating point number is 32 bits = 4 bytes AKA LITERAL FLITERAL ( f -- ) IMMEDIATE AKA CONSTANT FCONSTANT ( f -- ) AKA VARIABLE FVARIABLE ( -- ) AKA @ F@ ( addr -- f ) AKA ! F! ( f addr -- ) { -------------------------------------------------------------------- Supportive variables -------------------------------------------------------------------- } PRIVATE VARIABLE fpdpl \ for temporary storage of significand's decimal point location VARIABLE fpexp \ for temporary storage of an exponent VARIABLE fpsign \ for temporary storage of a sign (significand or exponent) 2VARIABLE fpstring \ for temporary storage of original string (addr len) { -------------------------------------------------------------------- Pack and unpack floting point numbers -------------------------------------------------------------------- } HEX : 8>32 ( n -- n') DUP 00000080 AND IF FFFFFF00 OR THEN ; \ sign-extend an 8 bit number to a 32 bit number : 24>32 ( n -- n') DUP 00800000 AND IF FF000000 OR THEN ; \ sign-extend a 24 bit number to a 32 bit number : PACK ( n-sig n-exp -- f) 0018 LSHIFT \ shift exponent 24 bits to the left [ee00|0000] SWAP 00FFFFFF AND \ clear top 8 bits of significand [00ss|ssss] OR ; \ combine both numbers to [eess|ssss] \ pack significand and exponent to a floating point number : UNPACK ( f -- n-sig n-exp) DUP 00FFFFFF AND \ remove exponent from fp number [00ss|ssss] 24>32 \ sign-extend significand SWAP 0018 RSHIFT \ shift exponent 24 bits to the right [0000|00ee] 8>32 ; \ sign-extend exponent \ unpack floating point number to significand and exponent DECIMAL { -------------------------------------------------------------------- Decimal multiplication and division -------------------------------------------------------------------- } : D10* ( d1 -- d2) 2DUP D2* D2* D+ D2* ; \ multiply double length number by 10 : D10/ ( d1 -- d2) S>D 10 UT/ ; \ divide double length number by 10 { -------------------------------------------------------------------- Numeric input words -------------------------------------------------------------------- } : DEFLATE ( ud -- ud') BEGIN 9.99999 2OVER D< WHILE D10/ 1 fpexp +! REPEAT ; \ deflate significand greater than 9.99999 and adjust exponent : INFLATE ( u -- u') BEGIN DUP 100000 < WHILE 10 * -1 fpexp +! REPEAT ; \ inflate significand smaller than 1.00000 and adjust exponent : NORMALIZE ( d-sig n-exp -- n-sig' n-exp') fpexp ! \ store original exponent 2DUP D0= \ check significand IF \ d-sig zero: return 0 0 ELSE \ d-sig nonzero: normalize DUP 0< >R \ save sign of significand DABS \ prepare for unsigned routines DEFLATE \ deflate significand if necessary D>S \ single length is now sufficient INFLATE \ inflate significand if necessary R> IF NEGATE THEN \ restore sign of significand fpexp @ \ fetch adjusted exponent THEN ; \ normalize significand and adjust exponent accordingly { -------------------------------------------------------------------- } PUBLIC : D>F ( d -- f) 5 NORMALIZE PACK ; \ convert double length number to float (ignoring DPL) : S>F ( n -- f) S>D D>F ; \ convert single length number to float { -------------------------------------------------------------------- } PRIVATE : ( addr len -- addr' len') OVER C@ DUP [CHAR] - = DUP fpsign ! SWAP [CHAR] + = OR IF 1 /STRING THEN ; \ handle leading plus or minus character : <.DIGITS> ( ud addr len -- ud' addr' len') OVER C@ [CHAR] . = IF 1 /STRING DUP >R >NUMBER R> OVER - fpdpl ! THEN ; \ handle decimal point and following digits : APPLY-SIGN ( ud addr len -- d addr len) fpsign @ IF 2SWAP DNEGATE 2SWAP THEN ; \ apply sign to converted number : ( 0. addr len -- d addr' len') -1 fpdpl ! DUP -EXIT >NUMBER DUP IF <.DIGITS> THEN APPLY-SIGN ; \ convert significand and return remaining string : ( 0. addr len -- d addr' len') DUP -EXIT OVER C@ [CHAR] E = IF 1 /STRING DUP IF >NUMBER APPLY-SIGN THEN THEN ; \ convert exponent and return remaining string : ( 0. addr len -- d addr' len') DUP -EXIT OVER C@ UPPER DUP [CHAR] D = SWAP [CHAR] E = OR IF 1 /STRING THEN DUP IF >NUMBER APPLY-SIGN THEN ; \ same as above, but accept {D, d, E, e, +, -} before exponent : FLOAT ( d-sig n-exp -- f) 5 + DPL @ 0 MAX - NORMALIZE PACK ; \ convert to float considering significand's DPL { -------------------------------------------------------------------- } : FCONVERT ( addr len -- f xt) \ float conversion successful; xt of FLITERAL ( addr len -- addr len 0) \ float conversion not successful 2DUP fpstring 2! \ save original string 0. 2SWAP \ d-sig addr' len' 0. 2SWAP \ d-sig d-exp addr" len" NIP \ d-sig d-exp len" IF \ d-sig d-exp; conversion failed 2DROP 2DROP fpstring 2@ 0 \ return original addr len 0 ELSE \ d-sig d-exp; conversion successful D>S \ d-sig n-exp fpdpl @ DPL ! \ restore significand's DPL FLOAT \ convert to floating point number ['] FLITERAL \ return f xt THEN ; \ convert string to floating point number \ DPL is based on significand conversion (-1 if no decimal point) : FNUMBER? ( xt -- xt) \ if xt on stack ( addr len 0 -- f xt ) \ float conversion successful; xt of FLITERAL ( addr len 0 -- addr len 0) \ float conversion not successful ?DUP ?EXIT \ exit if an xt is on stack BASE @ 10 = \ is base decimal? IF FCONVERT ELSE 0 THEN ; ' FNUMBER? NUMBER-CONVERSION FLOAT ( addr len -- 0 | f -1) -TRAILING 0. 2SWAP \ d-sig addr' len' 0. 2SWAP \ d-sig d-exp addr" len" NIP \ d-sig d-exp len" IF \ d-sig d-exp; conversion failed 2DROP 2DROP 0 \ return 0 ELSE \ d-sig d-exp; conversion successful D>S \ d-sig n-exp fpdpl @ DPL ! \ restore significand's DPL FLOAT -1 \ return f -1 THEN ; \ convert string to floating point number \ special case: blank string returns fp zero \ DPL is based on significand conversion (-1 if no decimal point) { -------------------------------------------------------------------- Numeric output words -------------------------------------------------------------------- } PRIVATE : DSHIFT ( d n -- d') DUP 0< IF NEGATE 0 DO D10/ LOOP ELSE 0 ?DO D10* LOOP THEN ; \ perform decimal right shift (n<0) or left shift (n>0) PUBLIC : F>D ( f -- d) UNPACK SWAP S>D ROT 5 - DSHIFT ; \ convert floating point number to double length number : F>S ( f -- n) F>D D>S ; \ convert floating point number to single length number { -------------------------------------------------------------------- } PRIVATE VARIABLE f#places \ number of digits to show behind the decimal point : (FD.) ( fd -- addr len) SWAP OVER DABS <# f#places @ 0 ?DO # LOOP [CHAR] . HOLD #S ROT SIGN #> ; \ convert fixed point double length number to formatted string : ?DECIMAL ( -- ) BASE @ 10 <> ABORT" base must be decimal" ; { -------------------------------------------------------------------- } PUBLIC : PLACES ( n -- ) 0 MAX f#places ! ; \ specify the number of digits behind the decimal point 5 PLACES \ default setting based on significand's number range { -------------------------------------------------------------------- } : (F.) ( f -- addr len) ?DECIMAL UNPACK fpexp ! S>D f#places @ 5 - fpexp @ + DSHIFT (FD.) ; \ convert f to a formatted string using fixed-point notation : F. ( f -- ) (F.) TYPE SPACE ; \ display f in fixed format : F.R ( f width -- ) >R (F.) R> OVER - SPACES TYPE ; \ display f right justified in a field of specified width { -------------------------------------------------------------------- } : (FS.) ( f -- addr len) ?DECIMAL UNPACK fpexp ! S>D f#places @ 5 - DSHIFT (FD.) PAD PLACE S" E" PAD APPEND fpexp @ (.) PAD APPEND PAD COUNT ; \ convert f to a formatted string using scientific notation : FS. ( f -- ) (FS.) TYPE SPACE ; \ display f in scientific format : FS.R ( f width -- ) >R (FS.) R> OVER - SPACES TYPE ; \ display f right justified in a field of specified width { -------------------------------------------------------------------- Basic stack operators -------------------------------------------------------------------- } AKA DROP FDROP ( f1 f2 -- f1) AKA DUP FDUP ( f -- f f ) AKA OVER FOVER ( f1 f2 -- f1 f2 f3) AKA SWAP FSWAP ( f1 f2 -- f2 f1) AKA ROT FROT ( f1 f2 f3 -- f2 f3 f1) AKA 2DROP F2DROP AKA 2DUP F2DUP AKA 2OVER F2OVER AKA 2SWAP F2SWAP { -------------------------------------------------------------------- Comparison /1/ -------------------------------------------------------------------- } AKA 0= F0= ( f -- ?) AKA 0<> F0<> ( f -- ?) AKA = F= ( f1 f2 -- ?) : F0< ( f -- ?) UNPACK DROP 0< ; : F0> ( f -- ?) UNPACK DROP 0> ; { -------------------------------------------------------------------- Basic arithmetic -------------------------------------------------------------------- } : FABS ( f -- f') UNPACK SWAP ABS SWAP PACK ; \ absolute value of floating point number : FNEGATE ( f -- f') UNPACK SWAP NEGATE SWAP PACK ; \ negate floating point number { -------------------------------------------------------------------- } : F+ ( f1 f2 -- f3) FDUP F0= IF FDROP EXIT THEN \ special case: f2=0 FSWAP FDUP F0= IF FDROP EXIT THEN \ special case: f1=0 UNPACK >R SWAP UNPACK R> \ s1 s2 e2 e1 2DUP MAX fpexp ! \ store greater exponent - \ s1 s2 (e2-e1) DUP 0> \ e2 greater exponent? IF >R SWAP R> THEN \ yes: s2 s1 (e2-e1) \ no: s1 s2 (e2-e1) ABS \ number of shifts DUP 5 > \ shifts would exceed significand's range? IF 2DROP fpexp @ \ return fp number that has greater exponent ELSE 0 ?DO 10 / LOOP \ perform up to five right shifts + \ add aligned significands S>D fpexp @ NORMALIZE THEN PACK ; \ add floating point numbers \ adjust significands according to greater exponent, add, \ and combine resulting significand with greater exponent : F- ( f1 f2 -- f3) FNEGATE F+ ; \ subtract floating point numbers { -------------------------------------------------------------------- } : F* ( f1 f2 -- f3) UNPACK >R SWAP UNPACK >R \ S: s2 s1 R: e1 e2 M* \ multiply significands 100000 M/ \ scale product (10^5) S>D R> R> + \ add exponents NORMALIZE PACK ; \ multiply floating point numbers \ multiply significands, add exponents : F/ ( f1 f2 -- f3) UNPACK >R SWAP UNPACK >R \ S: s2 s1 R: e1 e2 1000000 M* \ inflate numerator (10^6) ROT M/ \ divide significands S>D R> R> - 1- \ subtract exponents and adjust NORMALIZE PACK ; \ divide floating point numbers \ divide significands, subtract exponents { -------------------------------------------------------------------- } : F10* ( f -- f') FDUP F0<> IF UNPACK 1+ PACK THEN ; \ increase exponent by one (unless f=0) : F10/ ( f -- f') FDUP F0<> IF UNPACK 1- PACK THEN ; \ decrease exponent by one (unless f=0) { -------------------------------------------------------------------- Comparison /2/ -------------------------------------------------------------------- } : F< ( f1 f2 -- ?) F- F0< ; : F> ( f1 f2 -- ?) F- F0> ; \ compare floating point numbers : FMAX ( f1 f2 -- f3) F2DUP F< IF FSWAP THEN FDROP ; \ return greater number (more positive number) : FMIN ( f1 f2 -- f3) F2DUP F> IF FSWAP THEN FDROP ; \ return smaller number (more negative number) { -------------------------------------------------------------------- Square root -------------------------------------------------------------------- } PRIVATE 2VARIABLE arg \ for storage of the argument value : SQRT ( ud -- u) arg 2! 999999 \ x_0 starting value BEGIN \ x DUP \ x x arg 2@ ROT \ x a. x UM/MOD NIP \ x (a/x) 2DUP U> \ values not yet final? WHILE + 2/ REPEAT \ x_new D>S ; \ Newton iteration: x_new = (x + a/x) / 2 for a <= 999999² \ starting with x_0, x decreases steadily and a/x increases \ movement stops when x = a/x, i.e. x = sqrt(a) PUBLIC : FSQRT ( f -- f') FDUP F0= IF EXIT THEN UNPACK S>D 2 FM/MOD \ divide exponent by 2 (floored!) >R \ save new exponent IF 1000000 \ mod = 1, exponent was odd: 10^6 ELSE 100000 THEN \ mod = 0, exponent was even: 10^5 UM* \ inflate significand SQRT \ calc new significand (already normalized) R> \ return new exponent PACK ; \ square root of floating point number \ scaling factors (10^5 or 10^6) ensure that SQRT returns values \ in range 100000 ... 999999, i.e. normalized significands { -------------------------------------------------------------------- Decadic logarithm (by table interpolation) -------------------------------------------------------------------- } PRIVATE CREATE logtable ( 1.0) 0000000 , 0010723 , 0021189 , 0031408 , ( 1.1) 0041392 , 0051152 , 0060697 , 0070037 , ( 1.2) 0079181 , 0088136 , 0096910 , 0105510 , ( 1.3) 0113943 , 0122215 , 0130333 , 0138302 , ( 1.4) 0146128 , 0153814 , 0161368 , 0168792 , ( 1.5) 0176091 , 0183269 , 0190331 , 0197280 , ( 1.6) 0204119 , 0210853 , 0217483 , 0224014 , ( 1.7) 0230448 , 0236789 , 0243038 , 0249198 , ( 1.8) 0255272 , 0261262 , 0267171 , 0273001 , ( 1.9) 0278753 , 0284430 , 0290034 , 0295567 , ( 2.0) 0301029 , 0306425 , 0311753 , 0317018 , ( 2.1) 0322219 , 0327358 , 0332438 , 0337459 , ( 2.2) 0342422 , 0347330 , 0352182 , 0356981 , ( 2.3) 0361727 , 0366422 , 0371067 , 0375663 , ( 2.4) 0380211 , 0384711 , 0389166 , 0393575 , ( 2.5) 0397940 , 0402261 , 0406540 , 0410777 , ( 2.6) 0414973 , 0419129 , 0423245 , 0427323 , ( 2.7) 0431363 , 0435366 , 0439332 , 0443262 , ( 2.8) 0447158 , 0451018 , 0454844 , 0458637 , ( 2.9) 0462397 , 0466125 , 0469822 , 0473486 , ( 3.0) 0477121 , 0480725 , 0484299 , 0487845 , ( 3.1) 0491361 , 0494850 , 0498310 , 0501743 , ( 3.2) 0505149 , 0508529 , 0511883 , 0515211 , ( 3.3) 0518513 , 0521791 , 0525044 , 0528273 , ( 3.4) 0531478 , 0534660 , 0537819 , 0540954 , ( 3.5) 0544068 , 0547159 , 0550228 , 0553276 , ( 3.6) 0556302 , 0559308 , 0562292 , 0565257 , ( 3.7) 0568201 , 0571126 , 0574031 , 0576916 , ( 3.8) 0579783 , 0582631 , 0585460 , 0588271 , ( 3.9) 0591064 , 0593839 , 0596597 , 0599337 , ( 4.0) 0602059 , 0604765 , 0607455 , 0610127 , ( 4.1) 0612783 , 0615423 , 0618048 , 0620656 , ( 4.2) 0623249 , 0625826 , 0628388 , 0630936 , ( 4.3) 0633468 , 0635986 , 0638489 , 0640978 , ( 4.4) 0643452 , 0645913 , 0648360 , 0650793 , ( 4.5) 0653212 , 0655618 , 0658011 , 0660391 , ( 4.6) 0662757 , 0665111 , 0667452 , 0669781 , ( 4.7) 0672097 , 0674401 , 0676693 , 0678973 , ( 4.8) 0681241 , 0683497 , 0685741 , 0687974 , ( 4.9) 0690196 , 0692406 , 0694605 , 0696793 , ( 5.0) 0698970 , 0701136 , 0703291 , 0705436 , ( 5.1) 0707570 , 0709693 , 0711807 , 0713910 , ( 5.2) 0716003 , 0718086 , 0720159 , 0722222 , ( 5.3) 0724275 , 0726319 , 0728353 , 0730378 , ( 5.4) 0732393 , 0734399 , 0736396 , 0738384 , ( 5.5) 0740362 , 0742332 , 0744292 , 0746244 , ( 5.6) 0748188 , 0750122 , 0752048 , 0753965 , ( 5.7) 0755874 , 0757775 , 0759667 , 0761551 , ( 5.8) 0763427 , 0765295 , 0767155 , 0769007 , ( 5.9) 0770852 , 0772688 , 0774516 , 0776337 , ( 6.0) 0778151 , 0779957 , 0781755 , 0783546 , ( 6.1) 0785329 , 0787106 , 0788875 , 0790636 , ( 6.2) 0792391 , 0794139 , 0795880 , 0797613 , ( 6.3) 0799340 , 0801060 , 0802773 , 0804480 , ( 6.4) 0806179 , 0807873 , 0809559 , 0811239 , ( 6.5) 0812913 , 0814580 , 0816241 , 0817895 , ( 6.6) 0819543 , 0821185 , 0822821 , 0824451 , ( 6.7) 0826074 , 0827692 , 0829303 , 0830909 , ( 6.8) 0832508 , 0834102 , 0835690 , 0837272 , ( 6.9) 0838849 , 0840419 , 0841984 , 0843544 , ( 7.0) 0845098 , 0846646 , 0848189 , 0849726 , ( 7.1) 0851258 , 0852784 , 0854306 , 0855821 , ( 7.2) 0857332 , 0858837 , 0860338 , 0861832 , ( 7.3) 0863322 , 0864807 , 0866287 , 0867762 , ( 7.4) 0869231 , 0870696 , 0872156 , 0873611 , ( 7.5) 0875061 , 0876506 , 0877946 , 0879382 , ( 7.6) 0880813 , 0882239 , 0883661 , 0885078 , ( 7.7) 0886490 , 0887898 , 0889301 , 0890700 , ( 7.8) 0892094 , 0893484 , 0894869 , 0896250 , ( 7.9) 0897627 , 0898999 , 0900367 , 0901730 , ( 8.0) 0903089 , 0904445 , 0905795 , 0907142 , ( 8.1) 0908485 , 0909823 , 0911157 , 0912487 , ( 8.2) 0913813 , 0915135 , 0916453 , 0917768 , ( 8.3) 0919078 , 0920384 , 0921686 , 0922984 , ( 8.4) 0924279 , 0925569 , 0926856 , 0928139 , ( 8.5) 0929418 , 0930694 , 0931966 , 0933234 , ( 8.6) 0934498 , 0935759 , 0937016 , 0938269 , ( 8.7) 0939519 , 0940765 , 0942008 , 0943247 , ( 8.8) 0944482 , 0945714 , 0946943 , 0948168 , ( 8.9) 0949390 , 0950608 , 0951823 , 0953034 , ( 9.0) 0954242 , 0955447 , 0956648 , 0957846 , ( 9.1) 0959041 , 0960232 , 0961421 , 0962606 , ( 9.2) 0963787 , 0964966 , 0966141 , 0967313 , ( 9.3) 0968482 , 0969648 , 0970811 , 0971971 , ( 9.4) 0973127 , 0974281 , 0975431 , 0976579 , ( 9.5) 0977723 , 0978864 , 0980003 , 0981138 , ( 9.6) 0982271 , 0983400 , 0984527 , 0985650 , ( 9.7) 0986771 , 0987889 , 0989004 , 0990116 , ( 9.8) 0991226 , 0992332 , 0993436 , 0994537 , ( 9.9) 0995635 , 0996730 , 0997823 , 0998912 , ( 10.0) 1000000 , : LOG ( ud -- u) 2500 UM/MOD \ mod quot 40 - \ offset to start of table CELLS logtable + \ addr of left neighbor 2@ \ mod y2 y1 DUP >R - \ mod y2-y1 2500 */ \ dy R> + ; \ y1+dy \ interpolate logtable for 1.00000 <= ud <= 9.99999 \ table is given in steps of (x2-x1) = 0.02500 = 2500 \ offset is necessary because table starts at 1.0 (not 0.0) \ interpolation: if x = x1+dx, y = y1+dy \ then dy = dx * (y2-y1) / (x2-x1) where dx = mod PUBLIC : FLOG ( f -- f') UNPACK \ s e 1000000 * \ s e000000; 6 digits behind decimal point SWAP \ e000000 s S>D LOG \ e000000 0ssssss + S>D \ essssss; new significand -1 \ new exponent; back to 5 digits behind decimal point NORMALIZE PACK ; \ decadic logarithm of floating point number \ log(s*10^e) = e + log(s) { -------------------------------------------------------------------- Decadic exponential function (by table interpolation) -------------------------------------------------------------------- } PRIVATE CREATE alogtable ( 0.00) 100000 , 100577 , 101157 , 101741 , ( 0.01) 102329 , 102920 , 103514 , 104111 , ( 0.02) 104712 , 105317 , 105925 , 106536 , ( 0.03) 107151 , 107770 , 108392 , 109018 , ( 0.04) 109647 , 110280 , 110917 , 111557 , ( 0.05) 112201 , 112849 , 113501 , 114156 , ( 0.06) 114815 , 115478 , 116144 , 116815 , ( 0.07) 117489 , 118168 , 118850 , 119536 , ( 0.08) 120226 , 120920 , 121618 , 122320 , ( 0.09) 123026 , 123737 , 124451 , 125169 , ( 0.10) 125892 , 126619 , 127350 , 128085 , ( 0.11) 128824 , 129568 , 130316 , 131069 , ( 0.12) 131825 , 132586 , 133352 , 134121 , ( 0.13) 134896 , 135675 , 136458 , 137246 , ( 0.14) 138038 , 138835 , 139636 , 140442 , ( 0.15) 141253 , 142069 , 142889 , 143714 , ( 0.16) 144543 , 145378 , 146217 , 147061 , ( 0.17) 147910 , 148764 , 149623 , 150487 , ( 0.18) 151356 , 152229 , 153108 , 153992 , ( 0.19) 154881 , 155775 , 156675 , 157579 , ( 0.20) 158489 , 159404 , 160324 , 161250 , ( 0.21) 162181 , 163117 , 164058 , 165006 , ( 0.22) 165958 , 166916 , 167880 , 168849 , ( 0.23) 169824 , 170804 , 171790 , 172782 , ( 0.24) 173780 , 174783 , 175792 , 176807 , ( 0.25) 177827 , 178854 , 179887 , 180925 , ( 0.26) 181970 , 183020 , 184077 , 185139 , ( 0.27) 186208 , 187283 , 188364 , 189452 , ( 0.28) 190546 , 191646 , 192752 , 193865 , ( 0.29) 194984 , 196110 , 197242 , 198380 , ( 0.30) 199526 , 200678 , 201836 , 203001 , ( 0.31) 204173 , 205352 , 206538 , 207730 , ( 0.32) 208929 , 210135 , 211348 , 212569 , ( 0.33) 213796 , 215030 , 216271 , 217520 , ( 0.34) 218776 , 220039 , 221309 , 222587 , ( 0.35) 223872 , 225164 , 226464 , 227771 , ( 0.36) 229086 , 230409 , 231739 , 233077 , ( 0.37) 234422 , 235776 , 237137 , 238506 , ( 0.38) 239883 , 241268 , 242661 , 244061 , ( 0.39) 245470 , 246888 , 248313 , 249746 , ( 0.40) 251188 , 252638 , 254097 , 255564 , ( 0.41) 257039 , 258523 , 260015 , 261517 , ( 0.42) 263026 , 264545 , 266072 , 267608 , ( 0.43) 269153 , 270707 , 272270 , 273841 , ( 0.44) 275422 , 277012 , 278612 , 280220 , ( 0.45) 281838 , 283465 , 285101 , 286747 , ( 0.46) 288403 , 290068 , 291742 , 293426 , ( 0.47) 295120 , 296824 , 298538 , 300261 , ( 0.48) 301995 , 303738 , 305492 , 307255 , ( 0.49) 309029 , 310813 , 312607 , 314412 , ( 0.50) 316227 , 318053 , 319889 , 321736 , ( 0.51) 323593 , 325461 , 327340 , 329230 , ( 0.52) 331131 , 333042 , 334965 , 336899 , ( 0.53) 338844 , 340800 , 342767 , 344746 , ( 0.54) 346736 , 348738 , 350751 , 352776 , ( 0.55) 354813 , 356861 , 358921 , 360994 , ( 0.56) 363078 , 365174 , 367282 , 369402 , ( 0.57) 371535 , 373680 , 375837 , 378007 , ( 0.58) 380189 , 382384 , 384591 , 386812 , ( 0.59) 389045 , 391291 , 393550 , 395822 , ( 0.60) 398107 , 400405 , 402717 , 405041 , ( 0.61) 407380 , 409732 , 412097 , 414476 , ( 0.62) 416869 , 419275 , 421696 , 424130 , ( 0.63) 426579 , 429042 , 431519 , 434010 , ( 0.64) 436515 , 439035 , 441570 , 444119 , ( 0.65) 446683 , 449262 , 451855 , 454464 , ( 0.66) 457088 , 459726 , 462381 , 465050 , ( 0.67) 467735 , 470435 , 473151 , 475882 , ( 0.68) 478630 , 481393 , 484172 , 486967 , ( 0.69) 489778 , 492606 , 495450 , 498310 , ( 0.70) 501187 , 504080 , 506990 , 509917 , ( 0.71) 512861 , 515822 , 518800 , 521795 , ( 0.72) 524807 , 527837 , 530884 , 533949 , ( 0.73) 537031 , 540132 , 543250 , 546386 , ( 0.74) 549540 , 552713 , 555904 , 559113 , ( 0.75) 562341 , 565587 , 568852 , 572136 , ( 0.76) 575439 , 578761 , 582103 , 585463 , ( 0.77) 588843 , 592243 , 595662 , 599100 , ( 0.78) 602559 , 606038 , 609536 , 613055 , ( 0.79) 616595 , 620154 , 623734 , 627335 , ( 0.80) 630957 , 634599 , 638263 , 641948 , ( 0.81) 645654 , 649381 , 653130 , 656901 , ( 0.82) 660693 , 664507 , 668343 , 672202 , ( 0.83) 676082 , 679986 , 683911 , 687859 , ( 0.84) 691830 , 695824 , 699841 , 703882 , ( 0.85) 707945 , 712032 , 716143 , 720277 , ( 0.86) 724435 , 728618 , 732824 , 737055 , ( 0.87) 741310 , 745589 , 749894 , 754223 , ( 0.88) 758577 , 762956 , 767361 , 771791 , ( 0.89) 776247 , 780728 , 785235 , 789768 , ( 0.90) 794328 , 798913 , 803526 , 808164 , ( 0.91) 812830 , 817523 , 822242 , 826989 , ( 0.92) 831763 , 836565 , 841395 , 846252 , ( 0.93) 851138 , 856051 , 860993 , 865964 , ( 0.94) 870963 , 875991 , 881048 , 886135 , ( 0.95) 891250 , 896396 , 901571 , 906775 , ( 0.96) 912010 , 917275 , 922571 , 927897 , ( 0.97) 933254 , 938642 , 944060 , 949510 , ( 0.98) 954992 , 960505 , 966050 , 971627 , ( 0.99) 977237 , 982878 , 988553 , 994260 , ( 1.00) 1000000 , : ALOG ( u -- u') 250 /MOD \ mod quot CELLS alogtable + \ addr of left neighbor 2@ \ mod y2 y1 DUP >R - \ mod y2-y1 250 */ \ dy R> + ; \ y1+dy \ interpolate alogtable for 0.00000 <= ud <= 0.99999 \ table is given in steps of (x2-x1) = 0.00250 = 250 \ interpolation: if x = x1+dx, y = y1+dy \ then dy = dx * (y2-y1) / (x2-x1) where dx = mod PUBLIC : FALOG ( f -- f') UNPACK \ s e SWAP S>D ROT \ d-sig n-exp DSHIFT \ fd (fixed point double length number) 100000 FM/MOD \ mod quot SWAP ALOG \ calc new significand (already normalized) SWAP \ return new exponent PACK ; \ decadic exponential function of floating point number \ 10^(a.bbbbb) = 10^(a+0.bbbbb) = 10^(0.bbbbb)*10^a \ 100000 FM/MOD SWAP splits up the fixed point number \ 6.12345 is split into 6 and 0.12345 \ -6.12345 is split into -7 and 0.87655 \ floored division is required for negative f numbers { ==================================================================== } END-PACKAGE