extra/math/analysis/analysis.factor

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