extra/math/analysis/analysis.factor

   1 ! Copyright (C) 2008 Doug Coleman, Slava Pestov, Aaron Schaefer.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: combinators.short-circuit kernel math math.constants
   4 math.functions math.vectors sequences ;
   5 IN: math.analysis
   6
   7 <PRIVATE
   8
   9 ! http://www.rskey.org/gamma.htm  "Lanczos Approximation"
  10 ! n=6: error ~ 3 x 10^-11
  11
  12 CONSTANT: gamma-g6 5.15
  13
  14 CONSTANT: gamma-p6
  15     {
  16         2.50662827563479526904 225.525584619175212544 -268.295973841304927459
  17         80.9030806934622512966 -5.00757863970517583837 0.0114684895434781459556
  18     }
  19
  20 : gamma-z ( x n -- seq )
  21     [ + recip ] with map 1.0 0 pick set-nth ;
  22
  23 : (gamma-lanczos6) ( x -- log[gamma[x+1]] )
  24     #! log(gamma(x+1)
  25     [ 0.5 + dup gamma-g6 + [ log * ] keep - ]
  26     [ 6 gamma-z gamma-p6 v. log ] bi + ;
  27
  28 : gamma-lanczos6 ( x -- gamma[x] )
  29     #! gamma(x) = gamma(x+1) / x
  30     [ (gamma-lanczos6) exp ] keep / ;
  31
  32 : gammaln-lanczos6 ( x -- gammaln[x] )
  33     #! log(gamma(x)) = log(gamma(x+1)) - log(x)
  34     [ (gamma-lanczos6) ] keep log - ;
  35
  36 : gamma-neg ( gamma[abs[x]] x -- gamma[x] )
  37     dup pi * sin * * pi neg swap / ; inline
  38
  39 PRIVATE>
  40
  41 : gamma ( x -- y )
  42     #! gamma(x) = integral 0..inf [ t^(x-1) exp(-t) ] dt
  43     #! gamma(n+1) = n! for n > 0
  44     dup { [ 0.0 <= ] [ 1.0 mod zero? ] } 1&& [
  45         drop 1/0.
  46     ] [
  47         [ abs gamma-lanczos6 ] keep dup 0 > [ drop ] [ gamma-neg ] if
  48     ] if ;
  49
  50 : gammaln ( x -- gamma[x] )
  51     #! gammaln(x) is an alternative when gamma(x)'s range
  52     #! varies too widely
  53     dup 0 < [
  54         drop 1/0.
  55     ] [
  56         [ abs gammaln-lanczos6 ] keep dup 0 > [ drop ] [ gamma-neg ] if
  57     ] if ;
  58
  59 : nth-root ( n x -- y )
  60     swap recip ^ ;
  61
  62 ! Forth Scientific Library Algorithm #1
  63 !
  64 ! Evaluates the Real Exponential Integral,
  65 !     E1(x) = - Ei(-x) =   int_x^\infty exp^{-u}/u du      for x > 0
  66 ! using a rational approximation
  67 !
  68 ! Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
  69 ! 1980; Association for Computing Machinery Inc., New York,
  70 ! ISBN 0-89791-017-6
  71 !
  72 ! (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
  73 ! author to use this software for any application provided the
  74 ! copyright notice is preserved.
  75
  76 : exp-int ( x -- y )
  77     #! For real values of x only. Accurate to 7 decimals.
  78     dup 1.0 < [
  79         dup 0.00107857 * 0.00976004 -
  80         over *
  81         0.05519968 +
  82         over *
  83         0.24991055 -
  84         over *
  85         0.99999193 +
  86         over *
  87         0.57721566 -
  88         swap log -
  89     ] [
  90         dup 8.5733287401 +
  91         over *
  92         18.059016973 +
  93         over *
  94         8.6347608925 +
  95         over *
  96         0.2677737343 +
  97
  98         over
  99         dup 9.5733223454 +
 100         over *
 101         25.6329561486 +
 102         over *
 103         21.0996530827 +
 104         over *
 105         3.9584969228 +
 106
 107         nip
 108         /
 109         over /
 110         swap -1.0 * exp
 111         *
 112     ] if ;
 113
 114 ! James Stirling's approximation for N!:
 115 ! http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Numerical/Stirling/
 116
 117 : stirling-fact ( n -- fact )
 118     [ pi 2 * * sqrt ]
 119     [ [ e / ] keep ^ ]
 120     [ 12 * recip 1 + ] tri * * ;
 121