Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_specfun_1f1_real_scaled (s22bb)

## Purpose

nag_specfun_1f1_real_scaled (s22bb) returns a value for the confluent hypergeometric function 1F1 (a ; b ; x) ${}_{1}F_{1}\left(a;b;x\right)$ with real parameters a$a$, b$b$ and real argument x$x$ in the scaled form 1F1 (a ; b ; x) = mf × 2ms ${}_{1}F_{1}\left(a;b;x\right)={m}_{f}×{2}^{{m}_{s}}$. This function is sometimes also known as Kummer's function M(a,b,x)$M\left(a,b,x\right)$.

## Syntax

[frm, scm, ifail] = s22bb(ani, adr, bni, bdr, x)
[frm, scm, ifail] = nag_specfun_1f1_real_scaled(ani, adr, bni, bdr, x)

## Description

nag_specfun_1f1_real_scaled (s22bb) returns a value for the confluent hypergeometric function 1F1 (a ; b ; x) ${}_{1}F_{1}\left(a;b;x\right)$ with real parameters a$a$, b$b$ and x$x$ in the scaled form 1F1 (a ; b ; x) = mf × 2ms ${}_{1}F_{1}\left(a;b;x\right)={m}_{f}×{2}^{{m}_{s}}$, where mf${m}_{f}$ is the real scaled component and ms${m}_{s}$ is the integer power of two scaling. This function is unbounded or not uniquely defined for b$b$ equal to zero or a negative integer.
The confluent hypergeometric function is defined by the confluent series,
 ∞ 1F1(a ; b ; x) = M(a,b,x) = ∑ ( (a)s xs )/( (b)s s! ) = 1 + (a)/b x + ( a(a + 1) )/( b(b + 1) 2! )x2 + ⋯ s = 0
