[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

1136.0. "" by CLADA::KEATING (rocking the bay) Sun Oct 08 1989 14:47

T.RTitleUserPersonal
Name
DateLines
1136.1Clarification of the problemCLADA::KEATINGrocking the baySun Oct 08 1989 17:18219
1136.2AITG::DERAMODaniel V. {AITG,ZFC}:: D'EramoSun Oct 08 1989 18:5012
        You can bound the integral from x to oo of e^(-t^2)dt
        by noting that it is less than the integral from x to oo
        of e^(-xt)dt which is (e^(-x^2))/x.  I think by the time
        you get to x=6 the area under the tail of the curve is
        less than 10^-12.  So you can determine a point beyond
        which you might as well take the answer as being 1.0. 
        Don't forget to factor in the 2/sqrt(pi).  Don't your
        texts supply an approximation for large x?  Computing the
        sum of a sequence of large terms that alternate in sign is
        not a good numerical technique.
        
        Dan
1136.3use a continued fractionALLVAX::ROTHIf you plant ice you'll harvest windMon Oct 09 1989 11:4576
    If you want to estimate the error function or its complement for
    large arguments it's better to use a continued fraction expansion
    since this avoids the loss of accuracy that results from cancellation
    in the series expansion.  In fact, continued fraction expansions
    which match series expansions will usually have much wider regions
    of convergence than the series itself - not surprising since a
    continued fraction can have poles, while no convergent power series
    can.

    Look in AMS55, Handbook of Mathematical Functions for the necessary
    results.

    Here is a simple example program; if you combine it with a series
    subroutine for small X this will give accurate answers.  You can
    put a print statement in the loop to see how the continued fraction
    is converging...

    Sample output -

    1.0      0.1572992070502852    
    1.5      3.3894853524689278E-02
    2.0      4.6777349810472658E-03
    2.5      4.0695201744495897E-04
    3.0      2.2090496998585442E-05
    3.5      7.4309837234141274E-07
    4.0      1.5417257900280020E-08
    4.5      1.9661604415428875E-10
    5.0      1.5374597944280349E-12
    5.5      7.3578479179743983E-15
    6.0      2.1519736712498913E-17

    - Jim

	REAL*8 FUNCTION ERFC_CF(X)
C
C  Evaluate complementary error function via continued fraction expansion
C  for incomplete gamma function
C
C	erfc(x) = gamma(0.5, x^2)/sqrt(pi)
C
C  Convergence is good for x >= 2, useable for 1 < x < 2
C
C  Ref AMS55 chapters on gamma and error functions
C
	IMPLICIT REAL*8 (A-H, O-Z)

	PARAMETER (ITMAX = 1000, EPS = 1.0E-16)
	PARAMETER (PI = 3.14159265358979323846264338)

	X2 = X**2

	G_OLD = 0.0
	A0 = 1.0
	A1 = X2
	B0 = 0.0
	B1 = 1.0
	FAC = 1.0

	DO I=1,ITMAX
	  A0 = (A1+A0*(I-0.5))*FAC
	  B0 = (B1+B0*(I-0.5))*FAC
	  A1 = X2*A0+I*FAC*A1
	  B1 = X2*B0+I*FAC*B1
	  IF (A1.NE.0.0) THEN	! Protect against zero divide
	    FAC = 1.0/A1
	    G = B1*FAC
	    IF (DABS((G-G_OLD)/G).LT.EPS) GOTO 100
	    G_OLD = G
	  ENDIF
	ENDDO

100	ERFC_CF = DEXP(-X2)*G*X/DSQRT(PI)

	RETURN

	END
1136.4Look in CACM Collected AlgorithmsCOOKIE::PBERGHPeter Bergh, DTN 523-3007Mon Oct 09 1989 12:122
    Collected Algorithms from CACM has several routines to calculate this
    integral.  If memory serves, algorithm 209 is one of them.
1136.5Num. Rec. alsoAKQJ10::YARBROUGHI prefer PiMon Oct 09 1989 13:255
Or try Numerical Recipes, pp 163-4 for perhaps more recent routines. 
Problem: if two or more of these give different results for 6<X<10, say, 
how will you decide which is correct?

Lynn Yarbrough 
1136.6tyvm %^)CLADA::KEATINGrocking the bayTue Oct 10 1989 10:3428
	Thanks y'all for your helpful replies.

        Also, please excuse my not-all-that-it-could-be mathematical
        background. 

        Following Dan's reply, I realised that I had only used a value
        of 2/root(pi) from my calculator in the original program. My
        results now agree to 12 places of decimals with Jim's program up
        to values of x of between 2 and 3. I think then that I'll use
        this combination (the Taylor series expansion from my reply .1
        and Jim's continued fraction routine for x > 2). I may be able
        to get access to a set of IMSL routines, but licensing issues
        pervade. 

	Unfortunately I do not have access to Collected Algorithms or
	Numerical Recipes (except a photocopy of Chapter 12!).

        I would still be interested in an answer to Lynn's question: how
        DO I determine which method gives most accurate results? 

        Finally, and on a completely unrelated note, does anybody know
        whether and where there exists a Notes Conference on Chaos? I've
        been doing some light reading about it recently and it's a
        fascinating subject.

	Thanks again for all the help,

	   - Brian. 
1136.7Use the direct approachAKQJ10::YARBROUGHI prefer PiTue Oct 10 1989 12:2011
One way of insuring accuracy is to find a method that computes the 
complementary erf (erfc) directly rather than as 1-erf, which obviously 
leads to roundoff errors. Since both erf and erfc can be computed directly
from Gamma functions, the best approach would seem to be to find a good
Gamma routine with lots of precision. Such routines are available. In this 
regard I cannot recommend Numerical Recipes too highly.

There is a conference called MANDELBROJT - I think that's the right 
spelling - that deals with some chaos issues.

Lynn 
1136.8You may find what you want here ("toms" from netlib)VMSDEV::HALLYBThe Smart Money was on GoliathTue Oct 10 1989 13:24965
Disclaimer. This on-line distribution of algorithms is not an official
service of the Association for Computing Machinery (ACM).  This service is
being provided on a trial basis and may be discontinued in the future.  ACM
is not responsible for any transmission errors in algorithms received by
on-line distribution.  Official ACM versions of algorithms may be obtained 
from:	IMSL, Inc.
	2500 ParkWest Tower One
	2500 CityWest Blvd.
	Houston, TX 77042-3020
	(713) 782-6060

Columns 73-80 have been deleted, trailing blanks removed, and program
text mapped to lower case.  Assembly language programs have sometimes
been deleted and standard machine constant and error handling
routines removed from the individual algorithms and collected
into the libraries "core" and "slatec".  In a few cases, this trimming
was overzealous, and will be repaired when we can get back the original
files. We have tried to incorporate published Remarks; if we missed
something, please contact Eric Grosse, 201-582-5828, research!ehg or
ehg@research.att.com.

The material in this library is copyrighted by the ACM, which grants
general permission to distribute provided the copies are not made for
direct commercial advantage.  For details of the copyright and
dissemination agreement, consult the current issue of TOMS.

Careful! Anything free comes with no guarantee.
*** from netlib, Mon Oct  9 17:11:43 EDT 1989 ***

===== TOMS =====

Netlib has the Collected Algorithms published in TOMS.

To get, say, the Fornberg extraction of power series coefficients
(Algorithm 579), send netlib mail of the form:
    send 579 from toms

Some of the algorithms are enormous.  Please be reasonable when
using 1200-baud connections!

  The following is adapted from Algorithm 620, prepared by J. R. Rice
and R. J Hanson.  More recent items are added to the collection
and the index as soon as they are received from IMSL.  Items 1-492,
indexed in the original have been omitted here since they aren't online.
As of 12 Oct 1988, we have 493-655 and 660-661.

  these are the original references, titles, volume, issue and
page numbers, for acm-calgo algorithms published in acm transactions on
mathematical software (acm-toms).
the references listed below are of the following form:  the
volume and issue number is first.  then the page number of the
primary article about the algorithm.  then the title of the work
is given.  for example, the following text appears for algorithm
493:
 
493  1,2  p. 178
zeroes of a real polynomial
 
this should be interpreted as "acm-calgo algorithm 493 appeared in
acm-toms, vol. 1, number 2.  the reference 'zeroes of a real poly-
nomial' begins on page 178."
 
    the list of references is followed by a line of all dashes
and then by a list of keywords for the algorithms.  these
keywords have been prepared from an understanding of the content
of the algorithm.  the keywords given might not have appeared
with the primary publication.  the keyword omit is given for
algorithms which are considered to be trivial or obsolete for
any but historical purposes.
 
493  1,2  p. 178
zeroes of a real polynomial
494  1,3  p. 261
pdeone, solution of systems of partial differential equations
495  1,3  p. 264
solution of an over determined system of linear equations in the chebysh
norm
496  1,3  p. 271
the lz algorithm to solve the generalized eigenvalue problem for complex
matrices
497  1,4  p. 369
automatic integration of functional differential equations
98  1,4  p. 372
airy functions using chebyshev series approximations
499  2,1  p. 82
an efficient scanning technique
500  2,1  p. 87
minimization of unconstrained multivariate functions
501  2,1  p. 95
fortran translation of algorithm 409, discrete chebyshev curve fit
502  2,1  p. 98
dependence of solution of nonlinear systems on a parameter
503  2,2  p. 196
an automatic program for fredholm integral equations for the second kind
504  2,2  p. 200
gerk:global error estimation for ordinary differential equations
505  2,2  p. 204
a list insertion sort for keys with arbitrary key distribution
506  2,3  p. 275
hqr3 and exchng: fortran subroutines for calculating and ordering the
eigenvalues of a real upper hessenberg matrix
507  2,3  p. 281
procedures for quintic natural spline interpolation
508  2,4  p. 375
matrix bandwidth and profile reduction
509  2,4  p. 509
a hybrid profile reduction algorithm
510  2,4  p. 388
piecewise linear approximations to tabulated data
511  3,1  p. 93
cdc 6600 subroutines ibess and jbess for bessel functions i(sub)(greek
nu)(x),  j(sub)(greek nu)(x)  (greek nu)  .ge.o,  x  .ge.o.
512  3,1  p. 96
a normalized algorithm for the solution of positive definite symmetric
quindiagonal systems of linear equations
513  3,1  p. 104
analysis of in-situ transposition
514  3,2  p. 175
a new method of cubic curve fitting using local data
515  3,2  p. 180
generation of a vector from the lexicographical index
516  3,2  p. 183
an algorithm for obtaining confidence intervals and point estimates
based on ranks in the two-sample location problem
517  3,2  p. 186
a program for computing the conditions numbers of matrix eigenvalues
without computer eigenvector
518  3,3  p. 279
incomplete bessel function io: the von mises distribution
519  3,3  p. 285
three algorithms for computing kolmogorov-smirnov probabilities with
arbitrary boundaries and certification of algorithm 487
520  3,3  p. 295
an automatic revised simplex method for constrained resource network
scheduling
521  3,3  p. 301
repeated integrals of the coerror function
522  3,4  p. 404
esolve, congruence techniques for the exact solution of integer systems
of linear equations
523  3,4  p. 411
convex, a new convex hull algorithm for planar sets
524  4,1  p. 71
mp, a fortran multiple-precision arithmetic package
525  4,1  p. 82
adapt, adaptive smooth curve fitting
526  4,2  p. 160
bivariate interpolation and smooth surface fitting for irregularly
distributed data points
527  4,2  p. 165
a fortran implementation of the generalized marching algorithm
528  4,2  p. 177
framework for a portable library
529  4,2  p. 189
permutations to block triangular form
530  4,3  p. 286
an algorithm for computing the eigensystem of skew-symmetric matrices
and a class of symmetric matrices
531  4,3  p. 290
contour plotting
532  4,4  p. 388
software for roundoff analysis
533  4,4  p. 391
nspiv,  a fortran subroutine for sparse guassian elimination with partia
pivoting
534  4,4  p. 399
stint: stiff (differential equations) integrator
535  4,4  p. 404
the qz algorithm to solve the generalized eigenvalue problem for complex
matrices
536  5,1  p. 108
an efficient one-way enciphering algorithm
537  5,1  p. 112
characteristic values of mathieu's differential equations
538  5,1  p. 118
eigenvectors and eigenvalues of real generalized symmetric matrices by
simultaneous iteration
539  5,3  p. 324
basic linear algebraic subprogram for fortran usage
540  5,3  p. 326
pdecol, general collocation software for partial differential equations
541  5,3  p. 352
efficient fortran subprograms for the solution of separable elliptic
partial differential equations
542  5,4  p. 482
incomplete gamma functions
543  5,4  p. 490
fft9, fast solution of helmholtz-type partial differential equations
544  5,4  p. 494
l2a and l2b, weighted least squares solutions by modified gram-schmidt
with iterative refinement
545  5,4  p. 500
an optimized mass storage fft
546  5,4  p. 88
solveblok
547  6,1  p. 92
fortran routines for discrete cubic spline interpolation and smoothing
548  6,1  p. 104
solution of the assignment problem
549  6,1  p. 112
weierstrass' elliptic functions
550  6,1  p. 121
solid polyhedron measures
551  6,2  p. 228
a fortran subroutine for the l (sub 1) solution of overdetermined system
of linear equations
552  6,2  p. 231
solution of the constrained l (sub 1) linear approximation problem
553  6,2  p. 236
m3rk, an explicit time integrator for semidiscrete parabolic equations
554  6,2  p. 240
brentm, a fortran subroutine for the numerical solution of systems of
nonlinear equations
555  6,2  p. 252
chow-yorke algorithm for fixed points or zeros of c(super 2) maps
556  6,3  p. 420
exponential integrals
557  6,3  p. 429
pagp, a partitioning algorithm for (linear) goal programming problems
558  6,3  p. 430
a program for the multifacility location problem with rectilinear distan
by the minimum-cut approach
559  6,3  p. 432
the stationary point of a quadratic function subject to linear constrain
560  6,3  p. 437
jnf, an algorithm for numerical computation of the jordan normal form of
complex matrix
561  6,3  p. 444
fortran implementation of heap programs for efficient table maintenance
562  6,3  p. 450
shortest path lengths
563  6,4  p. 609
a program for linearly constrained discrete l (sub 1) problems
564  6,4  p. 615
a test problem generator for discrete linear l (sub 1) approximation pro
565  7,1  p. 126
pdetwo/psetm/gearb: solution for systems of two-dimensional nonlinear pa
differential equations
566  7,1  p. 136
fortran subroutines for testing unconstrained optimization software
567  7,1  p. 141
extended range arithmetic and normalized legendre polynomials
568  3,1  p. 162, acm-toplas, april, 1981.
pds-a portable directory system    (omitted from netlib because of size)
569  7,2  p. 223
colsys: collocation software for boundary value ode`s
570  7,2  p. 230
lopsi: a simultaneous iteration algorithm for real matrices
571  7,2  p. 233
statistics for von mises' and fisher's distribution of directions: i(sub
/i(sub 0)(x), i(sub 1.5)(x)/i(sub .5)(x)
572  7,2  p. 239
solution of the helmholtz equation for the dirichlet problem on general
bounded three dimensional regions
573  7,3  p. 369
nl2sol-an adaptive nonlinear least-squares algorithm
574  7,3  p. 384
shape-preserving osculatory quadratic splines
575  7,3  p. 387
permutations for a zero-free diagonal
576  7,3  p. 391
a fortran program for solving ax=b
577  7,3  p. 398
algorithms for incomplete elliptic integrals
578  7,4  p. 537
solution of  real linear equations in a paged virtual store
579  7,4  p. 542
cpsc, complex power series coefficients
580  7,4  p. 548
qrup: a set of fortran routines for updating qr factorizations
581  8,1  p. 84
an improved algorithm for computing the singular value decomposition
582  8,2  p. 190
the gibbs-poole-stockmeyer and gibbs-king algorithms for reordering spar
matrices
583  8,2  p. 195
lsqr: sparse linear equations and least-square problems
584  8,2  p. 210
cubtri-automatic cubature over a triangle
585  8,3  p. 290
a subroutine for the general interpolation and extrapolation problems
586  8,3  p. 302
itpack 2c: a fortran package for solving large sparse linear systems by
adaptive accelerated iterative methods
587  8,3  p. 323
two algorithms for the linearly constrained least squares problem
588  8,4  p. 369
fast hankel transforms using related and lagged convolutions
589  8,4  p. 371
sicedr: a fortran subroutine for improving the accuracy of computed
matrix eigenvalues
590  8,4  p. 376
dsubsp and exchqz: fortran subroutines for computing deflating subspaces
with specified spectrum
591  8,4  p. 383
a comprehensive, matrix-free algorithm for analysis of variance
592  9,1  p. 98
a fortran subroutine for computing the optimal estimate of f(x)
593  9,1  p. 117
a package for the helmholtz equation in nonrectangular planar regions
594  9,1  p. 125
software for relative error analysis
595  9,1  p. 131
an enumerative algorithm for finding hamiltonian circuits in a
directed graph
596  9,2  p. 236
pitcon: a program for a locally parametrized continuation process
597  9,2  p. 242
sequence of modified bessel functions of the first kind
598  9,2  p. 246
an algorithm to compute solvents of the matrix equation ax(super 2)+
bx+c=0
599  9,2  p. 255
sampling from gamma and poisson distributions
600  9,2  p. 258
translation of algorithm 507. procedures for quintic natural spline
interpolation
601  9,3  p. 344
a sparse matrix package-part ii: special cases
602  9,3  p. 355
hurry: an acceleration algorithm for scalar sequences and series
603  9,3  p. 376
colrow and arceco:fortran packages for solving almost block diagonal
linear systems by modified alternate row and column elimination
604  9,3  p. 381
a fortran program for the calculation of an extremal polynomial
605  9,4  p.391
pbasic-a verifier program for ansi minimal basic
606  9,4  p.418
nitpack-an interactive tree package
607  9,4  p.427
text exchange system: a transportable system for management and exchange
programs and other text
608  9,4  p.461
approximate solution of the quadratic assignment problem
609  9,4  p.480
a portable fortran subroutine for the bickley functions ki(sub n)(x)
610  9,4  p.494
a portable fortran subroutine for the derivation of the psi function
611  9,4  p.503
subroutines for unconstrained minimization using a model/trust-region
approach
612 10,1  p.17
triex: integration over a triangle using nonlinear extrapolation
613 10,1  p.108
minimum spanning tree for moderate integer weights
614 10,2  p.140
a fortran subroutine for integration in h(sub p) spaces
615 10,2  p.202
the best subset of parameters in least absolute value regression
616 10,3
fast computation of the hodges-lehman location estimator
617 10,3
dafne - a differential-equations algorithmf for nonlinear equations
618 10,3
fortran subroutines for estimating sparse jacobian matrices
619 10,3
automatic numerical inversion of the laplace transform
620 10,4
references and keywords for collected algorithms of the acm
621 10,4
software with low storage requirements for two-dimensional parabolic differential equations
622 10,4
a suimple macro processor
623 10,4
interpolation on the surface of a sphere
624 10,4
triangulation and interpolation of arbitrarily distributed points in the plane
625 10,4
a two-dimensional domain processor
626 10,4
tricp - a contour plot program for triangular meshes
627 11,1
a fortran subroutine for solving volterra integral equations
628 11,1
an algorithm for constructing canonical bases of polynomial ideals
629 11,2
an integral equation program for laplaces equation in three dimensions
630 11,2
bbvscg - a variable storage algorithm for function minimization
631 11,2
finding a bracketed zero by larkin's mehthod of rational interpolation
632 11,2
a program for the 01 multiple knapsack problem
633 11,2
an algorithm for linear dependency analysis of multivariate data
634 11,3
const and eval: routines for fitting multinomials in a least-squares sense
635 11,3
an algorithm for the solution of systems of complex linear equations in the l-infinity norm with constraints on the unknowns
636 11,4
fortran subroutines for estimating sparse hessian matrices
637 11,4
gencol: collocation on general domains with bicubic hermite polynomials
638 11,4
intcol and hermcol: collocation on rectangular domains wiht bicubic hermite polynomials
639 12,1
to integrate some infinite oscillating tails
640 12,1
efficient calculation of frequency response matrices from state space models
641 12,2
exact solution of general integer systems of linear equations
642 12,2
a fast procedure for calculating minimum corss-validation cubic smoothing splines
643 12,2
fexact: a fortran subroutine for fisher's exact test on unordered rxc contingency tables
644 12,3
bessel functions of a complex argument and nonnegative order
645 12,3
testing programs that compute the generalized inverse of a matrix
646 12,3
pdfind: find a positive definite linear combination of two real symmetric matrices
647 12,4
implementation and relative efficiency of quasirandom sequence generators
648 13,1
nsdtst and stdtst: routines for assessing the performance of initial value solvers
649 13,1
a package for computing trigonometric fourier coefficients based on lyness's algorithm
650 13,2
efficient square root implementaton on the 68000
651 13,3
algorithm hfft - high-order fast-direct solution of the helmholtz equation
652 13,3
hompack: a swuite of codes for globally convergent homotopy algorithms
653 13,3
translation for algorithm 539: pc-blas for 8087, 80287
654 13,3
fortran subroutines for computing the incomplete gamma function ratios and their inverse
655 13,4
iqpack: fortran subroutines for the weights of interpolatory quadratures

656 14,1
extended blas and test results

657 14,1
contours surfaces of a function of three variables

658 14,1
odessa: ODE solver with explicit simultaneous sensitivity analysis

659 14,1
Sobol quasirandom sequence generator

660 14,2
qshep2d: quadratic shepard method for bivariate interpolation of scattered data

661 14,2
qshep3d: quadratic shepard method for trivariate interpolation of scattered data

662 14,2
numerical inversion of Laplace transform by Weeks' method

663 14,2
blas (539) for Cyber 205

664 14,3
solving large banded linear systems using disk

665 14,4
machar: dynamically determine machine parameters

666 14,4
chabis: zeros of nonlinear equations

667 14,4
sigma: stochastic integration global minimization

668 14,4
h2pec: sampling from the hypergeometric distribution

669 15,1
brkf45:  first order, nonstiff initial value problems for ODE

670 15,1
a Runge-Kutta-Nystrom code

671 15,1
farb-e-2d: contour bicubics with area fill

672 15,2
interpolatory quadrature rule, preassigned nodes, general weight function

673 15,2
dynamic huffman coding

674 15,2
condition estimation, one-norm of a matrix

---------------------------------------------------------------------
the keywords are in the following format: algorithm number, algorithm
classification according to the modified share scheme and the
list of keywords (one per line after the first line).
 
 493  c2   polynomial zeros
 494  d3   partial differential equations
           method of lines
           ordinary differential equations
 495  f4   chebyshev solution
           linear system
           linear programming
           simplex method
 496  f2   eigenvalues
           generalized eigenvalue problem
 497  d2   functional differential equations
           numerical integration
           one-step and multistep methods
 498  s17  airy function
           chebyshev series approximation
           chebyshev coefficients
           asymptotic expansion
           taylor expansion
 499  z    pattern recognition
           partial differential equations
           finite differences
           laplace's equation
 500  e4   minimization
           optimization
 501  e2   approximation
           polynomial approximation
           exchange algorithm
           chebyshev approximation
           polynomial approximation
 502  c5   nonlinear equations
           differentiation with respect to a parameter
           one-parameter embedding
           dependence on parameter
 503  d5   linear integral equations
           nystrom method
 504  d    ordinary differential equations
           initial value problems
           global error estimation
           runge-kutta-fehlberg method
 505  s2   sorting
           searching
           linked lists
           data structure
           list operation
 506  f2   eigenvalues
           qr-algorithm
 507  e1   approximation
           interpolation
           spline approximation
           quintic natural spline
 508  f1   bandwidth reduction
           profile reduction
           sparse matrix
 509  f1   bandwidth reduction
           king algorithm
           profile reduction
           sparse matrix
           approximation
 510  e2   piecewise linear function
 511  s18  bessel function of the first kind
           modified function
           airy function
           uniform asymptotic expansion
 512  f4   linear  function
           normalized solution
           periodic quindiagonal
           symmetric positive definite
 513  f1   transposition in place
           matrix transposition
           permutation
 514  e2   interpolation
           cubic splines
           spline approximation
 515  g6   combinations
 516  f2   confidence interval
           illinois method
           point estimate
           regula falsi
           wilcoxon-mann-whitney rank test
 517  f2   eigenvalues
           condition number
 518  s14  incomplete bessel function
           von mises distribution
 519  s14  kolmogorov-smirnov probabilities
 520  h    resource allocation
           linear programming
 521  s15  integral of the coerror function
           miller's recurrence algorithm
 522  f4   symbolic and algebraic manipulation
           linear system
           congruence techniques
 523  z    partitioning
           sorting
 524  a1  multiple precision
           extended precision
           floating arithmetic
           elementary function
           euler's constant
           gamma function
           bessel function
           exponential integral
           logarithmic integral
           bernoulli numbers
           zeta function
 525  e2   spline approximation
           adaptive curve fitting
           hermite interpolation
 526 e1    bivariate interpolation
           piecewise polynomial interpolation
 527  d3   marching algorithms
           block tridiagonal matrix
           elliptic partial differential equations
 528  z    portability
           mathematical software libraries
           error handling
           storage management
           memory allocation
           machine dependencies
 529  f1   symmetric permutations
           block triangular form
           depth first search algorithm
           sparse matrix
 530  f2   eigenvalues
           eigenvectors
           skew-symmetric matrix
           symmetric matrix
           matrix with constant diagonal
 531  j6   contour plotting
 532  z    automatic roundoff analysis
           numerical stability
           numerical linear algebra
 533  f4   sparse matrix
           sparse matrix
           simultaneous linear equations
           partial pivoting algorithms
 534  d2   stiff differential equations
           stiffly stable methods
           composite multistep methods
           cyclic methods
           numerical integration
           ordinary differential equations
           initial value problems
           multistep formulas
           numerical integration program
 535  f2   eigenvalues
           generalized eigenvalue problem
 536  z    one-way security transformation
           password
           encipher
           decipher
           multiprecision integer arithmetic
 537  s22  mathieu's differential equation
           wave equation
           characteristic values
           eigenvalues
           separation constants
           mathieu function
           ordinary mathieu function
           modified mathieu function
           elliptic cylinder function
           hyperbolic cylinder function
 538  f2   eigenvalues
           eigenvectors
           sparse matrix
           diagonable matrix
           simultaneous iteration
 539  f1   linear algebra
           utilities
 540  d3   collocation methods
           partial differential equations
           numerical software
           method of lines
 541  d3   elliptic partial differential equations
           software
           linear systems
 542  s14  computation of incomplete gamma function
           taylor's series
           continued fractions
 543  d3   fast fourier transform
           fast helmholtz solver
           fast poisson solver
 544  f4   covariance matrix
           curve fitting
           iterative refinement
           least squares solution
           linear constraints
           overdetermined system of equations
           regression
           underdetermined system of equations
 545  c6   multidimensional fft
           fast fourier transform
           fft
           mass storage fft
           mass store sorting
           optimal sorting
 546  f4   almost block diagonal systems
           gaussian elimination
           spline approximation
           ordinary differential equations
 547  e3   discrete splines
           discrete cubic splines
           discrete natural splines
           interpolation
           smoothing
 548  h    assignment problem
           hungarian algorithm
 549  s21  weierstrass' elliptic function
 550  z    polyhedron
           graphics
           numerical integration
 551  f4   overdetermined system of linear equations
           discrete linear l(sub 1) approximation
           linear programming
           duel simplex algorithm
           triangular decomposition
 552  f4   constrained l(sub 1) approximation
           linear programming
           simplex method
 553  d3   parabolic partial differential equations
           semidiscretization
           explicit time integrator
 554  c5   nonlinear equations
           numerical solution
           brent's method
 555  c5   fixed point
           zero
           nonlinear systems
           homotopy method
           continuation method
           parameterized nonlinear systems
           zero curve of a homotopy map
           fixed points of nonlinear systems
           zeros of nonlinear systems
 556  s13  exponential integral
           miller algorithm
           confluent hypergeometric function
 557  h    goal program
           multiple objective optimization
           constraint partitioning
           simplex method
 558  h    multifacility
           optimal location
           rectilinear distance
           minimum cut
 559  e4   quadratic programming
           orthogonal decomposition
 560  f2   jordan normal form
           canonical form
           eigenvalues
           eigenvectors
           principal vectors
           block diagonal form
 561  z    heap
           table maintenance
 562  h    shortest path
           shortest route problem
 563  f4   numerical analysis
           overdetermined linear systems
           linear constraints
           discrete l(sub 1) approximation
 564  z    l(sub 1) approximation
           least absolute deviation
           problem generator
           test data
 565  d3   partial differential equations
           method of lines
           finite differences
           ordinary differential equations
 566  c5   performance testing
           systems of nonlinear equations
           nonlinear least squares
           unconstrained minimization
           optimization software
 567  a1   angular momentum
           extended-range arithmetic
           legendre polynomials
           overflow
           underflow
 568  z    file directory system
           unix
           ratfor
 569  d2   ordinary differential equations
           boundary-value problems
           collocation
           boundary spline
           mesh selection
           error estimates
           damped newton's method
           general-purpose code
 570  f2   eigenvalues
           eigenvectors
           simultaneous iteration
           real nonsymmetric matrix
           sparse matrix
 571  s14  direction statistics
           von mises distribution
           fisher distribution
           modified bessel function ratio
           continued fraction
           function inversion
           newton-raphson methods
 572  d3   helmholtz equation
           capacitance matrix
           fast poisson solver
           conjugate gradients
 573  e4   unconstrained optimization
           nonlinear least squares
           nonlinear regression
           quasi-newton methods
           secant methods
 574  e1   polynomial interpolation
           osculation
           shape preserving
           convexity preserving
           monotonicity preserving
           bernstein polynomial
 575  f1   nonsymmetric permutations
           maximum transversal
           maximum assignment
           block triangular form
           sparse matrix
 576  f4   simultaneous linear equations
           gaussian elimination
           new pivoting strategy
 577 s21   elliptic integral
           logarithms
           inverse circular function
           inverse hyperbolic function
           r-function
 578  f4   gaussian elimination
           paged virtual store
 579  d4   numerical differentiation
           taylor series coefficients
           analytic function
 580  f5   matrix factorization
           orthogonalization
 581  f1   singular value decomposition
 582  f4   matrix bandwidth
           matrix profile
           matrix wavefront
           banded matrix
           gibbs-poole-stockmeyer algorithm
           gibbs-king algorithm
 583  g2   analysis of variance
           conjugate-gradient method
           least squares solution
           simultaneous linear equations
           regression
           sparse matrix
 584  d1   quadrature rule
 585  e1   convergence acceleration
           extrapolation
           interpolation
           least squares approximation
           neville-aitken scheme
 586  f4   iterative methods
           numerical software
           sparse matrix
 587  f4   linear least squares solution
           equality constraints
           inequality constraints
           nonnegativity constraints
           inconsistent constraints
           covariance matrix
 588  z    hankel transforms of integer order
           bessel function of the first kind
           convolution integral
           linear digital filters
 589  f2   matrix eigensystems
           iterative method
           eigensystem improvement
 590  f2   generalized eigenvalue
           qz algorithm
 591  g2   linear models
           analysis of variance
           unbalanced data
           missing cells
           estimation
           hypothesis testing
           storage-efficient algorithm
           iterative improvement
           balanced data operators
           algebraic-model specification
           rank computations
           g-inverse solution
 592  e2   optimal estimation
           optimal interpolation
           perfect splines
 593  d3   helmholtz equation
           capacitance matrix
           fast poisson solvers
           conjugate gradients
           software package
 594  z    round-off analysis
           relative errors
           numerical stability
 595  h    hamiltonian circuit
           depth-first search
 596  d2   equilibrium problems
           underdetermined systems of equations
           solution manifolds of parameterized equations
           continuation methods
           local parameterization
           limit point computation
 597  s17  bessel function
 598  f2   matrix equations
           solvent
           newton's method
           qz algorithm
 599  g5   gamma distribution
           poisson distribution
           random numbers
           acceptance-rejection method
           simulation
 600  e2   approximation
           interpolation
           spline approximation
           quintic natural spline
 601  f1   sparse matrix
 602  c6   acceleration of convergence
           transformation of divergent series
           levin's 'u' transform
 603  f4   almost block and diagonal systems
           gaussian elimination
           alternate row and column elimination
           two-point boundary-value problems
 604  c1   remes algorithm
           extremal polynomial
           richardson iteration
 605  z    verifiers
           standard conformance
           basic programming language
 606  z    expert systems
           menu-driven applications
           computer-aided instruction techniques
 607  z    text exchange
           text management
           text organization
           text distribution
           text maintenance
 608  h    assignment problem
           quadratic assignment
           heuristic algorithms
           operations research
 609  s13  exponential integrals
           bickley functions
           integral of bessel function
           bessel function
 610  s14  psi functions
           gamma function
           derivatives of the gamma function
 611  e4   trust regions
           quasi-newton
           bfgs secant update
           finite differences
           reverse communication
 612  d1   quadrature rule
           two dimensional integration
           singular integrands
           epsilon algorithm
 613  h    minimum spanning tree
           shortest connection network
 614  d1   quadrature rule
           optimal quadrature rule
           singular integrands
 615  g2   regression
           least absolute value


% ==== Internet headers and postmarks (see DECWRL::"GATEWAY.DOC")
% Received: by decwrl.dec.com; id AA04761; Mon, 9 Oct 89 16:14:38 -0700
% Date: Mon, 9 Oct 89 17:11 EDT
1136.9ALLVAX::ROTHIf you plant ice you'll harvest windWed Oct 11 1989 11:26111
   Re .7 - I gave a complementary error function since that is what the
   base note requested...

   A check on the accuracy of calculations is to do the computation
   in different ways, comparing results.  For the error function, it's
   probably better to use a non-alternating series expansion because there
   will be no cancellation of terms; for the complementary error function,
   the continued fraction is quite accurate.  Another way is to use
   a series involving Laguerre polynomials (which can be calculated via
   a simple recurrance), which probably will execute faster than the continued
   fraction.  Then you could make one routine which used one of the series
   and finished in exactly the same way.

   The reason for using the recurrance for making the NBS tables you mentioned
   is probably that results used in getting some of the values can be used to
   make others, but this is not necessarily good for a routine that just
   returns a value for a random argument.

   The Recipes routines are generally quite good, but sometimes there
   are minor mistakes - caveat emptor.  For a description of incomplete
   gamma functions, there is a paper in the ACM Trans on Mathematical Software,
   "A Computational Procedure for Imcomplete Gamma Functions", W. Gautchi,
   TOMS Vol 5, No. 4, Dec 1979.

   Here is a comparison of an erf routine with an erfc routine, showing
   the number of terms in the series that had to be summed to machine
   accuracy.  The breakpoint is near x = 2, as recommended earlier.
   You could probably get a smidgen less roundoff error by summing the
   series terms in reverse order, but the results are pretty accurate as is.
   In particular, note that the errors are not biased in any direction -
   some are positive and some are negative. You could do calculations in
   H precision to really double check the answers.

   - Jim

    x      erf(x)	   terms     erfc(x)	    terms  1-(erf+erfc)
  ----  ------------------ ----- ------------------ ----- ------------------
  0.50  0.5204998778130465   13  0.4795001221869535  350  0.000000000000E+00
  1.00  0.8427007929497149   19  0.1572992070502851  101  0.000000000000E+00
  1.50  0.9661051464753107   25  0.0338948535246893   51  0.000000000000E+00
  2.00  0.9953222650189528   31  0.0046777349810473   32 -0.555111512313E-16
  2.50  0.9995930479825550   38  0.0004069520174450   23  0.138777878078E-16
  3.00  0.9999779095030014   45  0.0000220904969986   18  0.555111512313E-16
  3.50  0.9999992569016277   53  0.0000007430983723   15 -0.832667268469E-16
  4.00  0.9999999845827421   61  0.0000000154172579   13 -0.277555756156E-16
  4.50  0.9999999998033839   69  0.0000000001966160   12  0.555111512313E-16
  5.00  0.9999999999984625   78  0.0000000000015375   10  0.138777878078E-16

	REAL*8 FUNCTION ERF(X)
C
C  Evaluate error function via non-alternating series expansion of incomplete
C  Gamma function, based on a Kummer function
C
	IMPLICIT REAL*8 (A-H, O-Z)
	PARAMETER (PI = 3.14159265358979323846264338)
	PARAMETER (ITMAX = 1000)

	X2 = X**2
	AP = 0.5
	S = 2.0
	T = S

	DO I = 1,ITMAX
	  AP = AP+1.0
	  T = T*X2/AP
	  ST = S
	  S = S+T
	  IF (S.EQ.ST) GOTO 100
	ENDDO

100	ERF = DEXP(-X2)*X*S/DSQRT(PI)

	RETURN

	END

	REAL*8 FUNCTION ERFC(X)
C
C  Evaluate complementary error function via series expansion in terms of
C  Laguerre polynomials for incomplete gamma function
C
C	erfc(x) = gamma(0.5, x^2)/sqrt(pi)
C
	IMPLICIT REAL*8 (A-H, O-Z)

	PARAMETER (ITMAX = 1000)
	PARAMETER (PI = 3.14159265358979323846264338)

	XL0 = 0.0
	XL1 = 1.0

	S = 0.0
	P = 1.0

	X2 = X**2

	DO I = 1,ITMAX
	  XL2 = ((I-1.5)*(XL1-XL0)+(I+X2)*XL1)/I
	  ST = S
	  S = S+P/XL1/XL2
	  IF (S.EQ.ST) GOTO 100
	  XL0 = XL1
	  XL1 = XL2
	  P = P*(I-0.5)/(I+1)
	ENDDO

100	ERFC = DEXP(-X2)*S*X/DSQRT(PI)

	RETURN

	END
1136.10That'll do nicely, thank you!CLADA::KEATINGrocking the bayThu Oct 12 1989 11:0541
        Thanks again, Jim ... what you've given is perfect for what I want.
        That is (for the curious), I want to be able to calculate the Bit
        Error Rate (or probability of error) at the output of an optical
        receiver, which is given by 

                       BER = 0.5 * ERFC ( Q / SQRT (2.0) )

	where Q is effectively the signal to noise ratio at the Rx o/p.
	The region of interest is when Q is in the region of 3 to about
	10, as shown by the table below: outside of that range, the
	BER is so great or so small as to render its actual value
	irrelevant.

        I've sent mail to NETLIB looking for the TOMS coerror function
        algorithm, but have had nothing back since Tuesday, which is
        unusually slow for it. Thanks for that pointer, John.

	Lynn, thanks, I found the MANDELBROT conference (on node TURRIS).

	Cheers,
		- Brian.



                       Table of BER vs Q values
		       ------------------------

			 Q           BER         

		       2.000    0.2275013195D-01
		       3.000    0.1349898032D-02
		       4.000    0.3167124183D-04
		       5.000    0.2866515719D-06
		       6.000    0.9865876450D-09
		       7.000    0.1279812544D-11
		       8.000    0.6220960574D-15
		       9.000    0.1128588406D-18
		      10.000    0.7619853024D-23
		      11.000    0.1910659574D-27
		      12.000    0.1776482112D-32
		      13.000    0.0000000000D+00