extra/math/factorials/factorials.factor

   1 ! Copyright (C) 2013 John Benediktsson
   2 ! See http://factorcode.org/license.txt for BSD license
   3
   4 USING: combinators combinators.short-circuit inverse kernel
   5 math math.functions math.primes ranges sequences ;
   6
   7 IN: math.factorials
   8
   9 MEMO: factorial ( n -- n! )
  10     dup 1 > [ [1..b] product ] [ drop 1 ] if ;
  11
  12 ALIAS: n! factorial
  13
  14 : factorials ( n -- seq )
  15     1 swap [0..b] [ dup 1 > [ * ] [ drop ] if dup ] map nip ;
  16
  17 MEMO: double-factorial ( n -- n!! )
  18     dup [ even? ] [ 0 < ] bi [
  19         [ drop 1/0. ] [
  20             2 + -1 swap -2 <range> product recip
  21         ] if
  22     ] [
  23         2 3 ? swap 2 <range> product
  24     ] if ;
  25
  26 ALIAS: n!! double-factorial
  27
  28 : factorial/ ( n k -- n!/k! )
  29     {
  30         { [ dup 1 <= ] [ drop factorial ] }
  31         { [ over 1 <= ] [ nip factorial recip ] }
  32         [
  33             2dup < [ t ] [ swap f ] if
  34             [ (a..b] product ] dip [ recip ] when
  35         ]
  36     } cond ;
  37
  38 : rising-factorial ( x n -- x(n) )
  39     {
  40         { 1 [ ] }
  41         { 0 [ drop 0 ] }
  42         [
  43             dup 0 < [ neg [ + ] keep t ] [ f ] if
  44             [ dupd + [a..b) product ] dip
  45             [ recip ] when
  46         ]
  47     } case ;
  48
  49 ALIAS: pochhammer rising-factorial
  50
  51 : falling-factorial ( x n -- (x)n )
  52     {
  53         { 1 [ ] }
  54         { 0 [ drop 0 ] }
  55         [
  56             dup 0 < [ neg [ + ] keep t ] [ f ] if
  57             [ dupd - swap (a..b] product ] dip
  58             [ recip ] when
  59         ]
  60     } case ;
  61
  62 : factorial-power ( x n h -- (x)n(h) )
  63     {
  64         { 1 [ falling-factorial ] }
  65         { 0 [ ^ ] }
  66         [
  67             over 0 < [
  68                 [ [ nip + ] [ swap neg * + ] 3bi ] keep
  69                 <range> product recip
  70             ] [
  71                 neg [ [ dupd 1 - ] [ * ] bi* + ] keep
  72                 <range> product
  73             ] if
  74         ]
  75     } case ;
  76
  77 : primorial ( n -- p# )
  78     dup 0 > [ nprimes product ] [ drop 1 ] if ;
  79
  80 : multifactorial ( n k -- n!(k) )
  81     2dup >= [
  82         dupd [ - ] keep multifactorial *
  83     ] [ 2drop 1 ] if ; inline recursive
  84
  85 : quadruple-factorial ( n -- m )
  86     [ 2 * ] keep factorial/ ;
  87
  88 : super-factorial ( n -- m )
  89     dup 1 > [
  90         [1..b] [ factorial ] [ * ] map-reduce
  91     ] [ drop 1 ] if ;
  92
  93 : hyper-factorial ( n -- m )
  94     dup 1 > [
  95         [1..b] [ dup ^ ] [ * ] map-reduce
  96     ] [ drop 1 ] if ;
  97
  98 : alternating-factorial ( n -- m )
  99     dup 1 > [
 100         [ [1..b] ] keep even? '[
 101             [ factorial ] [ odd? _ = ] bi [ neg ] when
 102         ] map-sum
 103     ] [ drop 1 ] if ;
 104
 105 : exponential-factorial ( n -- m )
 106     dup 1 > [ [1..b] 1 [ swap ^ ] reduce ] [ drop 1 ] if ;
 107
 108 <PRIVATE
 109
 110 : -prime? ( n quot: ( n -- m ) -- ? )
 111     [ 1 1 [ pick over - 1 <= ] ] dip
 112     '[ drop [ 1 + ] _ bi ] until nip - abs 1 = ; inline
 113
 114 PRIVATE>
 115
 116 : factorial-prime? ( n -- ? )
 117     { [ prime? ] [ [ factorial ] -prime? ] } 1&& ;
 118
 119 : primorial-prime? ( n -- ? )
 120     { [ prime? ] [ 2 > ] [ [ primorial ] -prime? ] } 1&& ;
 121
 122 : reverse-factorial ( m -- n )
 123     1 1 [ 2over > ] [ 1 + [ * ] keep ] while [ = ] dip and ;
 124
 125 \ factorial [ reverse-factorial ] define-inverse