$F1 1 (a;b;x) = M(a,b,x) = ∑ s=0 ∞ (a)s xs (b)s s! = 1 + a b x + a(a+1) b(b+1) 2! x2 + ⋯$
where (a)s = 1 (a) (a + 1) (a + 2) (a + s1) ${\left(a\right)}_{s}=1\left(a\right)\left(a+1\right)\left(a+2\right)\dots \left(a+s-1\right)$ is the rising factorial of a$a$. M(a,b,x) $M\left(a,b,x\right)$ is a solution to the second order ODE (Kummer's Equation):
 x (d2M)/(dx2) + (b − x) (dM)/(dx) − a M = 0 . $x d2M dx2 + (b-x) dM dx - a M = 0 .$ (1)
Given the parameters and argument (a,b,x) $\left(a,b,x\right)$, this function determines a set of safe values {(αi,βi,ζi)i2} $\left\{\left({\alpha }_{i},{\beta }_{i},{\zeta }_{i}\right)\mid i\le 2\right\}$ and selects an appropriate algorithm to accurately evaluate the functions Mi (αi,βi,ζi) ${M}_{i}\left({\alpha }_{i},{\beta }_{i},{\zeta }_{i}\right)$. The result is then used to construct the solution to the original problem M(a,b,x) $M\left(a,b,x\right)$ using, where necessary, recurrence relations and/or continuation.
For improved precision in the final result, this function accepts a$a$ and b$b$ split into an integral and a decimal fractional component. Specifically a = ai + ar$a={a}_{i}+{a}_{r}$, where |ar|0.5$|{a}_{r}|\le 0.5$ and ai = aar${a}_{i}=a-{a}_{r}$ is integral. b$b$ is similarly deconstructed.
Additionally, an artificial bound, arbnd$\mathit{arbnd}$ is placed on the magnitudes of ai${a}_{i}$, bi${b}_{i}$ and x$x$ to minimize the occurrence of overflow in internal calculations. arbnd = 0.0001 × Imax $\mathit{arbnd}=0.0001×{I}_{\mathrm{max}}$, where Imax = x02bb${I}_{\mathrm{max}}=\mathbf{x02bb}$. It should, however, not be assumed that this function will produce an accurate result for all values of ai${a}_{i}$, bi${b}_{i}$ and x$x$ satisfying this criterion.
Please consult the NIST Digital Library of Mathematical Functions or the companion (2010) for a detailed discussion of the confluent hypergeometric function including special cases, transformations, relations and asymptotic approximations.

## References

NIST Handbook of Mathematical Functions (2010) (eds F W J Olver, D W Lozier, R F Boisvert, C W Clark) Cambridge University Press
Pearson J (2009) Computation of hypergeometric functions MSc Dissertation, Mathematical Institute, University of Oxford

## Parameters

### Compulsory Input Parameters

1:     ani – double scalar
ai${a}_{i}$, the nearest integer to a$a$, satisfying ai = aar${a}_{i}=a-{a}_{r}$.
Constraints:
• ani = ani${\mathbf{ani}}=⌊{\mathbf{ani}}⌋$;
• |ani|arbnd$|{\mathbf{ani}}|\le \mathit{arbnd}$.
ar${a}_{r}$, the signed decimal remainder satisfying ar = aai ${a}_{r}=a-{a}_{i}$ and |ar| 0.5$|{a}_{r}|\le 0.5$.
Constraint: |adr|0.5$|{\mathbf{adr}}|\le 0.5$.
Note: if |adr| < 100.0ε$|{\mathbf{adr}}|<100.0\epsilon$, ar = 0.0${a}_{r}=0.0$ will be used, where ε$\epsilon$ is the machine precision as returned by nag_machine_precision (x02aj).
3:     bni – double scalar
bi${b}_{i}$, the nearest integer to b$b$, satisfying bi = bbr${b}_{i}=b-{b}_{r}$.
Constraints:
• bni = bni${\mathbf{bni}}=⌊{\mathbf{bni}}⌋$;
• |bni|arbnd$|{\mathbf{bni}}|\le \mathit{arbnd}$;
• if bdr = 0.0${\mathbf{bdr}}=0.0$, bni > 0${\mathbf{bni}}>0$.
4:     bdr – double scalar
br${b}_{r}$, the signed decimal remainder satisfying br = bbi${b}_{r}=b-{b}_{i}$ and |br| 0.5$|{b}_{r}|\le 0.5$.
Constraint: |bdr|0.5$|{\mathbf{bdr}}|\le 0.5$.
Note: if |bdradr| < 100.0ε$|{\mathbf{bdr}}-{\mathbf{adr}}|<100.0\epsilon$, ar = br${a}_{r}={b}_{r}$ will be used, where ε$\epsilon$ is the machine precision as returned by nag_machine_precision (x02aj).
5:     x – double scalar
The argument x$x$ of the function.
Constraint: |x|arbnd$|{\mathbf{x}}|\le \mathit{arbnd}$.

None.

None.

### Output Parameters

1:     frm – double scalar
mf${m}_{f}$, the scaled real component of the solution satisfying mf = M(a,b,x) × 2ms${m}_{f}=M\left(a,b,x\right)×{2}^{-{m}_{s}}$.
Note: if overflow occurs upon completion, as indicated by ${\mathbf{ifail}}={\mathbf{2}}$, the value of mf${m}_{f}$ returned may still be correct. If overflow occurs in a subcalculation, as indicated by ${\mathbf{ifail}}={\mathbf{5}}$, this should not be assumed.
2:     scm – int64int32nag_int scalar
ms${m}_{s}$, the scaling power of two, satisfying ms = log2((M(a,b,x))/(mf))${m}_{s}={\mathrm{log}}_{2}\left(\frac{M\left(a,b,x\right)}{{m}_{f}}\right)$.
Note: if overflow occurs upon completion, as indicated by ${\mathbf{ifail}}={\mathbf{2}}$, then msImax${m}_{s}\ge {I}_{\mathrm{max}}$, where Imax${I}_{\mathrm{max}}$ is the largest representable integer (see nag_machine_integer_max (x02bb)). If overflow occurs during a subcalculation, as indicated by ${\mathbf{ifail}}={\mathbf{5}}$, ms${m}_{s}$ may or may not be greater than Imax${I}_{\mathrm{max}}$. In either case, ${\mathbf{scm}}=\mathbf{x02bb}$ will have been returned.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W ifail = 1${\mathbf{ifail}}=1$
Underflow occurred during the evaluation of M(a,b,x)$M\left(a,b,x\right)$.
The returned value may be inaccurate.
W ifail = 2${\mathbf{ifail}}=2$
On completion, overflow occurred in the evaluation of M(a,b,x)$M\left(a,b,x\right)$.
W ifail = 3${\mathbf{ifail}}=3$
All approximations have completed, and the final residual estimate indicates some precision may have been lost.
ifail = 4${\mathbf{ifail}}=4$
All approximations have completed, and the final residual estimate indicates no accuracy can be guaranteed.
ifail = 5${\mathbf{ifail}}=5$
Overflow occurred in a subcalculation of M(a,b,x)$M\left(a,b,x\right)$.
The answer may be completely incorrect.
ifail = 11${\mathbf{ifail}}=11$
Constraint: .
ifail = 13${\mathbf{ifail}}=13$
Constraint: ani = ani${\mathbf{ani}}=⌊{\mathbf{ani}}⌋$.
ifail = 21${\mathbf{ifail}}=21$
Constraint: |adr|0.5$|{\mathbf{adr}}|\le 0.5$.
ifail = 31${\mathbf{ifail}}=31$
Constraint: .
ifail = 32${\mathbf{ifail}}=32$
On entry.
M(a,b,x)$M\left(a,b,x\right)$ is undefined when b$b$ is zero or a negative integer.
ifail = 33${\mathbf{ifail}}=33$
Constraint: bni = bni${\mathbf{bni}}=⌊{\mathbf{bni}}⌋$.
ifail = 41${\mathbf{ifail}}=41$
Constraint: |bdr|0.5$|{\mathbf{bdr}}|\le 0.5$.
ifail = 51${\mathbf{ifail}}=51$
Constraint: .

## Accuracy

In general, if ${\mathbf{ifail}}={\mathbf{0}}$, the value of M$M$ may be assumed accurate, with the possible loss of one or two decimal places. Assuming the result does not under or overflow, an error estimate res$\mathit{res}$ is made internally using equation (1). If the magnitude of res$\mathit{res}$ is sufficiently large a nonzero ifail will be returned. Specifically,
 ${\mathbf{ifail}}={\mathbf{0}}$ res ≤ 1000ε$\mathit{res}\le 1000\epsilon$ ${\mathbf{ifail}}={\mathbf{3}}$ 1000 ε < res ≤ 0.1$1000\epsilon <\mathit{res}\le 0.1$ ${\mathbf{ifail}}={\mathbf{4}}$ res > 0.1$\mathit{res}>0.1$
A further estimate of the residual can be constructed using equation (1), and the differential identity,
 ( d M(a,b,x) )/(dx) = a/b M (a + 1,b + 1,x) , ( d2 M(a,b,x) )/(dx2) = (a(a + 1))/(b(b + 1)) M (a + 2,b + 2,x) .
$d M(a,b,x) dx = ab M (a+1,b+1,x) , d2 M(a,b,x) dx2 = a(a+1) b(b+1) M (a+2,b+2,x) .$
This estimate is however dependent upon the error involved in approximating M (a + 1,b + 1,x) $M\left(a+1,b+1,x\right)$ and M (a + 2,b + 2,x) $M\left(a+2,b+2,x\right)$.

The values returned in frm (mf${m}_{f}$) and scm (ms${m}_{s}$) may be used to explicitly evaluate M(a,b,x)$M\left(a,b,x\right)$, and may also be used to evaluate products and ratios of multiple values of M$M$ as follows,
 M(a,b,x) = mf × 2ms M (a1,b1,x1) × M (a2,b2,x2) = (mf1 × mf2) × 2(ms1 + ms2) ( M (a1,b1,x1) )/( M (a2,b2,x2) ) = (mf1)/(mf2) × 2(ms1 − ms2) ln|M(a,b,x)| = ln|mf| + ms × ln(2)
.
$M(a,b,x) = mf × 2ms M (a1,b1,x1) × M (a2,b2,x2) = ( mf1 × mf2 ) × 2 ( ms1 + ms2 ) M (a1,b1,x1) M (a2,b2,x2) = mf1 mf2 × 2 ( ms1 - ms2 ) ln| M (a,b,x) | = ln|mf| + ms × ln(2) .$

## Example

function nag_specfun_1f1_real_scaled_example
ai = -10;
bi = 30;
delta = 1e-4;
ar = delta;
br = -delta;
x = 25;
frmv = zeros(2, 1);
scmv = zeros(2, 1, 'int64');
fprintf('\n         a          b          x          frm    scm     M(a,b,x)\n');
[frmv(1), scmv(1), ifail] = nag_specfun_1f1_real_scaled(ai,  ar, bi,  br, x);
print_result(ai+ar, bi+br, x, frmv(1), scmv(1), ifail);
[frmv(2), scmv(2), ifail] = nag_specfun_1f1_real_scaled(ai, -ar, bi, -br, x);
print_result(ai-ar, bi-br, x, frmv(2), scmv(2), ifail);

% Calculate the product M1*M2
frm = prod(frmv);
scm = scmv(1)+scmv(2);

if scm < nag_machine_model_maxexp
scale = frm*double(nag_machine_model_base)^double(scm);
fprintf('\nSolution product %12.4e %6d %12.4e\n', frm, scm, scale);
else
fprintf('\nSolution product %12.4e %6d Not representable\n', frm, scm);
end

% Calculate the ratio M1/M2
if frmv(2) ~= 0
frm = frmv(1)/frmv(2);
scm = scmv(1) - scmv(2);
if scm < nag_machine_model_maxexp
scale = frm*double(nag_machine_model_base)^double(scm);
fprintf('\nSolution ratio   %12.4e %6d %12.4e\n', frm, scm, scale);
else
fprintf('\nSolution ratio   %12.4e %6d Not representable\n', frm, scm);
end
end

function print_result(a, b, x, frm, scm, ifail)
if ifail ==2 || ifail > 3
% Either the result has overflowed, no accuracy may be assumed,
% or an input error has been detected.
fprintf('%10.4f %10.4f %10.4f                      FAILED\n', a, b, x);
elseif scm < nag_machine_model_maxexp
scale = frm*double(nag_machine_model_base)^double(scm);
fprintf('%10.4f %10.4f %10.4f %12.4e %6d %12.4e\n', a, b, x, frm, scm, scale);
else
fprintf('%10.4f %10.4f %10.4f %12.4e %6d Not representable\n', a, b, x, frm, scm);
end

a          b          x          frm    scm     M(a,b,x)
-9.9999    29.9999    25.0000  -7.7329e-01    -15  -2.3599e-05
-10.0001    30.0001    25.0000  -7.7318e-01    -15  -2.3596e-05

Solution product   5.9789e-01    -30   5.5683e-10

Solution ratio     1.0001e+00      0   1.0001e+00

function s22bb_example
ai = -10;
bi = 30;
delta = 1e-4;
ar = delta;
br = -delta;
x = 25;
frmv = zeros(2, 1);
scmv = zeros(2, 1, 'int64');
fprintf('\n         a          b          x          frm    scm     M(a,b,x)\n');
[frmv(1), scmv(1), ifail] = s22bb(ai,  ar, bi,  br, x);
print_result(ai+ar, bi+br, x, frmv(1), scmv(1), ifail);
[frmv(2), scmv(2), ifail] = s22bb(ai, -ar, bi, -br, x);
print_result(ai-ar, bi-br, x, frmv(2), scmv(2), ifail);

% Calculate the product M1*M2
frm = prod(frmv);
scm = scmv(1)+scmv(2);

if scm < x02bl
scale = frm*double(x02bh)^double(scm);
fprintf('\nSolution product %12.4e %6d %12.4e\n', frm, scm, scale);
else
fprintf('\nSolution product %12.4e %6d Not representable\n', frm, scm);
end

% Calculate the ratio M1/M2
if frmv(2) ~= 0
frm = frmv(1)/frmv(2);
scm = scmv(1) - scmv(2);
if scm < x02bl
scale = frm*double(x02bh)^double(scm);
fprintf('\nSolution ratio   %12.4e %6d %12.4e\n', frm, scm, scale);
else
fprintf('\nSolution ratio   %12.4e %6d Not representable\n', frm, scm);
end
end

function print_result(a, b, x, frm, scm, ifail)
if ifail ==2 || ifail > 3
% Either the result has overflowed, no accuracy may be assumed,
% or an input error has been detected.
fprintf('%10.4f %10.4f %10.4f                      FAILED\n', a, b, x);
elseif scm < x02bl
scale = frm*double(x02bh)^double(scm);
fprintf('%10.4f %10.4f %10.4f %12.4e %6d %12.4e\n', a, b, x, frm, scm, scale);
else
fprintf('%10.4f %10.4f %10.4f %12.4e %6d Not representable\n', a, b, x, frm, scm);
end

a          b          x          frm    scm     M(a,b,x)
-9.9999    29.9999    25.0000  -7.7329e-01    -15  -2.3599e-05
-10.0001    30.0001    25.0000  -7.7318e-01    -15  -2.3596e-05

Solution product   5.9789e-01    -30   5.5683e-10

Solution ratio     1.0001e+00      0   1.0001e+00