Computing mean and population standard deviation in Forth.

stddev.fs

\ Computing mean and population standard deviation in Forth.
\ Tested on Gforth 0.7.0

16 constant FRAC \ number of bits behind fixed-point

\ convert from/to fixed point representation
: s>fix FRAC lshift ;
: fix>s FRAC rshift ;

\ fixed-point arithmetic
1 s>fix constant fix1 \ representation of 1 in fixed-point

: fix/ ( x y -- z ) swap s>fix swap / ;
: fix* ( x y -- z ) * fix>s ;

: fixinv ( x -- x ) fix1 swap fix/ ;

: fixsquare ( x -- x^2 )
  dup fix*
;

\ Convert fix to float (for easy output).
\ Of course the fact that this uses floating-point probably
\ makes the point of using fixed-point in the first place silly.
: fix>f ( x -- d )
  s>f [ fix1 s>f ] fliteral f/
;

\ Print fixpoint nicely.
: fix. ( x -- )
  fix>f f.
;

\ fixed-point square root using Newton's method
: newton-step ( c x -- c x ) 
  >r dup r@ fix/ r> + 1 rshift 
;

: fixsqrt ( c -- x ) 
  fix1 \ take 1 as initial guess
  begin 
    \ Note: due to rounding, this may end up cycling between two numbers.
    \ Therefore, we compare with the second-to-last step for termination.
    dup >r newton-step newton-step dup r> = 
  until 
  swap drop
;

: mean ( sum inv_len -- mean ) fix* ; \ trivial but useful

\ Gives *population* stddev.
: std ( sum2 inv_len mean -- std )
  \ std = √(sum2 * inv_len - mean²)
  fixsquare
  >r fix* r> 
  - fixsqrt
;

\ Computes both mean and population stddev
\ The abundance of stack manipulation words in the following
\ suggests the order of arguments and return values is non-optimal.
: mean_std ( sum2 sum inv_len -- mean std )
  dup >r mean r> swap ( sum2 inv_len mean )
  dup >r std r> swap
;

\ Example: 10 numbers
\ 28, 76, 15, 53, 24, 88, 75, 80, 72, 17
528 constant sum
35412 constant sum2
10 constant len

sum2 s>fix sum s>fix fix1 len / 
mean_std 
." std: " fix.  cr \ expected value: 27.447404248853843
." mean:" fix.  cr \ expected value: 52.8