\ 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. : <> { 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" *** <> 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" *** <> 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 ; \ <> 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. : { 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 <> node .den ! node .num 2! ; : simplify ( head -- ) ['] 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@ r> ; : new-ratio ( Dnum den -- node ) ratio alloc init-ratio ; : ( node -- ) dealloc ; : kill-ratio ( head -- ) each[ ]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 -- ) \ like U. except with no space afterward u>d <# #s #> type ; : { 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 dup if node .num 2@ d0< if [char] - else [char] + then emit then then ?dup if [char] / emit node .den @ u. else bl emit then ; : show-ratio ( head -- ) ." | " ['] each ." | " ; : { node -- } node .num 2@ dnegate node .num 2! ; : ratio-negate ( head -- ) \ negate all of the nodes in head-list ['] each ; : { node | sgn -- } node .num 2@ dup 0< to sgn dabs abort" *** 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 ['] 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. : { 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 ; : ratio-sqrt ( head -- ) \ square-root all of the nodes in head-list ['] 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 ; : { 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 <> node .den ! node .num 2! ; : ratio-sq ( head -- ) \ square all of the nodes in head-list ['] each ; : { 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 dest \ this helps to prevent overflow dest .num 2@ surr .num 2@ d* \ -- Dnum dest .den @ surr .den @ um* \ -- Dnum Dden <> dest .den ! dest .num 2! surr dealloc dest ; : product-ratio ( product-node head -- ) \ multiply all of head-list into product-node ['] each drop ; : <*ratio> ( multiplier node -- multiplier ) over 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 ! ; : { 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 <> dest .den ! dest .num 2! dest ; : sum-ratio ( sum-node head -- ) \ add all of head-list into the sum-node ['] each drop ; : <+ratio> ( addend node -- addend ) over 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 | AB 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 AB then exit then A=B exit then \ same signs nodeA .num 2@ nip 0< if AB then exit then nodeB .den @ 0= if nodeB .num 2@ nip 0< if A>B else AB then ; : ( new-node node -- new-node ? ) over ratio-compare A>B = ; : ratio-sort ( head -- new-head ) ['] 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 ; : { head | neg-mean surr result -- result-node } \ the variance numerator head mean dup 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 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 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 ;