\ ratio.4th --- written by Hugh Aguilar --- copyright (c) 2010, BSD license
\ Compiles under any ANS-Forth system (16-bit, 32-bit or 64-bit).
\ This is an implementation of lists of rational numbers
\ See the file RATIO.TXT for further documentation.
\ Author contact: hughaguilar96@yahoo.com
\ We need NOVICE.4TH and LIST.4TH already compiled.
marker ratio.4th
list
d field .num \ signed numerator the sign of the numerator is the sign of the ratio
w field .den \ unsigned denominator
constant ratio
\ The denominator can be zero, which indicates infinity; the sign of the numerator is still valid.
\ The numerator is a double to help support ratios larger than one.
: <<simplify>> { Lnum Hnum Lden Hden | sgn whole divisor -- Dnum den }
Hnum 0< to sgn
Lnum Hnum dabs to Hnum to Lnum
Lnum Hnum Lden Hden dgcd abort" *** <<SIMPLIFY>> overflow in DGCD ***"
begin dup even? while \ halve everything until divisor is odd
Lnum Hnum d2/ to Hnum to Lnum
Lden Hden d2/ to Hden to Lden
2/ repeat to divisor
divisor 1 = if \ if divisor is 1, we are done
Hden abort" *** <<SIMPLIFY>> overflow in denominator***"
Lnum Hnum sgn if dnegate then
Lden exit then
Hden 0= if \ if denominator is a single, reduce numerator
Lnum Hnum Lden um/mod to whole \ numerator is now single-precision
to Lnum 0 to Hnum then
Lden Hden divisor um/mod nip to Lden \ Hden is no longer valid
Lnum Hnum divisor um/mod nip u>d
whole Lden um* d+ sgn if dnegate then
Lden ;
\ <<SIMPLIFY>> first halves everything until the divisor is odd.
\ This might convert the denominator from a double to a single, which will help.
\ If the divisor is 1, then we are done; we can't simplify any further.
\ If the denominator is a single, the numerator is converted into a single.
\ This prevents overflow in the case that the divisor is small.
\ Finally both the numerator and denominator are divided by the divisor.
: <simplify> { node -- }
node .num 2@ d0= if exit then \ ignore zero
node .den @ 0= if exit then \ ignore infinite
node .num 2@ node .den @ u>d <<simplify>> node .den ! node .num 2! ;
: simplify ( head -- )
['] <simplify> each ;
: init-ratio ( Dnum den node -- node )
init-list >r
dup 0< abort" *** INIT-RATIO doesn't allow negative denominators ***"
r@ .den !
r@ .num 2!
r@ <simplify>
r> ;
: new-ratio ( Dnum den -- node )
ratio alloc
init-ratio ;
: <kill-ratio> ( node -- )
dealloc ;
: kill-ratio ( head -- )
each[ <kill-ratio> ]each ;
: add-int { n node | sgn -- }
n 0< to sgn
node .den @ n um* sgn if dnegate then
node .num 2@ d+ node .num 2! ;
min-signed constant | \ the bad number used as a bracket
\ The MIN-SIGNED value hopefully won't be used by accident, but this would cause a bug.
: make-ratio ( bad n... bad -- head ) \ fills numerators, sets denominators to 1
>r \ r: -- bad
nil begin over r@ <> while \ -- bad n... head
swap s>d 1 new-ratio swap link repeat
swap r> 2drop ; \ throw away the two bads
\ MAKE-RATIO only accepts single-precision numerators.
\ Our .NUM field is double-precision in order to support ratios greater than one.
: set-dens ( bad n... bad head -- new-head ) \ sets denominators
reverse dup >r
swap >r \ r: -- head bad
begin over r@ <> while
dup 0= abort" *** SET-DENS had too many parameters ***"
over 0< abort" *** SET-DENS doesn't allow negative denominators ***"
tuck .den ! .fore @ repeat drop \ throw away the head (NIL by now)
swap r> 2drop \ throw away the two bads
r> reverse
dup simplify ;
\ SET-DENS disallows negative numbers because this is most likely a mistake.
\ It is unlikely that the user was entering a gigantic unsigned number.
: add-ints ( bad n... bad head -- new-head ) \ adds whole numbers to numerators
reverse dup >r
swap >r \ r: -- head bad
begin over r@ <> while
dup 0= abort" *** ADD-INTS had too many parameters ***"
tuck add-int .fore @ repeat drop \ throw away the head (NIL by now)
swap r> 2drop \ throw away the two bads
r> reverse
dup simplify ;
: <u.> ( u -- ) \ like U. except with no space afterward
u>d <# #s #> type ;
: <show-ratio> { node -- }
node .den @ 0= if node .num @ 0< if ." -inf " else ." +inf " then exit then
node .num 2@ d0= if 0 . exit then
node .num 2@ d0< if [char] - emit then
node .num 2@ dabs node .den @ um/mod
?dup if <u.>
dup if
node .num 2@ d0< if [char] - else [char] + then emit
then
then
?dup if <u.> [char] / emit node .den @ u.
else bl emit then ;
: show-ratio ( head -- )
." | "
['] <show-ratio> each
." | " ;
: <ratio-negate> { node -- }
node .num 2@ dnegate node .num 2! ;
: ratio-negate ( head -- ) \ negate all of the nodes in head-list
['] <ratio-negate> each ;
: <ratio-invert> { node | sgn -- }
node .num 2@ dup 0< to sgn dabs
abort" *** <RATIO-INVERT> overflow ***" \ -- new-den
node .den @ u>d sgn if dnegate then
node .num 2!
node .den ! ;
: ratio-invert ( head -- ) \ invert all of the nodes in head-list
['] <ratio-invert> each ;
100000000. 2constant ratio-scalar
\ RATIO-SCALAR is small so that the ratio can later on be multiplied by another ratio
\ and it hopefully won't overflow. Overflow is the big problem with these ratios.
: <ratio-sqrt> { node -- }
node .num 2@ d0< abort" *** RATIO-SQRT can't accept negative arguments ***"
node .den @ 0= if exit then
ratio-scalar node .den @ um/mod nip dup um*
2dup node .num 2@ d* dsqrt u>d node .num 2!
node .den @ du* dsqrt node .den !
node <simplify> ;
: ratio-sqrt ( head -- ) \ square-root all of the nodes in head-list
['] <ratio-sqrt> each ;
: ratio>float { node -- float }
node .num 2@ d>f
node .den @ dup 0= abort" *** RATIO>FLOAT divide by zero ***"
s>d d>f f/ ;
: float>ratio ( float -- node )
ratio-scalar d>f f* f>d
ratio-scalar drop new-ratio ;
: <ratio-sq> { node -- }
node .den @ 0= if 1. node .num 2! exit then \ infinites become positive when squared
node .num 2@ 2dup d*
node .den @ dup um* \ -- Dnum Dden
<<simplify>> node .den ! node .num 2! ;
: ratio-sq ( head -- ) \ square all of the nodes in head-list
['] <ratio-sq> each ;
: <multiply-ratio> { dest node | surr -- dest }
node .den @ 0= dest .den @ 0= or if 0 dest .den !
dest .num 2@ nip node .num 2@ nip xor 0< if -1. else 1. then
dest .num 2! dest exit then
node clone-node to surr \ use a clone of NODE so we can modify it
surr .num 2@ dest .num 2@
surr .num 2! dest .num 2! \ exchange the numerators
surr <simplify> dest <simplify> \ this helps to prevent overflow
dest .num 2@ surr .num 2@ d* \ -- Dnum
dest .den @ surr .den @ um* \ -- Dnum Dden
<<simplify>> dest .den ! dest .num 2!
surr dealloc dest ;
: product-ratio ( product-node head -- ) \ multiply all of head-list into product-node
['] <multiply-ratio> each drop ;
: <*ratio> ( multiplier node -- multiplier )
over <multiply-ratio> drop ;
: *ratio ( multiplier head -- ) \ multiply multiplier into all of the nodes in head-list
['] <*ratio> each drop ;
: zero-ratio { node -- } \ makes node into a zero
0. node .num 2!
1 node .den ! ;
: <add-ratio> { dest node | dest-mult -- dest } \ adds node to dest
node .den @ 0= if
dest .den @ 0= if
node .num 2@ nip dest .num 2@ nip xor 0< if dest zero-ratio then
exit then
node .num @ dest .num !
node .den @ dest .den ! exit then
dest .den @ 0= if exit then
node .den @ dest .den @ gcd \ -- divisor
node .den @ over / to dest-mult
dest .den @ swap / \ -- Unode-mult
node .num 2@ rot du* \ -- Dnode-num
dest .num 2@ dest-mult du* \ -- Dnode-num Ddest-num
d+ \ -- Dnum
dest .den @ dest-mult um* \ -- Dnum Dden
<<simplify>> dest .den ! dest .num 2!
dest ;
: sum-ratio ( sum-node head -- ) \ add all of head-list into the sum-node
['] <add-ratio> each drop ;
: <+ratio> ( addend node -- addend )
over <add-ratio> drop ;
: +ratio ( addend head -- ) \ add addend to all of the nodes in head-list
['] <+ratio> each drop ;
: ratio-compare { nodeA nodeB -- A>B | A=B | A<B }
nodeA 0= if
nodeB 0= if A=B exit then
A<B exit then
nodeB 0= if A>B exit then
nodeA .den @ 0= if
nodeB .den @ 0= if
nodeA .num 2@ nip nodeB .num 2@ nip xor 0< if \ different signs
nodeA .num 2@ nip 0< if A<B else A>B then
exit then
A=B exit then \ same signs
nodeA .num 2@ nip 0< if A<B else A>B then
exit then
nodeB .den @ 0= if
nodeB .num 2@ nip 0< if A>B else A<B then
exit then
nodeA .den @ nodeB .den @ gcd \ -- divisor
nodeA .den @ over / \ -- divisor multB
nodeB .den @ rot / \ -- multB multA
nodeA .num 2@ rot du*
rot nodeB .num 2@ rot du* \ -- DnumA DnumB
2over 2over d= if 2drop 2drop A=B exit then
d< if A<B else A>B then ;
: <ratio-sort> ( new-node node -- new-node ? )
over ratio-compare A>B = ;
: ratio-sort ( head -- new-head )
['] <ratio-sort> sort-list ;
: mean { head | n result divisor -- result-node } \ the average of the head-list
head length to n
0. 1 new-ratio to result
1. n new-ratio to divisor
result head sum-ratio
result divisor product-ratio
divisor kill-ratio
result ;
: <variance> { head | neg-mean surr result -- result-node } \ the variance numerator
head mean dup <ratio-negate> to neg-mean
0. 1 new-ratio to result
head clone-list to surr
neg-mean surr +ratio surr ratio-sq
result surr sum-ratio
neg-mean kill-ratio surr kill-ratio
result ;
: sample-variance { head | mult -- result-node } \ variance of a sample
1. over length 1- new-ratio to mult
head <variance> dup mult product-ratio \ -- result-node
mult kill-ratio ;
: population-variance { head | mult -- result-node } \ variance of a complete population
1. over length new-ratio to mult
head <variance> dup mult product-ratio \ -- result-node
mult kill-ratio ;
: test ( -- head )
| 200 250 125 | make-ratio >r
| 15 75 55 | r> set-dens >r
| 2 0 3 | r> add-ints ;