Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Computing mean and population standard deviation in Forth.
\ 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.