\ -*- mode: forth -*-
\ $Id: moments.gforth,v 1.1 2001/06/03 12:10:25 doug Exp $
\ http://www.bagley.org/~doug/shootout/
\ from Anton Ertl

1024 constant max-line
create line max-line 2 + allot

: input-floats 
    >r 0e begin
    line max-line r@ read-line throw
    while
    line swap >float 0= abort" float expected"
    fdup f, f+
    repeat
    rdrop drop ;

: compute-loop 
    dup 0 d>f fdup { f: n } f/ { f: mean }
    0e fdup fdup fdup
    floats bounds do {  f: avg-deviation f: variance f: skew f: kurtosis }
    i f@ mean f- { f: deviation }
    deviation fabs avg-deviation f+ 
    deviation fdup f* fdup variance f+ 
    fswap deviation f* fdup skew f+ 
    fswap deviation f* kurtosis f+ 
    float +loop
    frot n 1e f- f/ to variance
    frot to avg-deviation 
    variance fsqrt { f: standard-deviation }
    variance f0<> if
    n variance fdup f* f* f/ 3e f-
    fswap n variance f*  standard-deviation f* f/ fswap
    endif
    fswap variance standard-deviation avg-deviation n f/ mean ;

 float- -1 floats ,

: partition { first last -- last-smaller first-larger }
    \ partition array addr1 u1 into all elements less than pivot and all
    \ others, addr1 u2 and addr3 u3 are the two partitions.
    \ lessthan-xt ( elemptr1 elemptr2 -- f ) compares the two elements
    first last + 1 rshift faligned f@ { f: pivot }
    first last begin 
    begin
        pivot dup f@ f< over first u> and
    while
        float-
    repeat
    swap begin 
        dup last u< over f@ pivot f< and
    while
        float+
    repeat
    2dup u>=
    while 
    dup f@ over f@ dup f! over f!
    float+ swap float-
    repeat ;

: quantile  recursive
    \ sorts the array [first,last] such that the contained part of
    \ [quant-low,quant-high] is the same as in the fully sorted array.
    { quant-low quant-high }
    begin { first last }
    first quant-high u< quant-low last u< and
    while
    first last partition 
    last quant-low quant-high quantile
    first swap
    repeat ;

: median { addr u -- rmedian }
    addr u 1- 2/ floats + addr u 2/ floats + 
    addr addr u 1- floats + 2over quantile
    f@ f@ f+ f2/ ;


: ff.  
  f$ dup >r 0<
  IF '0 emit ELSE scratch r@ min type r@ precision - zeros ENDIF
  '. emit
  r@ negate zeros
  scratch r> 0 max /string 0 max type ;

create nums \ s" moments.input" r/o open-file throw input-floats
stdin input-floats
nums here over - float /
." n:                  " dup 0 .r cr
compute-loop
nums here over - float / median  9 set-precision
." median:             " ff. cr
." mean:               " ff. cr
." average_deviation:  " ff. cr
." standard_deviation: " ff. cr 11 set-precision
." variance:           " ff. cr  7 set-precision
." skew:               " ff. cr
." kurtosis:           " ff. cr
bye