@n9 %% % WEB SYSTEM : fweb % PROGRAM : gamma.web % AUTHOR : Zdzislaw Meglicki [gustav@@blanc.ovpit.indiana.edu] % CREATION DATE : Mon Oct 5 15:18:03 1998 %% \Title{gamma.web} \author{Zdzis\l aw Meglicki} % web-mode edits the previous two lines when creating a new web. \raggedright \pretolerance=5000 \hfuzz=1in % begin Bottom of Contents Page macro \def\botofcontents{\vskip 0pt plus 1fil minus 1.5in {\bigskip\parskip6pt plus2pt \parindent20pt \vskip0.5in \noindent{\bf Abstract. }\it This document discusses the implementation of the incomplete gamma function. We use the quickly converging continued fraction expansion. } \vskip0.5in {\vfill\it This document illustrates the use of \FWEB\ for the P573 class. } \vfil \rightline{Definitely not for children.} } % end of Bottom of Contents Page macro @* The Continued Fraction Expansion. The purpose of this code is to evaluate \begin{equation} \Gamma(a, x) = \int_x^\infty e^{-t} t^{a-1}\,\textrm{d}t \end{equation} using the continued fraction expansion \begin{eqnarray} \label{e:gamma-fraction} \Gamma(a, x) &=& e^{-x}x^a \frac{1} {x + \frac{1 - a} {1 + \frac{1} {x + \frac{2-a} {1 + \frac{2} {x + \frac{3-a} {1 + \frac{3} {x + \frac{4 - a} {1 + \ldots}}}}}}}} \nonumber\\ &=& e^{-x}x^a \left(\frac{1}{x+}\, \frac{1-a}{1+} \, \frac{1}{x+} \, \frac{2-a}{1+} \frac{2}{x+}\, \frac{3-a}{1+} \, \frac{3}{x+} \, \cdots \, \frac{n-a}{1+} \, \frac{n}{x+}\, \cdots\right) \end{eqnarray} This expansion can be rewritten in terms of the following recursive relationship. Assume that for every finite expansion of order $n$ the fraction can be rewritten as \begin{equation} \label{e:fn} f_n(x) = b_0 + \frac{a_1}{b_1 +}\,\frac{a_2}{b_2 +}\, \cdots \, \frac{a_{n-1}}{b_{n-1}+}\, \frac{a_{n}}{b_{n}} = \frac{A_n}{B_n} \end{equation} It can be proven by induction (and inspected with Maxima) that \begin{eqnarray} \label{e:recurse} A_n &=& b_n A_{n-1} + a_n A_{n-2} \nonumber\\ B_n &=& b_n B_{n-1} + a_n B_{n-2} \end{eqnarray} where \begin{eqnarray} \label{e:initial} A_0 =& b_0\qquad\textrm{and}\qquad A_{-1} &= 1 \nonumber\\ B_0 =& 1\qquad\textrm{and}\qquad B_{-1} &= 0 \end{eqnarray} In our case we get $A_0 = b_0 = 0$, because there is no $b_0$ term in the $\Gamma(a, x)$ expansion. The other terms, i.e., $a_k$ and $b_k$, are listed in table \ref{t:assignments}. \begin{table}[h] \begin{center} \begin{tabular}{@|l@|l@|} % Observe these "at" characters in front % of vertical bars in the WEB input. % A vertical bar on its own marks a % verbatim mode in WEB. \hline $a_1$ & $1$ \\ $b_1$ & $x$ \\ \hline $a_2$ & $1 - a$ \\ $b_2$ & $1$ \\ $a_3$ & $1$ \\ $b_3$ & $x$ \\ \hline $a_4$ & $2 - a$ \\ $b_4$ & $1$ \\ $a_5$ & $2$ \\ $b_5$ & $x$ \\ \hline $\cdots$ & $\cdots$ \\ \hline $a_{2n}$ & $n - a$ \\ $b_{2n}$ & $1$ \\ $a_{2n + 1}$ & $n$ \\ $b_{2n + 1}$ & $x$ \\ \hline \end{tabular} \end{center} \caption{Assignments to $a_k$ and $b_k$ following equations (\ref{e:gamma-fraction}) and (\ref{e:fn}).} \label{t:assignments} \end{table} In this section we outline the mainlines for |function big_gamma(a, x)|: @a real(kind=long) function big_gamma(a, x) real (kind=long), intent(in) :: a, x @ if (x .le. 0.0_long) then if (x .lt. 0.0_long) then stop 'error: big_gamma called with x < 0' /* This is an error condition: we don't define the value of $\Gamma(a, x)$ for $x < 0$ */ else big_gamma = gamma(a) /* This is easy to see, because $\Gamma(a, 0) = \int_0^\infty e^{-t} t^{a-1}\,\textrm{d}t = \Gamma(a)$. */ end if else @ @ big_gamma = r_new * exp(-x + a log(x)) /* This is simply $e^{-x} x^a (\ldots)$ where the stuff within the brackets is stored on |r_new| */ end if end function big_gamma @* Internal Variables. Within the iteration loop we will evaluate first $A_{2n}$, and $B_{2n}$, and then $A_{2n+1}$ and $B_{2n+1}$ using equations (\ref{e:recurse}). The results will be stored on |a0|, and |b0|, and |a1|, and |b1| correspondingly. We will then evaluate $f_{2n+1} = A_{2n+1} / B_{2n+1}$ (cf.~equation~(\ref{e:fn})) and store the result on |r_new|, while the result from the previous evaluation will have been stored on |r_old|. The difference between |r_new| and |r_old| will be called |delta|: $\delta = f_{2n+1} - f_{2n-1}$. The iteration will stop when $\vert\delta\vert \le \vert\epsilon \times f_{2n+1}\vert$, where $\epsilon = 10^{-9}$. The integer |n| corresponds to index $n$. Table \ref{t:variables} sums up the meaning of variables used by function |big_gamma|. \begin{table}[h] \begin{center} \begin{tabular}{@|l@|l@|} % Watch those ``at''-bars again! \hline |n| & $n$ \\ |a0| & $A_{2n}$ \\ |b0| & $B_{2n}$ \\ |a1| & $A_{2n+1}$ \\ |b1| & $B_{2n+1}$ \\ |r_new| & $f_{2n+1} = A_{2n+1}/B_{2n+1}$ \\ |r_old| & $f_{2n-1} = A_{2n-1}/B_{2n-1}$ \\ |delta| & $\delta = f_{2n+1} - f_{2n-1}$ \\ |epsilon| & $\epsilon = 10^{-9}$ \\ \hline \end{tabular} \end{center} \caption{The meaning of variables and parameters used internally by function |big_gamma|} \label{t:variables} \end{table} @= real (kind=long), parameter :: epsilon = 1.0E-9 real (kind=long) :: a0, b0, a1, b1, r_old, r_new, delta integer :: n @* Initialize Iterations. In principle we can begin iteration from any $n$ assuming that we would have evaluated $A_{n-1}$, $B_{n-1}$, $A_{n-2}$, and $B_{n-2}$ by some other means. Table \ref{t:assignments} shows that proper sequencing in our case begins only from $a_2$ and $b_2$. Consequently, in order to commence the iteration, we have to evaluate $A_1$, and $B_1$ manually. Remember that equations (\ref{e:initial}) already give us $A_{-1} = 1$, $B_{-1} = 0$, $A_0 = b_0 = 0$, and $B_0 = 1$. Using equations (\ref{e:recurse}) we can easily find that \begin{eqnarray*} A_1 &=& b_1 A_0 + a_1 A_{-1} = x\cdot 0 + 1\cdot 1 = 1\\ B_1 &=& B_1 = b_1 B_0 + a_1 B_{-1} = x\cdot 1 + 1\cdot 0 = x \end{eqnarray*} Table \ref{t:initial} sums up the initial assignments. \begin{table}[h] \begin{center} \begin{tabular}{@|l@|l@|l@|} \hline $A_0$ & $0$ & |a0| \\ $B_0$ & $1$ & |b0| \\ $A_1$ & $1$ & |a1| \\ $B_1$ & $x$ & |b1| \\ $f_1$ & $A_1/B_1$ & |r_old| \\ $\delta$ & $A_1/B_1$ & |delta| \\ $n$ & $1$ & |n| \\ \hline \end{tabular} \end{center} \caption{Initial assignments} \label{t:initial} \end{table} @= a0 = 0.0_long b0 = 1.0_long a1 = 1.0_long b1 = x r_old = a1/b1 delta = r_old n = 1 @* Iterate. The iteration now corresponds closely to equations (\ref{e:recurse}). Remember that as we enter the loop, |a1| and |b1| store the previous values, i.e., $A_{2n-1}$ and $B_{2n-1}$, and |a0| and |b0| store even older values of $A_{2n-2}$ and $B_{2n-2}$. Using assignments from table \ref{t:assignments} in place of $a_{2n}$, $b_{2n}$, $a_{2n+1}$ and $b_{2n+1}$, we get the following: \begin{eqnarray*} A_{2n} &=& b_{2n} A_{2n-1} + a_{2n} A_{2n-2} = 1 \times A_{2n-1} + (n - a) \times A_{2n-2} \\ B_{2n} &=& b_{2n} B_{2n-1} + a_{2n} B_{2n-2} = 1 \times B_{2n-1} + (n - a) \times B_{2n-2} \end{eqnarray*} But now the roles change, and |a0| and |b0| contain the latest updated values, i.e., $A_{2n}$ and $B_{2n}$, whereas |a1| and |b1| still contain $A_{2n-1}$ and $B_{2n-1}$. \begin{eqnarray*} A_{2n+1} &=& b_{2n+1} A_{2n} + a_{2n+1} A_{2n-1} = x \times A_{2n} + n \times A_{2n-1} \\ B_{2n+1} &=& b_{2n+1} B_{2n} + a_{2n+1} B_{2n-1} = x \times B_{2n} + n \times B_{2n-1} \end{eqnarray*} Observe that {\em two\/} steps are performed within each iteration. After we have found $A_{2n+1}$ and $B_{2n+1}$ we evaluate $f_{2n+1} = A_{2n+1} / B_{2n+1}$, calculate the difference $f_{2n+1} - f_{2n-1} = A_{2n+1} / B_{2n+1} - A_{2n-1} / B_{2n-1}$, drop the old value $f_{2n-1}$ and store the new value $f_{2n+1}$ on |r_old|. The iteration stops when $\vert f_{2n+1} - f_{2n-1}\vert \le \vert\epsilon\times f_{2n+1}\vert$. @= do while (abs(delta) .gt. abs(epsilon * r_old)) a0 = a1 + (n - a) * a0 /* $A_{2n} = b_{2n} A_{2n-1} + a_{2n} A_{2n-2}$ */ b0 = b1 + (n - a) * b0 /* $B_{2n} = b_{2n} B_{2n-1} + a_{2n} B_{2n-2}$ */ a1 = x * a0 + n * a1 /* $A_{2n+1} = b_{2n+1} A_{2n} + a_{2n+1} A_{2n-1}$ */ b1 = x * b0 + n * b1 /* $B_{2n+1} = b_{2n+1} B_{2n} + a_{2n+1} B_{2n-1}$ */ r_new = a1/b1 /* $f_{2n+1} = A_{2n+1}/B_{2n+1}$ */ delta = r_new - r_old /* $\delta = f_{2n+1} - f_{2n-1}$ */ r_old = r_new n = n + 1 end do @* INDEX.