Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- \ 32 bit integer, fixed and floating point arithmetic.
- comment:
- The 32 bit integer and fixed point maths is by:-
- Prof Tim Hendtlass, Physics Department, Swinburne University of Technology,
- P.O.BOX 218 Hawthorn 3122 Australia. Voice 61 3 819 8863 Fax 61 3 818 3645.
- Please let me know of any bugs, fixes, add ons, add ins etc. The 32 bit
- floating point routines are slightly reworked from the originals by Martin
- Tracy. The cute little transcindental package is by Nathaniel Grossman.
- This code originally appeared in Forth Dimensions in an article entitled
- "Maths - Who needs it?" As well as some explanation, the article contains
- timing data and space requirments.
- My code is hereby placed in the public domain, I believe all the other code
- in this file to have also been placed in the public doimain by the authors.
- Happy Forthing. Tim. 18/8/92
- comment;
- \ *********************************************************
- \ * 32 bit Integer Arithmetic *
- \ *********************************************************
- \ Don't forget that the normal D+ and D- and DABS work for 32 bit integer
- \ and fixed point numbers too!
- \ *********************************
- \ * 32 bit Integer Multiplication *
- \ *********************************
- \ Unsigned double * unsigned double = unsigned double (No overflow check)
- : UD* \ unsigned 32 bit answer, no overflow check
- rot >r over >r >r over >r \ put c b a d on return stack
- um* \ b*d = part of 32 bit answer
- 2r> * 2r> * + + \ a*d+b*c= addition to top 16 bits
- ;
- \ Signed double * signed double = signed double (No overflow check)
- : D* \ signed, no overflow check
- dup>r dabs 2swap dup>r dabs \ #s +ve, keep info to work out final sign
- ud* \ get 32 bit answer
- 2r> xor ?dnegate \ work out and apply final sign
- ;
- \ Unsigned double * unsigned double = unsigned double (with overflow check)
- : UD*C \ unsigned, with overflow check
- dup>r rot dup>r >r over >r \ put a c c b on return stack
- >r swap dup>r \ put a d onto return stack
- um* \ b*d
- 0 2r> um* d+ 2r> um* d+ \ offset 16 bits, add on a*d+b*c
- 0 2r> um* d+ \ off another 16 bits, add on a*c
- or 0<> abort" D* overflow" \ check for overflow
- ;
- \ ***************************
- \ * 32 bit Integer Division *
- \ ***************************
- \ Use fast algorithm, remainder requires an additional 32 bit multiplication
- \ and sutraction.
- \ Unsigned double * unsigned single = unsigned triple
- : T* ( ud un -- ut )
- dup \ ud un un
- rot \ udl un un udh
- um* \ udl un high-ans
- 2>r \ udl un
- um* 0 2r> d+ \ low-ans then add on high-answer after ofsetting it 16 bits
- ;
- \ Unsigned triple / unsigned single = unsigned double
- : T/ ( ut un -- ud )
- >r r@ um/mod swap \ divisor to r, divide top two words, rem to top
- rot 0 r@ um/mod swap \
- rot r> um/mod swap drop
- 0 2swap swap d+
- ;
- \ Calculate ud * un1 / un2. Triple intermediate product.
- : U*/ ( ud un1 un2 -- ud2 )
- >r t* r> t/
- ;
- \ Unsigned 32 bit by 32 bit divide. No remainder.
- : UD/ ( ud1 ud2 -- ud3 )
- dup 0= \ top 16 bits of divisor = 0?
- if swap t/ \ make it a triple, do the division
- else
- dup 65536. rot 1+ um/mod >r \ work out scaling factor,copy to return stack
- drop r@ t* drop 2>r \ scale denominator, move to return stack
- dup 0 2r@ u*/ d- \ calculate (U-U0*W1/W0)
- 2r> r> -rot nip u*/ \ multiply by (D/W0)
- nip 0 \ /2^16 (use top 16bits only),make ans double
- then
- ;
- \ Divides two double numbers. All numbers are signed doubles.
- : D/MOD ( dn1 dn2 -- drem dquot )
- 2 pick over xor >r \ work out sign of answer
- dabs 2swap dabs 2swap \ convert numbers to positive
- 4dup ud/ 2dup 2>r \ do the division, save copy ans
- ud* d- \ calculate remainder
- 2r> r> ?dnegate \ retrieve answer,apply final sign
- ;
- : D/ ( dn1 dn2 -- dquot )
- 2 pick over xor >r \ work out sign of answer
- dabs 2swap dabs 2swap \ convert numbers to positive
- ud/ \ do the division
- r> ?dnegate \ retrieve answer,apply final sign
- ;
- \ *********************************************************
- \ * 32 bit Fixed Point Arithmetic *
- \ *********************************************************
- \ *****************************************
- \ * Defining the fixed point structure *
- \ *****************************************
- variable fdpl variable fscl
- : FPLACES ( -- n ) fdpl @ ; \ number of implied decimal places
- : FSCALE ( -- n ) fscl @ ; \ scaling factor we are using
- : FIXED ( n -- )
- 0 max 4 min fdpl ! \ clip to between 0 and 4 decimal places
- 1 fplaces 0 ?do 10 * loop fscl ! \ store scaling factor
- ;
- 3 FIXED \ default to three decimal places
- \ *****************************************
- \ * Outputting numbers *
- \ *****************************************
- : (F.) ( fn -- adr len ) \ prepare fixed point # ready to output
- tuck \ keep copy of top byte so we know sign
- dabs \ convert to positive number
- <# bl hold \ start conversion with a leading blank
- fplaces 0 ?do # loop \ convert places after decimal point
- ascii . hold \ put a decimal point in place
- #s \ convert integer part
- rot sign #> \ put sign in place, tidy stack
- ;
- : FIX. ( fn -- ) (f.) type ; \ print fixed point number
- : FIX.R ( fn p -- ) \ right justify in a field of p places
- >r (f.) \ convert number
- r> over - 0 ?do bl emit loop type \ pad with blanks as needed
- ;
- \ *****************************************
- \ * Inputting numbers *
- \ *****************************************
- : D10* ( d1 -- 10*d1) \ multiply a 32 bit number by 10
- d2* 2dup d2* d2* d+ \ 8*d+2*d=10*d
- ;
- \ Convert number to fixed point number - no check made for numbers too large
- \ Example 1234.5 FIX. To compile a fixed point number in a : definition use
- \ [ 1234.5 FIX ] DLITERAL
- : FIX ( dn -- fn )
- dpl @ 0< \ single or double number?
- if s>d 0 dpl ! then \ if single convert to double
- dpl @ fplaces <> \ # decimal places entered not fplaces?
- if dpl @ fplaces < \ too few places specified?
- if fplaces dpl @ ?do d10* loop \ yes, too few so scale the number up
- else abort" Too mnay decimal places" \ no, too many - we can't handle this
- then
- then
- ;
- \ *****************************************
- \ * 32 bit Fixed Point Multiply *
- \ *****************************************
- \ Multiply two fixed point numbers producing a fixed point result.
- : FIX* ( f1 f2 -- f1*f2 )
- rot 2dup xor >r \ sign of answer to return stack
- -rot dabs 2swap dabs \ make both numbers positive
- dup>r rot dup>r >r over >r \ put a c c b on return stack
- >r swap dup>r \ put a d onto return stack
- um* \ b*d
- 0 2r> um* d+ 2r> um* d+ \ offset 16 bits, add on a*d+b*c
- 2r> * + \ add on low byte of a*c
- fscale mu/mod \ divide ms32 bits, ans to R.
- 0<> abort" Fixed * Overflow!" >r \ unless overflow quotient to R
- fscale mu/mod rot drop \ divide remainder and last 16 bits
- r> + r> ?dnegate \ assemble final answer, negate if required
- ;
- \ *****************************************
- \ * 32 bit Fixed Point Divide *
- \ *****************************************
- : FIX/ ( f1 f2 -- fquot=f1/f2 ) \ Divide two numbers
- 2 pick over xor >r \ work out sign of answer and save
- dabs 2swap dabs 2swap \ make all numbers positive
- 2dup >r >r \ keep copy of divisor
- d/mod fscale 0 d* \ scale integer part of answer
- 2swap fscale 0 d* \ and then scale remainder
- r> r> d/ \ divide remainder by divisor
- d+ \ add fract part of ans
- r> ?dnegate \ put on final sign
- ;
- : FIX>INT ( fn -- dn ) \ get 32 bit integer from a fixed point number
- fscale s>d d/ \ just add drop to get a 16 bit integer
- ;
- \ *********************************************************
- \ * 32 bit Floating Point Arithmetic *
- \ * Based on Zen Math by Martin Tracy *
- \ *********************************************************
- \ Trim a double-number mantissa and an exponent of ten to a floating number.
- : TRIM ( dn n = f)
- >r \ exponent to return stack
- tuck dabs \ save copy of high word for sign, make double positive
- begin over 0< over 0<> or \ MSB low word set or top 16 bits no zero?
- \ if so too big to fit into 16bits when signed
- while
- 0 10 um/mod >r 10 um/mod nip r> \ divide 32 bit mantissa by 10
- r> 1+ >r \ and increase exponent
- repeat rot ?dnegate drop r> \ apply sign and final exponent
- ;
- \ ****************************************************
- \ * 32 bit Floating Point Addition and Subtraction *
- \ ****************************************************
- : F+
- rot 2dup - dup 0< \ work out difference in exponents
- if \ top number has the larger exponent
- negate rot >r nip >r swap r> \ keep larger (on return stack) and diff, swap mantissas
- else \ top has a smaller or equal exponent
- swap >r nip \ keep larger (on return stack) and diff
- then
- >r s>d r> dup 0 \ convert mantissa to be shift to double
- ?do >r d10* r> 1- \ multiply mantissa by 10, decrement exponent
- over abs 6553 > \ would a *10 cause overflow of these 16 bits?
- if leave then \ prematurely terminate loop if so
- loop
- r> over + >r \ calculate final exponent
- if rot drop \
- else rot s>d d+
- then r> trim \ get final exponent and trim
- ;
- : FNEGATE >r negate r> ;
- : F- fnegate f+ ; \ add negative of the top value
- \ ****************************************
- \ * 32 bit Floating Point Multiplication *
- \ ****************************************
- : F* ( f1 f2 -- f3 )
- rot + >r \ calc exp of answer,save on return stack
- 2dup xor >r \ save xor of mantissas too (sign of answer)
- abs swap abs um* \ make mantissas positive and multiply
- r> ?dnegate r> trim \ apply sign and then get exponent and trim
- ;
- \ ***********************************
- \ * 32 bit Floating Point Division *
- \ ***********************************
- : F/
- over 0= abort" d/0 error!" \ check for divide by zero
- rot swap - >r \ get exponent of answer, put on return stack
- 2dup xor -rot \ get sign of answer, tuck down on stack
- abs dup 6553 min rot abs 0 \
- begin 2dup d10* nip 3 pick < \
- while d10* r> 1- >r \
- repeat 2swap drop um/mod \ now do the division
- nip 0 rot ?dnegate r> trim \ lose remainder, apply sign get exp and trim
- ;
- \ ******************************************
- \ * 32 bit Floating Point Input and Output *
- \ ******************************************
- \ Numbers to be floated must include a decimal point when entered.
- \ DPL contains the number of digits entered after the decimal point.
- : FLOAT ( n -- f) \ float the last entered number.
- dpl @ negate trim
- ;
- : F. ( f --) \ print a floating number in fixed format.
- >r dup abs 0
- <# r@ 0 max 0 ?do ascii 0 hold loop
- r@ 0<
- if r@ negate 0 max 0 ?do # loop ascii . hold
- then r> drop #s rot sign
- #> type space
- ;
- \ *********************************************************
- \ * 32 bit Floating Point Transcindental Functions *
- \ * Based on Zen sliderule by Nathaniel Grossman *
- \ *********************************************************
- comment:
- These functions are calculated to an accuracy of about one or two units in the
- third decimal place. Domains of the functions - scale input if needs be to
- keep within range. The CORDIC algorithm is used so these routines are not fast.
- SQRT 0.03 - 2.42
- LN 0.10 - 9.58
- SIN,COS,TAN -1.74 - 1.74
- ATAN -infinity - +infinity
- SINH,COSH,TANH -1.13 - 1.13
- ATANH -0.81 - 0.81
- comment;
- : F2* 2dup f+ ;
- : F/2 swap s>d 2dup d2* d2* d+ rot 1- trim ; \ mantissa * 5, dec exp by 1
- : F/2^N ( r n --- r/2^n) 0 ?do f/2 loop ; \ divide by 2, n times
- \ Convience renaming of existing words. Either these or original names work.
- ' 2DROP alias FDROP
- ' 2DUP alias FDUP
- ' 2OVER alias FOVER
- ' 2SWAP alias FSWAP
- ' 2ROT alias FROT
- ' 2! alias F!
- ' 2@ alias F@
- ' 2VARIABLE alias FVARIABLE
- \ Extra floating point words
- : F0< DROP 0< ;
- : F, FLOAT , , ;
- : FCONSTANT FLOAT 2CONSTANT ;
- \ Constants, Variables, Deferred Words and Arrays
- 4 CONSTANT F#BYTES VARIABLE MODE_FLAG
- 0.6073 FCONSTANT 1/K VARIABLE DELTA_FLAG
- 1.2076 FCONSTANT 1/K' VARIABLE NDX
- 0.0000 FCONSTANT F0 FVARIABLE F-BIN
- 1.0000 FCONSTANT F1 FVARIABLE XX
- 0.2500 FCONSTANT F1/4 FVARIABLE YY
- FVARIABLE ZZ
- DEFER EPS \ it will be either +EPS OR -EPS
- DEFER !STACK \ it will be R-STORE or V-STORE
- DEFER DO-IT \ it will be either ROT'ING or VEC'ING
- : FARRAY \ array building word
- CREATE \ ( -- ) compile time stack
- DOES> \ ( n --- adr ) run time stack
- swap f#bytes * + f@ \ calc address of nth entry
- ;
- FARRAY +EPS .7854 F, .4636 F, .2450 F, .1244 F, .06242 F,
- .03124 F, .05162 F, .007812 F, .003906 F, .001953 F,
- .0009766 F, .0004883 F, .0002441 F, .0001221 F, .00006104 F,
- FARRAY -EPS .0000 F, .5493 F, .2554 F, .1257 F, .06258 F,
- .03126 F, .01563 F, .007813 F, .003906 F, .001953 F,
- .0009766 F, .0004883 F, .0002441 F, .0001221 F, .00006104 F,
- : DELTA_SIGN delta_flag @ if fnegate then ;
- : MODE_SIGN mode_flag @ if fnegate then ;
- : R-DELTA= ( r --- r ) fdup f0< delta_flag ! ;
- : V-DELTA= ( r --- r ) fdup f0< not delta_flag ! ;
- : R-STORE ( x y z --- ) r-delta= zz f! yy f! xx f! ;
- : V-STORE ( x y z --- ) zz f! v-delta= yy f! xx f! ;
- : NEW_Z zz f@ ndx @ eps delta_sign f- ;
- : NEW_X xx f@ yy f@ ndx @ f/2^n delta_sign mode_sign f+ ;
- : NEW_Y yy f@ xx f@ ndx @ f/2^n delta_sign f+ ;
- : ROT'ING new_x new_y new_z r-store ;
- : VEC'ING new_x new_y new_z v-store ;
- : MODE=+1 -1 mode_flag ! ['] +eps is eps ;
- : MODE=-1 0 mode_flag ! ['] -eps is eps ;
- : CORDIC ( xstart ystart zstart --- xend yend zend )
- !stack mode_flag @ dup 0= >r
- if 0 ndx ! do-it then
- 4 1 do i ndx ! do-it loop
- r@ if 4 ndx ! do-it then
- 14 4 do i ndx ! do-it loop
- r> if 13 ndx ! do-it then
- 14 ndx ! do-it xx f@ yy f@ zz f@
- ;
- : FCOS&SIN ( r --- cos{r} sin{r} f )
- 1/k fswap f0 fswap ['] r-store is !stack
- mode=+1 ['] rot'ing is do-it cordic
- ;
- : FCOS ( r --- cos{r}) fcos&sin fdrop fdrop ;
- : FSIN ( r --- sin{r}) fcos&sin fdrop fswap fdrop ;
- : FTAN ( r --- tan{r}) fcos&sin fdrop fswap f/ ;
- : FCOSH&SINH ( r --- cosh{r} sinh{r} f )
- 1/k' fswap f0 fswap ['] r-store is !stack
- mode=-1 ['] rot'ing is do-it cordic
- ;
- : FCOSH ( r --- cosh{r}) fcosh&sinh fdrop fdrop ;
- : FSINH ( r --- sinh{r}) fcosh&sinh fdrop fswap fdrop ;
- : FTANH ( r --- tanh{r}) fcosh&sinh fdrop fswap f/ ;
- : FALN ( r --- exp{r}) fcosh&sinh fdrop f+ ;
- : FLN ( r --- ln{r} )
- fdup f1 f+ fswap f1 f- f0 ['] v-store is !stack
- mode=-1 ['] vec'ing is do-it cordic
- fswap fdrop fswap fdrop f2*
- ;
- : FSQRT ( r --- sqrt{r} )
- fdup f1/4 f+ fswap f1/4 f- f0 ['] v-store is !stack
- mode=-1 ['] vec'ing is do-it cordic
- fdrop fdrop 1/k' f*
- ;
- : R>P ( x y --- {x^2 + y^2}^1/2 arctan{y/x} )
- 1/k f* fswap 1/k f* fswap f0 ['] v-store is !stack
- mode=+1 ['] vec'ing is do-it cordic fswap fdrop
- ;
- : P>R ( radius angle --- x y )
- fover fswap fcos&sin fdrop frot f* frot frot f* fswap
- ;
- : FATAN ( r --- arctan{r}) f1 fswap r>p fswap fdrop ;
- : FATANH ( r --- argtanh{r}) f1 fswap f0 mode=-1
- ['] vec'ing is do-it cordic fswap fdrop fswap fdrop
- ;
Advertisement
Add Comment
Please, Sign In to add comment
