|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 52992 (0xcf00)
Types: TextFile
Names: »matprocinp«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦09b4e9619⟧ »thcømat«
└─⟦this⟧
*se "*
"pl 297,4,272,0,0"
"lw180""pn 0,1""ps0"
"sb #,6"
H. C. Ørsted Institute"nl"
Chemistry Laboratory III"nl"
Universitetsparken 5"nl"
DK-2100 København Ø"nl14"
"ct"
external algol6 procedures
"nl2"
RC4000 Mathematical Library
"nl7"
'matprocman'
"nl4"
Rene Moss
"nl17"
"qr"1981-02-27"nl2"rapport no.80-1
"pn 2,1""rj"
"ps0"
Content."nl"
--------"nl2"
adapint"nl"
besselik"nl"
besseljy"nl"
beta"nl"
bisec"nl"
chstiffp"nl"
decompose / solve "nl"
detgauss "nl"
diag"nl"
eberlein"nl"
fft"nl"
fftipow"nl"
fftir"nl"
fit "nl"
gamma"nl"
householder "nl"
invertsym"nl"
jacobi "nl"
lobachevski"nl"
matrixmult"nl"
min1a"nl"
min1b"nl"
min2a"nl"
minimum "nl"
minl"nl"
minn2d"nl"
polfit1"nl"
polfitw"nl"
prik"nl"
pzero "nl"
ran, see rand"nl"
rana"nl"
rani"nl"
ranl"nl"
rann, see rana"nl"
ranr"nl"
rand"nl"
rkfifth "nl"
solve, see decompose"nl"
specc"nl"
sperc"nl"
spercre"nl"
spln3"nl"
spln3int"nl"
spln3inv"nl"
spln3smth"nl"
spln3val / spln3dif / spln3dif2"nl"
symdet"nl"
symin"nl"
symsol"nl"
w3j "nl"
w3l "nl"
zero1 "nl"
"qr"1981 februar 27
"rj"
"ps20"
Submitting procedures to the library."nl"
-------------------------------------"nl2"
Users are invited to submit mathematical procedures of general interest
to the library. At least the proceduretext, suitable test-documentation
and reference to a detailed description of the method
should be handed to
"lm40""nl"lektor O.J.Heilmann"nl"
Kemisk Laboratorium III
"lm0""nl2"
who then decide whether the procedure should be included
in the library.
"nl5"
References"nl"
----------"nl2"
The description of the procedures is rather short and follows
the conventions from the ALGOL6 user manual. Details on the
method, precision, core-store and time requirement and the
test-method may be found in the reference, which is
given for each procedure."nl"
The reference can be
"lm40""nl2"
An algorithme from an international journal."nl"
An institute rapport."nl"
A manual."nl"
"lm0""nl"
RC mathematical procedure descriptions may be ordered
from
"lm20""nl2"
A/S Regnecentralen af 1979"nl"
Documentation Department"nl"
Hovedvejen 9"nl"
2600 Glostrup"nl"
Tlf. (02)965366
"lm0""nl2"
The complete set has the document order No: SDO 018.
"nl3"
Formel parameter list"nl"
--------------------- "nl2"
The format of a call of a procedure with name 'procname' can be
obtained by calling the program writestd:
"nl2"######writestd name.procname
"nl2"
"ps0"
"qr"1980-9-9
RCSL NO 55-D48"rj""ct""nl4"
real procedure adapint"rj""nl2"
calculates a definite integral.
"nl2"
call: adapint(a,b,f,x,delta);
"lm 37"
"nl2"
"mt 1, adapint ("
return value, real). The approximation to the integral of f(x)
obtained by the procedure.
"nl2"
"mt 1, a ("
call value, real, integer or long). One of the integration limits.
"nl2"
"mt 1, b ("
call value, real, integer or long). The other integration limit.
"nl2"
"mt 1, f ("
call value, real). The function f(x) given as a real
arithmetric expression in x.
"nl2"
"mt 1, x ("
call and return value, real). The independent variable used
in the expression for f(x). x need not be initialized. At exit
x=sign(e-Ia*delta)*e where e is an estimate of the absolute
error and Ia is the integral of abs(f(x)).
x<=0 indicates a succesfull integration.
"nl2"
"mt 1, eps ("
call value, real, integer or long). The permitted error relative to Ia.
"lm0""nl3"
language: ALGOL; text in adapinttxt"nl2"
Maintenance level: RC
"nl2"References:
RC Mathematical Library: RCSL NO 55-D48.
"ps0"
"qr"1980-10-1
RCSL NO 53-M3"rj""ct""nl4"
procedure besselik"rj""nl2"
calculates the modified Besselfunctions of integer order up to n.
"nl2"
call: besselik(n,x,I,K);
"lm22""nl2"
"mt 1, n ("
call value, integer). Maximum
order of the Bessel functions.
"nl2"
"mt 1, x ("
call value, integer, real or long). The argument, must be >0.
"nl2"
"mt 1, I ("
return value, array(0:n)). The values of the modified Bessel
functions: I(j) = Ij(x), j=0,..,n.
"nl2"
"mt 1, K ("
return value, array(0:n)). The values of the modified Bessel
functions: K(j) = Kj(x), j=0,..,n.
"lm0""rj""nl3"
language: ALGOL; text in besseliktxt"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 53-M3.
"ps0"
"qr"1980-10-1
RCSL N0 53-M2"rj""ct""nl4"
procedure besseljy"rj""nl2"
calculates the Bessel functions of integer order up to n.
"nl2"
call: besseljy(n,x,J,Y);
"lm22""nl2"
"mt 1, n ("
call value, integer). Maximum
order of the Bessel functions.
"nl2"
"mt 1, x ("
call value, integer, real or long). The argument, must be >0.
"nl2"
"mt 1, J ("
return value, array(0:n)). The values of the Bessel
functions: J(j) = Jj(x), j=0,..,n.
"nl2"
"mt 1, Y ("
return value, array(0:n)). The values of the Bessel
functions: Y(j) = Yj(x), j=0,..,n.
"lm0""rj""nl3"
language: ALGOL; text in besseljytxt"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 53-M2.
"ps0"
"qr"1980-10-1
RCSL NO 53-M8"rj""ct""nl4"
real procedure beta"rj""nl2"
calculates the beta function.
"nl2"
call: beta(x,y);
"lm30""nl2"
"mt 1, beta ("
return value, real). The value of the beta function. For values of
the arguments which results in over- or underflow the value
is undefined.
"nl2"
"mt 1, x ("
call value, integer, long or real). Argument.
"nl2"
"mt 1, y ("
call value, integer, long or real). Argument.
"lm0""rj""nl3"
language: ALGOL; text in betatxt"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 53-M8.
"ps0"
"qr"1980-9-9
Ole J. Heilmann"rj""ct""nl4"
boolean procedure bisec"rj""nl2"
finds a zero for f(x) in the closed interval a,b provided
f(a)*f(b)<=0. The method is bisection.
"nl2"
call: bisec(f,x,a,b,eps);
"lm 32"
"nl2"
"mt 1, bisec ("
return value, boolean). true if a zero was found.
"nl2"
"mt 1, f ("
call value, real). f should be a real arithmetric expression in x.
"nl2"
"mt 1, x ("
call and return value, real). x is the argument in the expression
for f(x). The value of x at return specifies the zero.
"nl2"
"mt 1, a ("
call value, real, integer or long). One of the limit of the interval
where a zero is searched for.
"nl2"
"mt 1, b ("
call value, real, integer or long). The other limit.
"nl2"
"mt 1, eps ("
call value, real, integer or long). eps specifies the precision
with which the zero
is determined.
"lm0""nl3"
language: ALGOL; text in bisectxt"nl2"
Maintenance level: user
"ps0"
"qr"1980-9-11
Preben Graae Sørensen"rj""ct""nl2"
integer procedure chstiffp "rj""nl2"
integrate the system dX/dt = F(X,t), X=Xi(t), F=Fi(X,t),
i=1,..,n by Calahans semiimplicit A-stable method.
"nl2"
call: chstiffp(n,fjac,xin,xout);
"nl2""lm40"
"mt 1, chstiffp ("
return value, integer).
chstiffp<0 error condition from
decompose/solve.
"nl2"
"mt 1, n ("
call value, integer). Number of variables."nl2"
"mt 1, fjac ("
call value, procedure). Computation of the functions
F and Jacobian A=dF/dx with following format:
"sj"
procedure fjac(bjac,X,t,F,A);
boolean bjac; real t; array X,F,A;
call: X array(1:n); vector x.
t time t.
return: F array(1:n); values F(x,t).
A array(1:n,1:n); if bjac is true the Jacobian
matrix A else unchanged.
bjac result ?.
"nl""rj"
"mt 1, xin ("
call value, procedure).
Procedure called at the begining of the integration,
initialising the parameters in chstiff.
"sj"
procedure xin(X,tb,dt,te,h,eps);
real tb,dt,te,h,eps; array X;
call: X array(1:n); initial values of X.
tb initial value of t.
dt interval between calls of xout.
te end value of t.
h maximal integrationstep.
eps truncation error.
"rj""nl"
"mt 1, xout ("
call value, procedure). Gives acces to intermediate results.
"sj"
procedure xout(X,t,h,s);
integer s; real t,h; array X;
call: s gives information on the cause of the call
of xout: s=0 t-tb = m*dt
s>0 maximum of X(abs(s))
s<0 minimum of X(abs(s)).
t value of t.
h last value of integration step. h<0 indicate
trouble with automatic steplength selection.
X values of X.
"rj""lm 0""nl2"
language: ALGOL; text in chstifftxt"nl2"
Maintenance level: user.
"nl""sj"
References: P.G.Sørensen, rapport 80-2, Kem.Lab.III
"ps0"
"qr"1980-10-1
RCSL NO 55-D60"rj""ct""nl2"
boolean procedure decompose
procedure solve"rj""nl2"
decompose performs a triangular decomposition of an
arbitrary non-singular matrix of size n*n and a set of n linear
equations can then be solved by solve.
"nl2"
call: decompose(A,P,mode);
"lm42""nl2"
"mt 1, decompose ("
return value, boolean). True if the matrix A is non-singular,
otherwise false.
"nl2"
"mt 1, A ("
call and return value, array with n*n elements). At entry the matrix
to be decomposed, on exit the resulting upper and lower
triangular matrices. In case of a one dimensional array, the elements
must be stored row-wise: A(i*(n-1)+j) = A(i,j).
"nl2"
"mt 1, P ("
return value, integer array with n elements. NB! the size of P defines n).
Contains information about
row-exchanges for use in solve.
"nl2"
"mt 1, mode ("
call value, integer). The parameter governs the precision of the
calculation of the inner-products:
"lm63""nl"
"mt 1, mode=0:" normal floating point mode"nl"
"mt 1, mode=1:" intermediate variables with
45 bits mantisse and 24 bits exponent.
"lm0""rj""nl2"
call: solve(A,P,mode,B);
"lm32""nl2"
"mt 1, A ("
call value, array of n*n elements). The decomposed coefficients-matrix
as proced by decompose.
"nl2"
"mt 1, P ("
call value, integer array with n elements). The informations on row-exchanges
from decompose.
"nl2"
"mt 1, mode ("
call value, integer). as decompose.
"nl2"
"mt 1, B ("
call and return value, array with n elements). On entry the right hand side
of the equations to be solved, on exit the n solutions.
"lm0""nl""sj"
error messages:
real lin. eq. 1 Overflow/underflow outside the
innerproduct procedure.
real lin. eq. 2 Overflow/underflow in the
innerproduct procedure.
decomp 1 lin. eq. 2 A do not have n**2 elements.
solve 2 lin. eq. 2 Wrong content of P.
decomp 3 lin. eq. 2 mode<0 or mode>1.
solve 4 lin. eq. 2 B do not have n elements.
"rj""nl"
language: SLANG"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 55-D60.
"ps0"
"qr"1981-2-20
Preben Graae Sørensen"rj""ct""nl4"
real procedure detgauss"rj""nl2"
solves simultaneously m systems of linear equations with n unknowns
i.e. the matrix equation AX = B where A is a n*n matrix and X,B are
n*m matrices. The determinant of A is calculated."nl2"
It is RECOMMENDED to use decompose/solve to solve linear
equations (especially m=1).
"nl2"
call: detgauss(n,m,M,EXIT);
"lm40""nl2"
"mt 1, detgauss ("
return value, real). The determinant of A.
"nl2"
"mt 1, n ("
call value, integer). The number of unknowns in one set of
linear equations.
"nl2"
"mt 1, m ("
call value, integer). The number of right hand sides, i.e. the
number of linear systems to solve.
"lm50""nl"
"mt 1, m=0"
implies that only the determinant will be computed.
"nl"
"mt 1, m=n"
and B=(the unit matrix) implies that the inverse matrix of A
will be calculated.
"lm40""nl2"
"mt 1, M ("
call and return value, array(1:n,1:n+m)). Contains on entry the
coefficient matrix A and the right hand sides B:
M(i,j)=A(i,j) i=1,..,n j=1,..,n;"nl" M(i,j+n)=B(i,j) i=1,..,n j=1,..,m;
"nl"
on exit the solution matrix X is stored in M:
"nl"
X(i,j)=M(i,j+n) =1,..,n j=1,..,m.
"nl2"
"mt 1, EXIT ("
label). Label to which detgauss jumps when the system is singular.
"lm0""nl3"
language: ALGOL; text in detgausstxt
"nl2"
Maintenance level: user
"nl2"
Reference: GIER system library, Regnecentralen 1963.
"ps0"
"qr"1980-9-8
Knud H. Kjær"rj""ct""nl4"
boolean procedure diag"rj""nl2"
finds the eigenvectors and the corresponding eigenvalues of a
symmetric matrix.
"nl2"
call: diag(n,z,d);
"lm 31"
"nl2"
"mt 1, diag ("
return value, boolean). False if more than 30 iterations
are needed to find one eigenvalue.
"nl2"
"mt 1, n ("
call and return value, integer). The order of the matrix.
"nl2"
"mt 1, d ("
return value, array(1:n)). Contains the eigenvalues in ascending order.
"nl2"
"mt 1, z ("
call and return value, array(1:n,1:n)). At call the coefficient matrix
to be diagonalized; only the lower triangle of the array is
actually needed in the transformations."nl"
At return it contains the normalized eigenvectors column by column of
the call matrix.
"lm0""nl3"
language: ALGOL; text in diagtxta
"nl2"
Maintenance level: user
"nl2"References:"lm30"
Martin, Reinsch & Wilkinson: Householder tridiagonalization of a
symmetric matrix (tred2). Num.Mat. 11 181-195 (1968)."nl"
Bowler, Martin & , Reinsch & Wilkinson: QR and QL algoritms for
symmetric matrices (tql2). Num.Mat. 11 293-306 (1968).
"lm0"
"ps0"
"qr"1980-10-1
RCSL NO 55-D57"rj""ct""nl3"
procedure eberlein
"rj""nl2"
solves the eigenvalue problem for a real matrix by means of a
sequence of Jacobi-like transformations. If this iteration
process does not converge within tmx number of transformations
or converges to a matrix that is not of 1*1 and 2*2
blockdiagonal form, no solution to the eigenproblem is found,
and new is false.
"nl2"
call: eberlein(n,A,T,tmx,new,res);
"lm27""nl2"
"mt 1, n ("
call value, integer). The order of the matrix A.
"nl2"
"mt 1, A ("
call and return value, matrix(1:n,1:n)).
"lm53""nl"
"mt 1, At entry:"
The matrix for which the eigenproblem is to be solved
(or the transformed matrix from a previous unsuccesfull call of eberlein
to be improved if possible by a new sequence of transformations,
new = false).
"lm50""nl"
"mt 1, At exit:"
A contains the result of tmx transformations on
the matrix."nl" If a solution to the eigenproblem is found
"sj"
A(j,j) = real eigenvalues,
A(k,k) +i*A(k,k+1) = complex eigenvalues and
A(k+1,k+1)+i*A(k+1,k) = the conjugate eigenvalues;
unless res=true in which case
A(j,1) = real eigenvalues,
A(k,1) +i*A(k,2) = complex eigenvalues and
A(k+1,1)+i*A(k+1,2) = conjugate eigenvalues.
"rj""lm27""nl"
"mt 1, T ("
call and return value, array(1:n,1:n)). A transformation matrix.
If res is false and tmx>0 at entry then T (from a previous
call of eberlein) is used
for the next transformations. If the eigenproblem is solved T contains the
eigenvectors as columns.
"nl2"
"mt 1, tmx ("
call and return value, integer). At entry abs(tmx) is the maximum
number of transformations performed. If tmx<0 then T is unaltered.
At exit tmx gives the number of transformations performed.
"nl2"
"mt 1, new ("
call and return value, boolean). At entry new = false indicates a
repeated call of eberlein (using old A and T) to improve the
convergence. At exit new is true if convergence occurs.
"nl2"
"mt 1, res ("
call value, boolean). If res is true the eigenvalues, if found,
will be placed in the two first columns of A.
"lm0""rj""nl2"
language: ALGOL; text in eberleintxt"nl2"
Maintenance: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 55-D57.
"ps0"
"qr"1980-9-9
RCSL NO 55-D47"rj""ct""nl4"
procedure fft"rj""nl2"
calculates the fourier-sum of 2**m complex input data of type real.
"nl2"
call: fft(A,B,m,analysis);
"lm 40"
"nl2"
"mt 1, A ("
call and return value, array(0:N-1)).
At call the real part of the N=2**m inputdata; at
return the real part of the N fourier-sum data.
"nl2"
"mt 1, B ("
as A). The imaginary parts corresponding to A.
"nl2"
"mt 1, m ("
call value, integer). N=2**m gives the number
of inputdata. 0<=m<24 but the maximal m depends on the available
space in corestore to the arrays A and B.
"nl2"
"mt 1, analysis ("
call value, integer).
analysis = +1 or -1; determines the sign in the fourier-sum exponents.
"lm0""nl3"
language: ALGOL; text in ffttxt"nl2"
Maintenance level: RC
"nl2"References:
RC Mathematical Library. RCSL NO 55-D47.
"ps0"
"qr"1979-2-23
Rene Moss"rj""ct""nl4"
procedure fftipow"rj""nl2"
finds the power spectrum of 2**m=N real input data of type integer
and returns with N/2+1 values of type real.
"nl2"
call: fftipow(m,A)
"lm 22""nl2"
"mt 1, m ("
call value, integer). N=2**m gives the number of inputdata.
2<=m<=12.
"nl2"
"mt 1, A ("
call and return value, array with N//2+1 elements).
"lm46""nl2"
"mt 1, at call:"
N inputdata r(j),j=1,...,N of type integer, i.e.
"nl"fi:=4*p+2+j; A.fi):=r(j).
"lm51""nl2"
"mt 1, at return:"
N/2+1 power spectrum values at frequency n,
a(n), n=0,...,n/2
"sj"a(0)=A(0), a(2q) =A(q)
a(2q-1)=A(N/2+1-q), q=1,...,N/4
"lm49"
"mt 1, example:"m=3 n=4 N=8
input: ! -!r1!r2!r3!r4!r5!r6!r7!r8! -!
output: ! a0 ! a2 ! a4 ! a3 ! a1 !
"lm0"
"sj"error messages:
index <m> fftipow illegal value for m
"rj""nl2"
language: SLANG; text in fftipowtxt"nl2"
Maintenance level: user
"nl2"
Reference: R.Moss, rapport 54 1975, Kem.Lab. III.
"ps0"
"qr"1979-2-22
Rene Moss"rj""ct""nl4"
procedure fftir"rj""nl2"
finds the fourier-sum of 2**m=N real input data of type integer
and delivers N/2+1 complex numbers also of type integer, i.e.
complex conjugate results are only represented once.
"nl2"
call: fftir(m,A)
"lm 22""nl2"
"mt 1, m ("
call value, integer). N=2**m gives the number of inputdata.
2<=m<=13.
"nl2"
"mt 1, A ("
call and return value, real or integer array with N+2 integer elements).
"lm46""nl2"
"mt 1, at call:"
N inputdata of type integer:"nl"r(j),j=1,...,N, i.e.
"nl"fi:=4*p+2+j; R.fi:=I(h+j):=r(j).
"lm51""nl2"
"mt 1, at return:"
N/2+1 complex outputdata:"nl"
c(n)=a(n)+i*b(n), n=0,...,N/2 (b(0)=b(N/2)=0!) i.e.
"sj"fr:=4*p-2+2*n; a(n):=I(h+n):=R.fr;
fc:=4*p+2*N-2*n; b(n):=I(h+N+1-n):=R.fc;
"lm45"
"mt 1, example:"m=3 n=4 N=8
input: ! -!r1!r2!r3!r4!r5!r6!r7!r8! -!
output: !a0!a1!a2!a3!a4! 0!b3!b2!b1! 0!
"lm0""nl"
"sj"error messages:
index <m> fftir illegal value for m
param fftir wrong type of A"rj""nl2"
language: SLANG; text in fftitxt"nl2"
Maintenance level: user
"nl2"
Reference: R.Moss, rapport 54 1975, Kem.Lab.III.
"ps0"
"qr"1981-2-20
RCSL NO 31-D129"rj""ct""nl4"
procedure fit"rj""nl2"
computes the coefficients of a weighted least-square
polynomial-approximation to a set of points. The degree of
the polynomial is selected within specified limits.
"nl2"
The procedure IS NOT RECOMMENDED; refer to the polfit proceedures.
"ps0"
"qr"1980-10-1
RCSL NO 55-D58"rj""ct""nl4"
real procedure gamma
"rj""nl2"
calculates the gamma function.
"nl2"
call: gamma(z);
"lm32""nl2"
"mt 1, gamma ("
return value, real). The value of the gammafunction of argument z.
"nl2"
"mt 1, z ("
call value, integer, long or real). The argument. Values
outside -301<z<301 will give overflow/underflow, and gamma is undefined.
"lm0""rj""nl3"
language: ALGOL; text in gammatxt"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 55-D58.
"ps0"
"qr"1981-2-20
RCSL NO 53-M7"rj""ct""nl4"
procedure householder"rj""nl2"
calculates eigenvalues and if wanted the corresponding eigenvectors
of a real symmetric matrix M of order n.
"nl2"
call: householder(n,l,h,A,EV,X,eps);
"lm27""nl2"
"mt 1, n ("
call value, integer). The order of the matrix M.
"nl2"
"mt 1, l ("
call value, integer). 1<=l<=n, the eigenvalues from the l'th smallest
and up to the h'th eigenvalue included are calculated.
"nl2"
"mt 1, h ("
call value, integer). l<=h<=n, see l.
"nl2"
"mt 1, A ("
call value, array(1:n*(n+1)/2). A must contain
the elements in the lower triangle of M taken row
be row, i.e. A(i*i(-1)//2+j) = M(i,j), j=1,..,n i=1,..,j.
"nl2"
"mt 1, EV ("
return value, array(l:h)). The calculated eigenvalues in increasing
order.
"nl2"
"mt 1, X ("
return value, array(l:h,1:n+2)). Contains the eigenvectors, if wanted,
corresponding to the eigenvalues. (for each k X(k,n+1)=X(k,n+2)=0).
"nl2"
"mt 1, eps ("
call and return value, real). if eps>0.0 the eigenvectors are calculated.
The absolut value of eps at call determines the precision with which the eigenvalues
EV(k) are calculated. Roughly error(EV(k)) < 4'-10*abs(EV(k))+eps.
At return eps contains an errorbound for any eigenvalue.
"lm0""nl3"
language: ALGOL"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 53-M7.
"ps0"
"qr"1980-10-1
RCSL NO 53-M5"rj""ct""nl4"
boolean procedure invertsym
"rj""nl2"
inverts a symmetrical matrix M of order n. If M is singular a generalized inverse matrix
of M is returned, and the procedure is false.
"nl2"
call: invertsym(n,A);
"lm42""nl2"
"mt 1, invertsym ("
return value, boolean). false if the matrix M is singular, else true.
"nl2"
"mt 1, n ("
call value, integer). The order of the matrix M.
"nl2"
"mt 1, A ("
call and return value, array(1:n*(n+1)/2).
At call the matrix M, A(i*(i-1)//2+j) = M(i,j), j=1,..,n
i=1,..,j i.e. the elements in the lower triangle taken row by row;"nl"
At return A will contain the inverse of M stored in the same way.
"lm0""rj""nl3"
language: ALGOL; text in invertsytxt"nl2"
Maintenance level: RC"nl2"
Reference: RC Mathematical Library. RCSL NO 53-M5.
"ps0"
"qr"1981-2-20
RCSL NO 55-D61"rj""ct""nl4"
real procedure jacobi"rj""nl2"
calculates all eigenvalues and, if wanted, the corresponding
eigenvectors of a symmetric matrix.
"nl2"
The procedure IS NOT RECOMMENDED; refer to householder.
"ps0"
"qr"1980-10-1
Rene Moss"rj""ct""nl4"
real procedure lobachevski
"rj""nl2"
calculates the Lobachevskiy function.
"nl2"
call: lobachevski(x);
"lm 48""nl2"
"mt 1, lobachevski ("
return value, real). The Lobachevskiy function of argument x.
"nl2"
"mt 1, x ("
call value, integer, long or real). The argument. abs(x)<=pi/2.
"lm0""nl3"
"sj"error message:
***call x < -pi/2 or x > pi/2
"rj""nl"
language: ALGOL; text in lobachevtxt"nl2"
Maintenance level: user"nl2"
Reference: R. Moss, rapport 45 oktober 1973, KemLabIII.
"ps0"
"qr"1980-10-1
Niels-Holger Pedersen"rj""ct""nl4"
procedure matrixmult"rj""nl2"
calculates a matrix-product using extended precision internally.
"nl2"
call: matrixmult(p,q,r,A,B,C);
"lm 22""nl2"
"mt 1, p ("
call value, integer). "nl2"
"mt 1, q ("
call value, integer). "nl2"
"mt 1, r ("
call value, integer). "nl2"
"mt 1, A ("
call value, real array with p*q elements).
"nl2"
"mt 1, B ("
call value, real array with q*r elements).
"nl2"
"mt 1, C ("
return value, real array with p*r elements). C contains the matrixproduct
of A and B.
"lm0""nl2""sj"
error messages:
index <tal> matrixmult, wrong number of elements in some arrays.
param matrixmult, p,q or r <= 0.
"rj""nl"
language: SLANG; text in priktxt"nl2"
Maintenance level: user"nl2"
"ps0"
"qr"1980-9-9
Ole J. Heilmann"rj""ct""nl4"
real procedure min1a"rj""nl2"
finds a minimum for a function f(x) of one variable. The method
is a search where x is changed additively.
"nl2"
call: min1a(f,x,step1,eps,minstep,d2);
"lm 37"
"nl2"
"mt 1, min1a ("
return value, real). The value of f(x) at the minimum.
"nl2"
"mt 1, f ("
call value, real). f should be an arithmetric expression in the parameter x
giving the function f(x).
"nl2"
"mt 1, x ("
call and return value, real). x is the argument in the expression
for f(x). At call the value of x where the search should start, at
return the value of x at the minimum or at the stoppoint.
"nl2"
"mt 1, step1 ("
call and return value, real). At call the initial step; at return the final step.
"nl2"
"mt 1, eps ("
call value, real, integer or long).
The search is stopped for step1**2*d2<eps*abs(f). (see the reference)
"nl2"
"mt 1, minstep ("
as f). (see the reference)
call value, real, integer or long). The search is stopped for
abs(step1)<minstep. (see the reference)
"nl2"
"mt 1, d2 ("
return value, real).
The second derivative of f(x) at the minimum or the stoppoint.
"lm0""nl3"
"rj"
language: ALGOL; text in min1atxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, rapport 51 juli 19774, KemLabIII.
"ps0"
"qr"1980-9-9
Ole J. Heilmann"rj""ct""nl4"
real procedure min1b"rj""nl2"
finds a minimum for a function f(x) of one variable. This procedure should only be used
if the sign of the variable at the point where the function
has a minimum is known in advance. The method is a search where x>0 and x is
changed multiplicatively.
"nl2"
call: min1b(f,x,step1,eps,minstep,d2);
"lm 37"
"nl2"
"mt 1, min1b ("
return value, real). The value of f(x) at the minimum.
"nl2"
"mt 1, f ("
call value, real). f should be an real arithmetric expression in the parameter
x giving the function f(x).
"nl2"
"mt 1, x ("
call and return value, real). x is the argument in the expression
for f(x). At call the value of x where the search should start, at
return the value of x at the minimum or at the stoppoint.
"nl2"
"mt 1, step1 ("
call and return value, real). At call the initial steplength factor;
at return the final steplength factor.
"nl2"
"mt 1, eps ("
call value, real, integer or long).
The search is stopped for (step1-1)**2*d2<eps*abs(f). (see the reference)
"nl2"
"mt 1, minstep ("
call value, real, integer or long). The search is stopped for
1/minstep < step1 < minstep, (see the reference).
"nl2"
"mt 1, d2 ("
return value, real).
The second derivative of f(x) at the minimum or the stoppoint.
"lm0""nl3"
"rj"
language: ALGOL; text in min1btxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, rapport 51 juli 1974, KemLabIII.
"ps0"
"qr"1980-9-9
Ole J. Heilmann"rj""ct""nl4"
real procedure min2a"rj""nl2"
finds a minimum for a function f(x,y) of two variables.
"nl2"
call: min2a(f,x,y,ST,minstep,eps,dd);
"lm 37"
"nl2"
"mt 2, min1a ("
return value, real). The value of f(x,y) at the minimum.
"nl2"
"mt 1, f ("
call value, real). f should be a real arithmetric expression in the parameters
x and y giving the function f(x,y).
"nl2"
"mt 1, x ("
call and return value, real). x is an argument in the expression
for f(x,y). At call the value of x where the search should start, at
return the value of x at the minimum or the stoppoint.
"nl2"
"mt 1, y ("
call and return value, real). y is an argument in the expression
for f(x,y). At call the value of y where the search should start; at
return the value of y at the minimum or the stoppoint.
"nl2"
"mt 1, ST ("
call and return value, array(1:2)).
At call ST(1) should contain the initial steplength in the
x-direction and ST(2) the initial steplength in the y-direction;
At return ST(1) and ST(2) contains the final steplengths.
"nl2"
"mt 1, minstep ("
call value, real, integer or long). The search is stopped for
the steplength<minstep. (see the reference)
"nl2"
"mt 1, eps ("
call value, real, integer or long).
The search is stopped for dd<eps*abs(f). (see the reference)
"nl2"
"mt 1, dd ("
call and return value, real).
At entry the number of iterations with the same lowest value of f(x,y),
before one starts to treat equal values of f(x,y) as higher values (i.e.
a non negative integer).
At return dd contains an average second order difference.
"lm0""nl3"
"rj"
language: ALGOL; text in min2atxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, rapport 51 juli 1974, KemLabIII.
"ps0"
"qr"1981-2-20
RCSL NO "rj""ct""nl4"
procedure minimum"rj""nl2"
The procedure IS NOT RECOMMENDED; refer to the other min-procedures.
"ps0"
"qr"1980-9-9
Ole J. Heilmann"rj""ct""nl4"
real procedure minl"rj""nl2"
finds a minimum for a function f(l) of one variable. The method
is a search where l is changed additiv. l is of type long. The search
is stopped when the step is zero.
"nl2"
call: minl(f,l,stepl);
"lm 32"
"nl2"
"mt 1, minl ("
return value, real). The value of f(l) at the minimum.
"nl2"
"mt 1, f ("
call value, real). f should be an arithmetric expression in the parameter l
giving the function f(l).
"nl2"
"mt 1, l ("
call and return value, long). l is the argument in the expression
for f(l). At call the value of l where the search should start, at
return the value of l at the minimum or at the stoppoint.
"nl2"
"mt 1, stepl ("
call value, real, integer or long). The initial steplength.
"lm0""nl3"
"rj"
language: ALGOL; text in minltxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, rapport 51 juli 1974, KemLabIII.
"ps0"
"qr"1980-10-1
Ole J. Heilmann"rj""ct""nl4"
real procedure minn2d"rj""nl2"
finds a minimum for a function of n variables f(X)
using essentially
a gradient search method. The procedure requires the value of
the function as well as of the first and second derivatives.
"nl2"
call: minn2d(n,f,X,eps,minimum);
"lm 37""nl2"
"mt 1, minn2d ("
return value, real). The value of the function f(X) at the stoppoint.
"nl2"
"mt 1, f ("
call and return value, real procedure)."nl"
A procedure to be called with the values of X which then returns
with the value of the function f(X), the values of the first derivatives
and the second derivatives. The procedure must have the format:
"sj"
real procedure f(X,A,B);
array X,A,B;
call: X array(1:n); The vector X.
return: f The value of f(X).
A array(1:n*(n+1)/2); A(i*(i-1)/2+j)=dB(Xi)/dXj
j=1,..,n i=1,..,j.
B array(1:n); B(i)=df(X)/dXi, i=1,..,n.
"rj""nl"
"mt 1, X ("
call and return value, array(1:n)). X contains the n arguments in the
function f(X). At call the values of X where the search should start,
at return the values of X at the minimum or the stoppoint.
"nl2"
"mt 1, eps ("
call and return value, real). The search is stopped when the length of
B < eps and minimum==true. At return eps contains the actuel
length of B."nl"
The search may stop with eps > eps,start. This could be the
case if a minimum point is found (i.e. minimum true)
within machine precision or if the procedure
has given up.
"nl2"
"mt 1, minimum ("
return value, boolean). minimum is true if the matrix A is positive
definit at the stoppoint i.e. the stoppoint is (close to) a mininum.
If minimum is false the stoppoint must be examined to
determine whether it is a minimum, a saddle point, a maximum
or a wrong stoppoint (cfr. eps).
"lm0""nl3"
language: ALGOL; text in minn2dtxt"nl2"
Maintenance level: user"nl2"
"ps0"
"qr"1980-8-9
Ole J. Heilmann"rj""ct""nl4"
real procedure polfit1
"rj""nl2"
From a table of corresponding values of x and y the procedures compute
the coefficients of the polynomial (of specified degree) in x which
minimizes the square sum of the deviations between corresponding
values of y and the polynomial. The procedures also compute standard
deviations on the coefficients.
"lm 0""nl2"
call: polfit1(X,R,n,C,S,r);
"lm 37"
"nl2"
"mt 1, polfit1 ("
return value, real).
The square sum of the deviations divided by (if n>r+1 then
(n-r-1) else 1) (see however error reaction).
"nl2"
"mt 1, X ("
call and return value, array(1:n)). At call the
j'th x value. At return unchanged.
"nl2"
"mt 1, R ("
call and return value, array(1:n)). At call the j'th
y value; at return
the difference between the original values
and the values of the approximating polynomial.
"nl2"
"mt 1, n ("
call value, integer). The number of points.
"nl2"
"mt 1, C ("
return value, array(0:r)).
C(j) will hold the coefficient of
x**j (j=0,...,r) in the approximating polynomial. (see however
error reaction).
"nl2"
"mt 1, S ("
return value, array(0:r)). S(j) will hold
the standard deviations of the corresponding
coefficient assuming that a standard deviation on the y-values is 1 if
n<=r+1 and the sqrt of the value of polfit1 if n>r+1 (see however
error reaction).
"nl2"
"mt 1, r ("
call value, integer). the degree of the approximating polynomial.
"lm0""nl3"
error reaction:"nl"
if r<0 then 'polfit1' has the value -1 and
none of the arrays are touched."nl"
if r>=0 and n<0 then polfit1 has the value -1 and the arrays C and
S are changed to C(0):=...C(r):=0 and S(0):=...S(r):=maxreal. R and X
are not touched."nl"
if r+1>n then C(n):=...C(r):=0 and S(n):=...S(r):=maxreal and the
procedure acts otherwise as if n=r+1.
"nl3"
language: ALGOL; text in polfit1txt"nl2"
Maintenance level: user
"ps0"
"qr"1980-8-9
Ole J. Heilmann"rj""ct""nl4"
real procedure polfitw
"rj""nl2"
From a table of corresponding values of x and y the procedures compute
the coefficients of the polynomial (of specified degree) in x which
minimizes the square sum of the deviations between corresponding
values of y and the polynomial. The procedures also compute standard
deviations on the coefficients. The procedure polfitw includes a
weighting of the polynomial values IN AN UNCOMMON FASHION.
"lm 0""nl2"
call: polfitw(X,R,F,n,C,S,r);
"lm 37"
"nl2"
"mt 1, polfitw ("
return value, real).
The square sum of the deviations divided by (if n>r+1 then
(n-r-1) else 1) (see however error reaction).
"nl2"
"mt 1, X ("
call and return value, array(1:n)). At call the
j'th x value. At return unchanged.
"nl2"
"mt 1, R ("
call and return value, array(1:n)). At call the j'th
y value; at return
the difference between the original values
and the values of the approximating polynomial.
"nl2"
"mt 1, F ("
call value, array(1:n)).
The value of the approximating polynomial for x=X(j) will be
multiplied with F(j) before being compared with R(j).
"nl2"
"mt 1, n ("
call value, integer). The number of points.
"nl2"
"mt 1, C ("
return value, array(0:r)).
C(j) will hold the coefficient of
x**j (j=0,...,r) in the approximating polynomial. (see however
error reaction).
"nl2"
"mt 1, S ("
return value, array(0:r)). S(j) will hold
the standard deviations of the corresponding
coefficient assuming that a standard deviation on the y-values is 1 if
n<=r+1 and the sqrt of the value of polfitw if n>r+1 (see however
error reaction).
"nl2"
"mt 1, r ("
call value, integer). the degree of the approximating polynomial.
"lm0""nl3"
error reaction:"nl"
if r<0 then 'polfitw' has the value -1 and
none of the arrays are touched."nl"
if r>=0 and n<0 then polfitw has the value -1 and the arrays C and
S are changed to C(0):=...C(r):=0 and S(0):=...S(r):=maxreal. R and X
are not touched."nl"
if r+1>n then C(n):=...C(r):=0 and S(n):=...S(r):=maxreal and the
procedure acts otherwise as if n=r+1.
"nl3"
language: ALGOL; text in polfitwtxt"nl2"
Maintenance level: user
"ps0"
"qr"1980-10-1
Niels-Holger Pedersen"rj""ct""nl4"
real procedure prik"rj""nl2"
calculates the scalarproduct of two vectors using
extended precision internally.
"nl2"
call: prik(n,A,B);
"lm 30""nl2"
"mt 1, prik ("
return value, real). The value of the scalarproduct of A and B.
"nl2"
"mt 1, n ("
call value, integer). The number of elements in A and B.
"nl2"
"mt 1, A ("
call value, real array with n elements). "nl2"
"mt 1, B ("
call value, real array with n elements).
"lm0""nl2""sj"
error message:
index <tal> matrixmult, wrong number of elements in some arrays.
param matrixmult, n <= 0.
"rj""nl"
language: SLANG; text in priktxt"nl2"
Maintenance level: user"nl2"
"ps0"
"qr"1981-2-20
RCSL NO 53-M4"rj""ct""nl4"
boolean procedure pzero"rj""nl2"
calculates all roots of a polynomial of order 2,3,4.
"nl2"
The procedure IS NOT RECOMMENDED; refer to sperc.
"ps0"
"qr"1980-10-1
Ole J. Heilmann"rj""ct""nl4"
########procedure rana
integer procedure rani
###long procedure ranl
###real procedure ranr"rj""nl2"
generates pseudo-random numbers with period 2**35 according to
the formula"np10"
x(m+1) = ((x(m)+1)*5**15) mod 2**35, m = 1,2,..."nl"
The generator is controlled by the own long variable rann,
which before each call of one of the procedures is supposed to
contain 2*x(m). After a call of 'rani', 'ranl' or 'ranr', rann will
contain 2*x(m+1) while it after a call of 'rana' will contain
2*x(m+n).
At program start rann is initialized to 2 (i.e. x(1)=1), but another
value can be assigned to rann.
"nl2"
call: rana(n,A);
"lm30""nl2"
"mt 1, n ("
call value, integer). The number of random numbers to be generated.
"nl2"
"mt 1, A ("
return value, any type array with at least n elements).
If A is of type integer, long or real, the first n elements
will be filled with what would have been the results of n
repeated calls of rani, ranl or ranr respectively. If A is of type
boolean, the n first elements will be filled with the first
12 bits of x(m+1),x(m+2),...,x(m+n).
"nl2"
"mt 1,call: rani;""nl2"
"mt 1, rani ("
return value, integer). rani = (x(m+1)shift(-11))extract 24 "nl"
########################i.e. -2**24 <= rani <= 2**24.
"nl2"
"mt 1,call: ranl;""nl2"
"mt 1, ranl ("
return value, long). ranl = x(m+1), "nl"
#####################i.e. 0 <= ranl <= 2*35.
"nl2"
"mt 1,call: ranr;""nl2"
"mt 1, ranr ("
return value, real). ranr = x(m+1)*2**(-35),"nl"
#####################i.e. 0 <= ranr < 1.
"lm0""nl3"
"sj"error messages:
(only from rana).
index <n> n < 0 or A do not have n elements.
param the second parameter is not an array (or zone).
"rj""nl"
language: SLANG; text in ranatxt."nl2"
Maintenance level: user"nl2"
Reference: O.J.Heilmann, rapport 52 juli 1974, KemLabIII.
"ps0"
"qr"1980-10-1
Ole J. Heilmann"rj""ct""nl4"
real procedure rand"rj""nl2"
generates a pseudo-random number. rand is completely equivalent to
ranr except that the random generator is controlled by the own long
variable ran, making the two generators independent.
"nl2"
call: rand;
"lm27""nl2"
"mt 1, rand("
return value, real). rand = x(m+1)*2**(-35), "nl"
i.e. 0<= rand < 1.
"lm0""nl3"
language: SLANG; text in rantxt"nl2"
Maintenance level: user"nl2"
Reference: O.J.Heilmann, rapport 52 juli 1974, KemLabIII.
"ps0"
"qr"1981-2-20
Preben Graae Sørensen"rj""ct""nl4"
procedure rkfifth"rj""nl2"
integrates a system of first order differential equations
dY/dx = F(x,Y) = Fi(x,Y), Y=Yi(x), i=1,..,n from x0,Y0(x0) to b,Y(b) by a
fifth order runge-kutta method.
"nl2"
call: rkfifth(f,x,Y,b,eps,n,fi);
"lm27""nl2"
"mt 1, f ("
procedure), that calculates the derivatives F=Fi(x,Y) i=1,..,n
with the following format:"sj"
procedure f(x,Y,F);
real x; array Y,F;
call: x the independent variable
Y array(1:n); the dependent variables
return: F array(1:n); the derivatives.
"rj"
"nl"
"mt 1, x ("
call and return value, real). The independent variable; at call x=x0.
"nl2"
"mt 1, Y ("
array(1:n)). The dependent variables; at call Y=Y0.
"nl2"
"mt 1, b ("
call value , real, integer or long). The final value of x.
"nl2"
"mt 1, eps ("
call and return value, real). Affects the steplength i.e. the precision of the integration.
"nl2"
"mt 1, n ("
call value, integer). The number of equations and dependent variables.
"nl2"
"mt 1, fi ("
call value, boolean). If fi is false, the initial stplength is set equal
to the final steplength in a former call of the procedure (or 0), otherwise the
steplengh is set equal to b-x.
"lm0""nl3"
language: ALGOL; text in rkfifthtxt
"nl2"
Maintenance level: user
"nl2"
Reference: J.A.Zonnenveld, Automatic numerical integration.
"ps0"
"qr"1980-9-11
Rene Moss"rj""ct""nl4"
procedure specc"rj""nl2"
Solves a Polynomium Equation af degree n with Complex Coefficient.
(for equations having real coefficient see sperc and spercre).
"nl2"
call: specc(AR,n,AI);
"lm 25"
"nl2"
"mt 1, AR ("
call and return value, array(0:n))."nl"
At call: AR(j)=real part of coefficient to z**(n-j) j=0,..,n"nl"
at return: The real part of the zeros = AR(j) j=1,..,n.
"nl2"
"mt 1, n ("
call value, integer). The degree of the equation.
"nl2"
"mt 1, AI ("
call and return value, array(0:n))."nl"
At call: AI(j)=imaginary part of coefficient to z**(n-j) j=0,..,n;
"nl"
at return: The imaginary part of the zeros AI(j) j=1,..,n.
"lm0""nl3"
language: ALGOL; text in sperctxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann & R.Moss, rapport 39, 1973, Kem. Lab. III.
"ps0"
"qr"1980-9-11
Rene Moss"rj""ct""nl4"
procedure sperc"rj""nl2"
Solves a Polynomium Equation af degree n with Real Coefficient
(cfr. spercre).
"nl2"
call: sperc(AR,n,AI);
"lm 25"
"nl2"
"mt 1, AR ("
call and return value, array(0:n))."nl"
At call: AR(j)= coefficient to z**(n-j) j=0,..,n"nl"
at return: The real part of the zeros = AR(j) j=1,..,n.
"nl2"
"mt 1, n ("
call value, integer). The degree of the equation.
"nl2"
"mt 1, AI ("
return value, array(0:n))."nl"
At return: The imaginary part of the zeros AI(j) j=1,..,n.
"lm0""nl3"
language: ALGOL; text in sperctxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann & R.Moss rapport 39, 1973, Kem. Lab. III.
"ps0"
"qr"1980-9-11
Rene Moss"rj""ct""nl4"
procedure spercre"rj""nl2"
Solves a Polynomium Equation af degree n with Real Coefficient
, but the search for REal zeros is favoured, and so it should be
used when some of the zeros are expected or known to be real. (cfr. sperc).
"nl2"
call: spercre(AR,n,AI);
"lm 25"
"nl2"
"mt 1, AR ("
call and return value, array(0:n))."nl"
At call: AR(j)= coefficient to z**(n-j) j=0,..,n"nl"
at return: The real part of the zeros = AR(j) j=1,..,n.
"nl2"
"mt 1, n ("
call value, integer). The degree of the equation.
"nl2"
"mt 1, AI ("
return value, array(0:n))."nl"
At return: The imaginary part of the zeros AI(j) j=1,..,n.
"lm0""nl3"
language: ALGOL; text in sperctxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann & R.Moss, rapport 39, 1973, Kem. Lab. III.
"ps0"
"qr"1978-9-22"nl"Rene Moss"rj""ct""nl4"
procedure spln3"rj""nl2"
determines a cubic spline through a set of points:"np10"
X(i), Y(i) i=1,...,n;
X(i) < X(j) for i<j."nl"
spln3 calculates the 2. derivative of the
spline M(i) at the x-coordinates X(i); X,Y,M define the spline
completely. Different boundary-conditions are possible as specified in
TYPE."nl"The function and the derivatives are evaluated
by spln3val, spln3dif, spln3dif2 and spln3int and the inverse
function is given by spln3inv."nl2"
call: spln3(TYPE,X,Y,M,n)
"nl2""lm31"
"mt 1, TYPE ("
call value, array(0:4)) TYPE(0) selects the type of end-condition,
while the end-condition is specified in TYPE(1:4) as indicated:
"lm36""nl2"
"mt 1,TYPE(0) = 0 : "
general end-conditions:
"nl"
M(1)+TYPE(2)*M(2)########=#TYPE(1)"nl"
#####TYPE(3)*M(n-1)+M(n)#=#TYPE(4)"nl1"
"mt 1, 1 : "
2. derivative = 0 at X(1) and X(n)."nl1"
"mt 1, 2 : "
The endslopes at X(1) and X(n) are
specified in TYPE(1) respectively TYPE(4)."nl1"
"mt 1, 3 : "
2. derivative = 0 (i.e. has support at
x0,xp where
X(1)-x0=3(X(2)-X(1)), x0<X(1)"nl"
######xp-X(n)=3(X(n)-X(n-1)), xp>X(n).
This condition is recommended to give a smooth run-out."nl2"
"lm31""mt 1, X ("
call value, array (1:n)) the x-coordinates in increasing order."nl2"
"mt 1, Y ("
call value, array (1:n)) the corresponding y-coordinates."nl2"
"mt 1, M ("
return value, array (1:n)) the calculated spline parameters, the
second derivatives."nl2"
"mt 1, n ("
call value, integer) the number of points abs(n)>=2. A call of spln3 with
n>=2 initiatize spln3inv and spln3int but this is suppressed by
a call with n<=-2."lm0""nl3"
"sj"error message:
***spln3 n<2 , called with abs(n)<2.
***spln3 illegal X(i) , X(i)>=X(j) for some i<j.
"rj"
"nl1"
Language: ALGOL; text in spln3txt
"nl2"
Maintenance level: user
"ps0"
"qr"1978-9-22
Rene Moss"rj""ct""nl4"
real procedure spln3int"rj""nl2"
finds the integral of a cubic spline from X(1) to x. The integral from
X(1) to X(i) where X(i)<x<X(i+1) is stored in an own variable, and initialized
by a call of spln3."nl2"
call: spln3int(x,X,Y,M,n)
"nl2""lm41"
"mt 1, spln3int ("
return value,real) the integral from X(1) to x."nl2"
"mt 1, x ("
call value, real) the upper limit."nl2"
"mt 1, X ("
call value, array) see spln3"nl2"
"mt 1, Y ("
call value, array) see spln3"nl2"
"mt 1, M ("
call value, array) see spln3"nl2"
"mt 1, n ("
call value, integer) see spln3"lm0""nl3"
"sj"error message:
***spln3int illegal x = <value> , x<X(1) or x>X(n)
"rj""nl1"
Language: ALGOL; text in spln3txt"nl2"
Maintenance level: user
"ps0"
"qr"1978-9-20"nl"Rene Moss"rj""nl4""ct"
boolean procedure spln3inv"rj""nl2"
calculates the inverse function to the cubic spline. Given
the ordinate y the procedure looks for the smallest
inverse value in the interval from x to X(n). The interval
is open to the left if x equals the latest determined
inverse value to y else the interval is closed. If there exist a solution,
this is delivered in x and the procedure is set to true."nl2"
call: spln3inv(x,y,X,Y,M,n)"nl2""lm41"
"mt 1, spln3inv ("
return value, boolean) true, if an inverse value exists in
the interval from x to X(n)."nl2"
"mt 1, x ("
call and return value, real) at call the lower interval,
at return the inverse value if it exists. X(1)<x<X(n)
otherwise an alarm is given."nl2"
"mt 1, y ("
call value, real) the ordinate."nl2"
"mt 1, X ("
call value, array) see spln3."nl2"
"mt 1, Y ("
call value, array) see spln3."nl2"
"mt 1, M ("
call value, array) see spln3."nl2"
"mt 1, n ("
call value, integer) the number of splinepoints. For n<0
the procedure is initialized before it is called with n=abs(n);
n=0 just initialize i.e. the knowledge of the latest
solution is lost. After a normal call of spln3, spln3inv
is initialized."lm0""nl3"
"sj"error message:
***spln3inv illegal x = <value> , x<X(1) or x>X(n)
"rj"
"nl1"
Language: ALGOL; text in spln3txt
"nl2"
Maintenance level: user
"ps0"
"qr"1978-9-22
Rene Moss"rj""nl4""ct"
boolean procedure spln3smth"rj""nl2"
changes the Y-coordinates of a set of points such that
a spline drawn through these new points is smoother,
i.e. the square deviation divided by the variance(SD(i)**2)
#+#the curvature is minimized."nl2"
call: spln3smth(X,YS,SD,n)"nl2"
"lm43"
"mt 1, spln3smth ("
return value, boolean) true if smoothing was possible
and false otherwise. In this case YS is unchanged."nl2"
"mt 1, X ("
call value, array(1:n)). x-coordinates, X(i),
in increasing order."nl2"
"mt 1, Y ("
call and return value, array(1:n)) the y-coordinates,
Y(i)."nl2"
"mt 1, SD ("
call value, array(1:n)) the standard deviation
of the points, SD(i)."nl2"
"mt 1, n ("
call value integer) the number of points. 3<n<29
on RC4000 at the moment. Depends on size of corestore."lm0"
"nl3"
"sj"error message:
***spln3smth n<4 , n<4
***spln3 illegal X(i) , X(i)>=X(j) for some i<j.
"rj""nl1"
Language: ALGOL; text in spln3txt"nl2"
Maintenance level: user"ps0"
"qr"1978-9-22
Rene Moss"rj""nl4""ct"
real procedure spln3val
real procedure spln3dif
real procedure spln3dif2"rj""nl2"
find the value, 1st derivative, 2nd derivative of a cubic
spline at x."nl2"
call: spln3val(x,X,Y,M,n)"nl"
"lm15"
spln3dif(x,X,Y,M,n)"nl"
spln3dif2(x,X,Y,M,n)"nl2"
"lm43"
"mt 1, spln3val ("
return value, real) the value at x."nl2"
"mt 1, spln3dif ("
return value, real) 1st derivative at x."nl2"
"mt 1, spln3dif2 ("
return value, real) 2nd derivative at x."nl2"
"mt 1, x ("
call value, real) X(1)<=x<=X(n)."nl2"
"mt 1, X ("
call value, array) see spln3."nl2"
"mt 1, Y ("
call value, array) see spln3."nl2"
"mt 1, M ("
call value, array) see spln3."nl2"
"mt 1, n ("
call value, integer) see spln3."lm0""nl3"
"sj"error message:
***spln3val illegal x = <value> , x<X(1) or x>X(n)
***spln3dif illegal x = <value> , -
***spln3dif2 illegal x = <value> , -
"rj""nl1"
Language: ALGOL"nl2"
Maintenance level: user
"ps0"
"qr"1980-11-9
Rene Moss"rj""ct""nl4"
boolean procedure symdet"rj""nl2"
calculates the determinant and the cholesky decomposition of a
symmetric, positiv definit matrix, A.
"nl2"
call: symdet(n,a,d1,d2);
"lm 35"
"nl2"
"mt 1, symdet ("
return value, boolean). false, if det(A)=0.
"nl2"
"mt 1, n ("
call value, integer). The order of the matrix.
"nl2"
"mt 1, a ("
call and return value, array(1:n*(n+1)//2))."nl"
At call the matrix A, a(i*(i-1)//2+j) = A(i,j), j=1,..,n i=1,..,j
i.e. the elements in the lower triangle taken row by row,"nl"
at return the cholesky decomposition of A to be used
by the procedure symsol.
"nl2"
"mt 1, d1 ("
return value, real).
"nl2"
"mt 1, d2 ("
return value, integer). Det(A) = d1*2**d2.
"lm0""nl3"
language: ALGOL; text in sympdefmtxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, J.Rotne & R.Moss, rapport 44 1973, KemLabIII.
"ps0"
"qr"1980-11-9
Rene Moss"rj""ct""nl4"
boolean procedure symin"rj""nl2"
calculates the inverse of a
symmetric, positiv definit matrix, A.
"nl2"
call: symin(n,a);
"lm 32"
"nl2"
"mt 1, symin ("
return value, boolean). false, if det(A)=0.
"nl2"
"mt 1, n ("
call value, integer). The order of the matrix.
"nl2"
"mt 1, a ("
call and return value, array(1:n*(n+1)//2))."nl"
At call the matrix A, a(i*(i-1)//2+j) = A(i,j), j=1,..,n i=1,..,j
i.e. the elements in the lower triangle taken row by row,"nl"
at return the lower triangular part of the inverse matrix.
"lm0""nl3"
language: ALGOL; text in sympdefmtxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, J.Rotne & R.Moss, rapport 44 1973, KemLabIII.
"ps0"
"qr"1980-11-9
Rene Moss"rj""ct""nl4"
procedure symsol"rj""nl2"
solves a set of linear equations A*X=B after the cholesky
decomposition of the coefficient
matrix A has been calculated by symdet.
"nl2"
call: symsol(n,a,x,b);
"lm 35"
"nl2"
"mt 1, n ("
call value, integer). The order of the matrix.
"nl2"
"mt 1, a ("
call value, array(1:n*(n+1)//2)).
the cholesky decomposition of A
calculated by the procedure symdet.
"nl2"
"mt 1, x ("
return value, array(1:n)). The solution X to the equations.
"nl2"
"mt 1, b ("
call value, array(1:n)). The vector B in the equations.
"lm0""nl3"
Exampel: if symdet(n,A,d1,d2) then symsol(n,A,x,b);
"nl2"
language: ALGOL; text in sympdefmtxt"nl2"
Maintenance level: user
"nl2"References:
O.J.Heilmann, J.Rotne & R.Moss, rapport 44 1973, KemLabIII.
"ps0"
"qr"1981-2-20
S.E.Harnung"rj""ct""nl4"
real procedure w3j"nl2"
real procedure w3l"sj""nl2"
j1,j2,j3
w3j calculates Wigners 3-j symbol: m1,m2,m3 ,
l1,l2,l3
while w3l calculates the 3-l symbol: t1,t2,t3 ,
corresponding to real spherical solid harmonics.
"rj""nl"
call: w3j(tj1,tj2,tj3,tm1,tm2,tm3);
"lm27""nl2"
"mt 1, w3j ("
return value, real). The value of the function.
"nl2"
"mt 1, tj1 ("
call and return value, integer). tj1:=2*j1 and subject to the condition
j1+j2+j3+1 < 51, beside the usual triangle inequalities.
"nl2"
"mt 1, etc."
"nl3"
"mt 1,call: w3l(t"
l1,tl2,tl3,t1,t2,t3);
"nl2"
"mt 1, tl1 "
as tj1.
"nl2"
"mt 1, etc."
"nl2"
"mt 1, t1 ("
as tj1 ). For t1>=0 a cosinus function is used, otherwise
a sinus function.
"nl2"
"mt 1, t2 "
as t1.
"nl2"
"mt 1, t3"
as t1.
"lm0""nl3"
"sj"error message:
param w3j segm 1 ji<0, ji<abs(mi) i=1,2,3 or
triangle inequality violated.
"rj""nl"
language: SLANG; text in w3jtxt"nl2"
Maintenance level: user
"lm27""nl2"
"mt 1,Reference:"
M. Rotenburg, R.Bivins, N.Metropolis and J.K.Wooten jr.:
The 3-j and 6-j symbols. The M.I.T. Press 1959."nl"
S.E.Harnung and C.E.Schæffer. The 3-l symbol in Structure and
Bonding vol.12 202-295, 1972.
"lm0"
"ps0"
"qr"1981-2-20
RCSL NO 53-M1"rj""ct""nl4"
boolean procedure zero1"rj""nl2"
evaluates a zero of an arbitrary real function.
"nl2"
The procedure IS NOT RECOMMENDED; refer to bisec.
"nl"
"ef"
▶EOF◀