extra/math/extras/extras.factor

   1 ! Copyright (C) 2012 John Benediktsson
   2 ! See http://factorcode.org/license.txt for BSD license
   3
   4 USING: combinators.short-circuit grouping kernel math
   5 math.combinatorics math.functions math.order math.primes
   6 math.ranges math.statistics math.vectors memoize sequences ;
   7
   8 IN: math.extras
   9
  10 <PRIVATE
  11
  12 DEFER: sterling
  13
  14 : (sterling) ( n k -- x )
  15     [ [ 1 - ] bi@ sterling ]
  16     [ [ 1 - ] dip sterling ]
  17     [ nip * + ] 2tri ;
  18
  19 PRIVATE>
  20
  21 MEMO: sterling ( n k -- x )
  22     2dup { [ = ] [ nip 1 = ] } 2||
  23     [ 2drop 1 ] [ (sterling) ] if ;
  24
  25 <PRIVATE
  26
  27 DEFER: bernoulli
  28
  29 : (bernoulli) ( p -- n )
  30     [ iota ] [ 1 + ] bi [
  31         0 [ [ nCk ] [ bernoulli * ] bi + ] with reduce
  32     ] keep recip neg * ;
  33
  34 PRIVATE>
  35
  36 MEMO: bernoulli ( p -- n )
  37     [ 1 ] [ (bernoulli) ] if-zero ;
  38
  39 : chi2 ( actual expected -- n )
  40     0 [ dup 0 > [ [ - sq ] keep / + ] [ 2drop ] if ] 2reduce ;
  41
  42 <PRIVATE
  43
  44 : df-check ( df -- )
  45     even? [ "odd degrees of freedom" throw ] unless ;
  46
  47 : (chi2P) ( chi/2 df/2 -- p )
  48     [1,b) dupd n/v cum-product swap neg e^ [ v*n sum ] keep + ;
  49
  50 PRIVATE>
  51
  52 : chi2P ( chi df -- p )
  53     dup df-check [ 2.0 / ] [ 2 /i ] bi* (chi2P) 1.0 min ;
  54
  55 <PRIVATE
  56
  57 : check-jacobi ( m -- m )
  58     dup { [ integer? ] [ 0 > ] [ odd? ] } 1&&
  59     [ "modulus must be odd positive integer" throw ] unless ;
  60
  61 : mod' ( x y -- n )
  62     [ mod ] keep over zero? [ drop ] [
  63         2dup [ sgn ] bi@ = [ drop ] [ + ] if
  64     ] if ;
  65
  66 PRIVATE>
  67
  68 : jacobi ( a m -- n )
  69     check-jacobi [ mod' ] keep 1
  70     [ pick zero? ] [
  71         [ pick even? ] [
  72             [ 2 / ] 2dip
  73             over 8 mod' { 3 5 } member? [ neg ] when
  74         ] while swapd
  75         2over [ 4 mod' 3 = ] both? [ neg ] when
  76         [ [ mod' ] keep ] dip
  77     ] until [ nip 1 = ] dip 0 ? ;
  78
  79 <PRIVATE
  80
  81 : check-legendere ( m -- m )
  82     dup prime? [ "modulus must be prime positive integer" throw ] unless ;
  83
  84 PRIVATE>
  85
  86 : legendere ( a m -- n )
  87     check-legendere jacobi ;
  88
  89 : moving-average ( seq n -- newseq )
  90     <clumps> [ mean ] map ;
  91
  92 : exponential-moving-average ( seq a -- newseq )
  93     [ 1 ] 2dip [ [ dupd swap - ] dip * + dup ] curry map nip ;
  94
  95 : moving-median ( u n -- v )
  96     <clumps> [ median ] map ;
  97
  98 : nonzero ( seq -- seq' )
  99     [ zero? not ] filter ;