info   MACRO
) File of macros for various items of mathematical computing
) Many use features available only in MacAnova4.11 or later
) Some use features available only in release 2 of MacAnova 4.13 or later
)) Version of 030512
)
) Macros related to matrices
) blockdiag    Construct block diagonal matrix bith given blocks
) kronecker    Compute the Kronecker product of two matrices
) matsqrt      Compute upper or lower triangular or symmetric matrix B
)              such that B' %*% B = A, for given positive definite A
) moorepenrose Compute the Moore-Penrose inverse of a matrix
) qrdcomp      Compute the QR decomposition of a matrix
)
) Macros for working with fully complex forms of complex matrices ) A
) and B
) cdiag        diag(A)
) cmatmultc    matrix product A %*% B
) ctranspose   A'
) cjtranspose  conj(A)'
) ceigen       eigenvalues and eigenvectors of A
) ctrace       trace(A)
) csolve       inverse of A
) csubscr      simulated A[], A[i], A[i,], or A[i,j]
)
) Macros related to optimization
) minimizer    Optimization by choice of several quasi-Newton optimizers:
)              Broyden-Fletcher-Shanno, Davidon-Fletcher-Powell,
)              and Broyden methods
) bfs          Optimization by Broyden-Fletcher-Shanno method using
)              golden section linear search; uses macro minimizer()
) dfp          Optimization by Davidon-Fletcher-Powell method using
)              golden section linear search; uses macro minimizer()
) broyden      Optimization by a method due to Broyden with no linear
)              search; uses macro minimizer()
) neldermead   Optimization by Nelder-Mead simplex direct search
)              method with optional quadratic "polish"
) levmar       Levenberg-Marquart non-linear least squares macro
)              based on a Fortran program of K. M. Brown
) _cgrad       Used by levmar() to compute gradient, Jacobian, Hessian
) _lmout       Used by levmar() to print out info on each iteration
)
) Macros related to polynomials and series in powers of x
) chebcoefs    Compute the coefficients of the expansion of a polynomial
)              in Chebysev polynomials
) economize    "Economize" a power series expansion on [-1, 1]
) invchebcoefs Compute the coefficients powers of x in a finite
)              series involving Chebysev polynomials
) invertseries Macro to find coefficients of y^j of solution to
)              y = sum(a[i]*x^i,i=1,n)
) orthopoly    Compute standard orthogonal polynomials by recursion
)
) Macros related to special sequences and functions
) binom        Compute binomial coefficients
) continfrac   Compute continued fraction b0+a1/(b1+a2/(b2+a3/(b3+a4/... )))
) factorial    Compute x!
) i0           Compute modified Bessel function of the first kind I0(x)
) i1           Compute modified Bessel function of the first kind I1(x)
)
) Other math related macros
) factors      Macro to compute prime factors of each element in a
)              vector of positive integers
) printfactors Macro to print output from factors
) partitions   Macro to compute partitions of integers n
)
) Other macros
) mathhelp     print help and usage on macros in this file
)
)) 001213 added subtopics
)) 011112 Minor changes to match help
)) 011116 dfp(), bfs() and broyden() recognize keyword 'minit',
))   and optionally print or record partial results under the control of
))   keywords 'printwhen' and 'recordwhen'
)) 011118 modified golden section search in dfp(), bfs() and
))   replaced keyword 'epsilon' by 'criteria' to implement  different
))   convergence criteria in dfp(), bfs() and broyden()
)) 011203 neldermead() now has optional argument to specify data or
))   other fixed quantities; also some printed output has been modified
))   and at least one bug fixed.
)) 011210 added minimizer() which combines updates formerly in dfp(),
))   bfs() and broyden().  Replaced dfp(), bfs() and broyden() by
))   versions that use minimizer()
)) 011216 cleaned up quadratic polish in neldermead() and made other
))   changes
)) 011217 added keyword 'epsilon' to moorepenrose()
)) 030318 added macros cdiag(), ceigen(), cmatmultc(), cjtranspose(),
))        csolve(), csubscr(), ctrace(), cjtranspose()
))        for working with complex matrices
)) 030319 Updated cmatmultc() and deleted cmatmultcj()
)) 030320 Added levmar() and supporting macros (moved from arima.mac)
))        rewrote help for levmar()
)) 030401 Modified help for levmar()
)) 030814 added subtopic titles
%info%

matsqrt        MACRO DOLLARS
) Macro to compute the matrix square root of a positive semi-
) definite matrix.
) usage:
)  b <- matsqrt(a)
)  b <- matsqrt(a, lower:T)
)  b <- matsqrt(a, symmetric:T)
)   a        REAL square positive semi-definite matrix with no MISSING
)            values
)   b        REAL square matrix satisfying a = b' %*% b = a
)            With no keywords, b is upper triangular with non-negative
)            diagonal, identicial to cholesky(b)
)            With lower:T, b is lower triangular with non-negative
)            diagonal elements
)            With symmetric:T, b is symmetric positive semidefinite
)            satisfying a = b %*% b
))Version 000824
#b <- $S(a [,symmetric:T or lower:T]) a positive semidefinite symmetric
@a <- argvalue($1,"argument 1","square real nonmissing")
@symmetric <- keyvalue($K,"symmet*","TF",default:F)
@lower <- keyvalue($K,"lower","TF",default:F)
if(@symmetric && @lower){
	error("you can't use 'lower:T' with 'symmetric:T'")
}
if (delete(@symmetric,return:T)){
	@eigs <- eigen(delete(@a, return:T))
	if (min(@eigs[1]) < 0){
		error("argument 1 not positive definite")
	}
	@sqrt <- @eigs$vectors %*% dmat(sqrt(@eigs$values)) %C% @eigs$vectors
	delete(@eigs)
}else{
	if (@lower){
		@J <- run(nrows(@a),1)
		@a <- @a[@J,@J]
	}
	@sqrt <- cholesky(delete(@a,return:T))
	if (@lower){
		@sqrt <- @sqrt[@J,@J]
		delete(@J)
	}
}
delete(@sqrt,return:T)
%matsqrt%

moorepenrose  macro  dollars
) Macro to compute the Moore-Penrose inverse of a matrix
) Usage:
)   b <- moorepenrose(a [, epsilon:eps), a an m by n REAL matrix
)     a         m by n REAL matrix with no MISSING values
)     eps       small positive scalar [2^-51 = 4.4409e-16]
)
)     b         n by m REAL matrix, the Moore-Penrose inverse
)               of a.  When m >= n and a is full rank
)               b is solve(a' %*% a, a')
)
)   b is computed from the singular value decomposition, setting to 0
)   all singular values < eps*max(singular values)
) C. Bingham, kb@stat.umn.edu
))011217 revised to use single matrix multiplication instead of
))       doing sum of outer products in a loop; added keyword 'epsilon'
# b <- $S(a [, epsilon:eps]), REAL matrix a with no MISSING values, scalar eps > 0
if ($v != 1 || $k > 1) {
	error("usage: b <- $S(a [,epsilon:eps), REAL matrix a with no missing values")
}
@a <- matrix(argvalue($01,"argument","nonmissing real matrix"))
@eps <- keyvalue($K, "eps*", "positive number",default:2^-51)
if (isscalar(@a)) {
	@b <- if (@a == 0) {
		0
	} else {
		1/@a
	}
	delete(@a)
} else {
	@nr <- nrows(@a)
	@nc <- ncols(@a)
	@a <- svd(@a,all:T)
	@values <- @a$values
	@left <- @a$leftvectors
	@right <- @a$rightvectors
	delete(@a)

	if (@values[1] > 0) {
		@J <- @values/@values[1] > @eps
		@b <- (@right[,@J]/@values[@J]') %C% @left[,@J]
		delete(@J)
	} else {
		@b <- matrix(rep(0,@nr*@nc),@nc)
	}
	delete(@values,@left,@right,@nr,@nc, @eps)
}
delete(@b,return:T)
%moorepenrose%

blockdmat  MACRO DOLLARS
) Macro to return block diagonal matrix
) Usage
)  d <- blockdmat(a1,a2,...,ak)
)   a1,...,ak    matrices all of the same type, REAL, LOGICAL, CHARACTER
)   d            (nrows(a1)+...+nrows(ak)) by (ncols(a1)+...+ncols(ak))
)                matrix, with a1, a2, ...,ak as blocks on the diagonal
) Version of 000824
# usage: $S(a1,a2,...,ak), all arguments matrices of the same type
if ($k > 0 || $v == 0){
	error("usage: $S(a1,a2,...,ak), all arguments matrices of the same type",\
		macroname:F)
}
@args <- structure($0)
for(@i, 1, $v){
	@arg <- @args[@i]
	if (!ismatrix(@arg)){
		error(paste("argument",@i,"is not a matrix"))
	}
	@type <- if (isreal(@arg)){
		1
	} elseif (islogic(@arg)){
		2
	} else {
		3
	}
	if (@i == 1) {
		@type1 <- @type
	} elseif (@type != @type1) {
		error(paste("argument",@i,"has different type from argument 1"))
	}
}
@M <- sum(vector(nrows(@args)))
@N <- sum(vector(ncols(@args)))
@c <- structure(0,F,"")[@type]

@result <- matrix(rep(@c,@N*@M),@M)
delete(@c,@M,@N)
@placer <- @placec <- 0
for(@i, 1, $v){
	@ai <- @args[@i]
	@nri <- nrows(@ai)
	@nci <- ncols(@ai)
	@result[run(@nri) + @placer,run(@nci) + @placec] <- @ai
	@placer <-+ @nri
	@placec <-+ @nci
}
delete(@args,@i,@placer,@placec,@ai,@nri,@nci)
delete(@result, return:T)
%blockdmat%

kronecker     MACRO DOLLARS
) Macro to compute kronecker product of two matrices
) Usage:
)  r <- kronecker(a,b)
)   a, b       REAL matrices of with no MISSING values
)   r          REAL ma*mb by na*nb, where na = nrows(a),
)              ma = ncols(a), mb = nrows(b), nb = ncols(b).
)              r consists of ma*na blocks of the form a[i,j]*b,
)              i = 1,..,ma, j = 1,...,na
))000220 stripped $$, use argvalue()
# r <- $S(a,b), a and b REAL matrices with no missing values
if($v != 2){ #r <- kronecker(a,b)
	error("$S must have two matrix arguments")
}
@a <- matrix(argvalue($1, "argument 1", "real matrix nonmissing"))
@b <- matrix(argvalue($2, "argument 2", "real matrix nonmissing"))

@ma <- nrows(@a)
@na <- ncols(@a)
@mb <- nrows(@b)
@nb <- ncols(@b)
@tmp <- matrix(rep(0,@ma*length(@b)),@mb)
@result <- matrix(rep(0,length(@a)*length(@b)),@ma*@mb)
@I <- run(@mb)
for(@i,run(@ma)){
	@J <- run(@nb)
	for(@j,run(@na)){
		@result[@I,@J] <- @a[@i,@j]*@b
		@J <-+ @nb
	}
	@I <-+ @mb
}
delete(@a,@b,@ma,@na,@mb,@nb,@i,@j,@I,@J,@tmp)
delete(@result,return:T)
%kronecker%

qrdcomp       MACRO DOLLARS
) Front end for qr() to produce both Q and R
) Usage:
)  result <- qrdcomp(x)
)  result <- qrdcomp(qr(x))
)  result <- qrdcomp(x, pivot:T)
)  result <- qrdcomp(qr(x,pivot:T))
)   x        REAL matrix with no MISSING values
)  Without 'pivot:T'
)   result   structure(q:Q, r:R) such that x = Q %*% R, computed
)            without pivoting
)  With 'pivot:T'
)   result   structure(q:Q, r:R, pivot:Pivots) such that
)            x[,Pivots] = Q %*% R, computed with pivoting
)            without pivoting
) The columns of Q are orthonormal and R is upper triangular
)
)) Version of 960918 allowed output from qr as argument
)) Version of 961001 bug fix (now works with square matrix)
)) Version of 981028 PEF (now works with rows < cols & permutes cols of r
)) when pivot is T, so that q %*% r is x)
)) Version of 981117 (permuting cols of r when pivot is T disabled)
)) Version of 000411 stripped $$, reorganized to save memory
#$S(x,[[pivot:] T (default) or F]) or $S(qr(x [,pivot]))
@qr <- $1
if (isstruc(@qr)){
	@names <- compnames(@qr)
	@ok <- if (ncomps(@qr) < 2) {F} else {
		@names[1] == "qr" && @names[2] == "qraux" }
	if (!@ok){
		error("argument is a structure that was not computed by qr()")
	}
	@pvt <- if (ncomps(@qr) < 3){F}else{@names[3] == "pivot"}
}else{
	@qr <- argvalue(@qr,"argument","real nonmissing matrix")
	@pvt <- if($v > 1){
		argvalue($02,"pivot","TF")
	}else{
		keyvalue($K,"pivot","TF",default:F)
	}
	@qr <- qr(matrix(@qr), @pvt)
}

 #extract components of @qr and then delete it
@qrqr <- @qr$qr
if(@pvt){
	@pivot <- @qr$pivot
}
@p <- min(dim(@qrqr))
@prows <- run(@p)
@qraux <- @qr$qraux[@prows]
delete(@qr)

@n <- nrows(@qrqr)
@r <- triupper(@qrqr[@prows,]) # first @p rows
@v <- trilower(delete(@qrqr,return:T))

 # set diagonal of @v to @qraux
@v[hconcat(@prows,@prows)] <- @qraux
delete(@prows)

 # compute q
@q <- padto(dmat(@p,1),@n)
for(@i,@p,1){
	@u <- @v[,@i]
	@q <-- @u %*% ((@u %c% @q)/@qraux[@i]);;
}
delete(@u,@v,@p,@n, @i,@qraux)

@qr <- if (delete(@pvt,return:T)){
	structure(r:delete(@r,return:T), q:delete(@q,return:T),\
		pivot:delete(@pivot,return:T))
}else{
	structure(r:delete(@r,return:T), q:delete(@q,return:T))
}
delete(@qr,return:T)
%qrdcomp%

minimizer  macro dollars
) Macro to use one of Broyden-Fletcher-Shanno, Davidon-Fletcher-Powell
) or Broyden method to minimize a function
) Based on Dahlquist and Bjorck, Numerical methods, Prentice Hall, 1974
) p. 441-444
) usage:
)   minimizer(x0, fun [,params] [,method:M] [,h:invhes] [,goldsteps:m]\
)            [,criteria:vector(nsigx,nsigfun,dg)]\
)            [, maxit:maxit] [,minit:minit]\
)            [,printwhen:d1] [,recordwhen:d2]))
)    x0          REAL length k vector of starting values
)    fun         macro to evaluate function (f <- fun(x,0 [,params]))
)                and gradient (g <- fun(x,1 [,params]))
)    params      optional variable with information for fun
)    M           CHARACTER scalar, one of "bfs", "dfp" or "broyden"
)                ["bfs"]
)    invhes      symmetric k by k REAL matrix, a starting approximation
)                to the inverse of the Hessian matrix [dmat(k,1)]
)    m           positive integer = number of cycles using golden mean
)                line search, default [5]; ignored for M = "broyden"
)    maxit > 0   integer = maximum number of iterations [30]
)    minit >= 0  integer = minimum number of iterations [0]
)    nsigx       integer; desired significant digits in coefficients [5]
)    nsigfun     integer; desired significant digits in minimum value [8]
)    dg          REAL convergence criterion for ||gradient|| [-1]
)    d1 >= 0     integer [0].  When d1 > 0, at iterations d1, 2*d1, ...
)                partial results are printed
)    d2 >= 0     integer [0].  When d2 > 0, at iterations d2, 2*d2, ...,
)                partial results are recorded in side effect variable
)                BFSRECORD, DJPRECORD or BROYDNRECORD in the form
)                = structure(xvals, funvals, gradient) where
)                xvals[i,], funvals[i], gradient[i,] are values of x,
)                f(x) and gradient(x) at iteration i*d2.
) Return value:
)   structure(x:minimum, f:minimizedValue, gradient:gradient,\
)      h:invhessian, iterations:number_of_interations,status:critno)
) C. Bingham (kb@stat.umn.edu)
)) 011112 Component for interations is now 'iterations' to match help
)) 011116 added keywords minit, printwhen, recordwhen
)) 011118 cleaned up golden section search; removed keyword 'epsilon';
))        added keyword 'criteria'
)) Version 011125, based on bfs(), dfp() and broyden() which have now
))        been rewritten to use minimizer()
)) 031215 status is now -2 when MISSING value found during golden search
#$S(start, fun [, params] [,method:M] [,h:invhes] [, goldsteps:m] \
#   [,criteria:vector(nsigx,nsigfun,dg)] [, maxit:maxit] [, minit:minit] \
#   [,printwhen:d1] [,recordwhen:d2])
if ($v < 2 || $v > 3) {
	error("usage:$S(start, fun [, params] [,method:M] [,h:invhes] \\
    [, goldsteps:m] [,criteria:vector(nsigx,nsigfun,dg)]  \\
    [, maxit:maxit] [, minit:minit] [,printwhen:d1] [,recordwhen:d2])",\
	macroname:F)
}	
@x <- argvalue($1, "starting values", "nonmissing real vector")
@fun <- argvalue($2, "$2", "macro")
 #@fun should have calling sequence @fun(x,0[,params]) or @fun(x,1[,params])
@params <- if ($v == 3) {
	argvalue($03,"data or parameters")
} else {
	NULL
}

@n <- length(@x)
@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(NoTaKeY:NULL)
}

@h <- matrix(keyvalue(@keys,"h","real nonmissing square",default:dmat(@n,1)))
if (nrows(@h) != @n) {
	error(paste("value for 'h' not",@n,"by",@n))
}
@method <- keyvalue(@keys,"meth*","string",default:"bfs")
@gold <- keyvalue(@keys,"gold*","count")
@maxit <- keyvalue(@keys,"maxit*","positive count",default:30)
@minit <- keyvalue(@keys,"minit*","nonneg count",default:0)
@print <- keyvalue(@keys,"printwhen","count",default:0)
@record <- keyvalue(@keys,"recordwhen","count",default:0)
@silent <- keyvalue(@keys,"silent","TF",default:F)
@crit <- keyvalue(@keys,"crit*","nonmissing real vector",\
		default:vector(5, 8, -1))

@imeth <- match(@method,vector("bfs","dfp","broyden"),0)
if (@imeth == 0) {
	error(paste("Unrecognized minimizing method",@method))
}
if (isnull(@gold)) {
	@gold <- (@imeth != 3) * 5
} elseif (@imeth == 3) {
	if (@gold > 0 && !@silent) {
		print(paste("WARNING: goldstep > 0 ignored with method:\"",@method,\
					"\"",sep:""))
	}
	@gold <- 0
}

if (max(@crit) <= 0){
	error("You must have at least one positive criterion value")
}
if (isscalar(@crit)){
	@crit <- vector(@crit,8, -1)
} elseif (length(@crit) == 2){
	@crit <- vector(@crit, -1)
}
if (!isreal(@crit[run(2)],integer:T)){
	error("crit[1] and crit[2] must be integers")
}

@relcoef <- (@crit[1] > 0) * 10^-@crit[1]

@relcrit <- (@crit[2] > 0) * 10^-@crit[2]

@dg <- @crit[3]
delete(@keys,@crit)


if (@record > 0) {
	# set up space for recorded info
	@funvals <- rep(0,ceiling(@maxit/@record))
	@xvals <- @gradients <- @funvals + @x' # create matrices
	@nrec <- 0
}

if (anymissing(@fun(@x,-1,@params))){ # initialize
	error("problem initializing $2()")
}
@f <- @fun(@x, 0,@params)
@g <- if (!ismissing(@f)) {
	@fun(@x, 1,@params)
} else {
	NULL
}
if (anymissing(vector(@f,@g))) {
	error("function or gradient not defined at starting values")
}
if (@gold > 0){
	@r <- (sqrt(5) - 1)/2
}
@critno <- 0
for(@iter, 1, @maxit) {
	@oldf <- @f
	@oldg <- @g
	@d <- vector(-@h %*% @g)
	if (@gold == 0){
		@f <- @fun(@x + @d, 0, @params)
		@lambda <- 1
	} else {
		# do linear search to find @lambda
		@l <- 2.0
		@x2 <- @r*@l
		@x1 <- @r*@x2
		@left <- 0
		@phi1 <- @fun(@x + @x1*@d, 0,@params)
		if (!ismissing(@phi1)) {
			@phi2 <- @fun(@x + @x2*@d, 0,@params)
		}
		for(@j,1,@gold) {
			if (anymissing(vector(@phi1,@phi2))) {
				break
			}
			@l <-* @r
			if (@phi1 < @phi2) { #discard upper segment
				@x2 <- @x1
				@phi2 <- @phi1
				@x1 <- @left + @r^2*@l
				@phi1 <- @fun(@x + @x1*@d, 0,@params)
			} else { #phi1 >= phi2; discard lower segment
				@left <- @x1
				@x1 <- @x2
				@phi1 <- @phi2
				@x2 <- @left + @r*@l
				@phi2 <- @fun(@x + @x2*@d, 0,@params)
			}
		}
		
		if (anymissing(vector(@phi1,@phi2))) {
			@critno <- -2
			@lambda <- 1
			@f <- ?
			break
		}
		if (@phi1 < @phi2) {
			@lambda <- @x1
			@f <- @phi1
		} else {
			@lambda <- @x2
			@f <- @phi2
		}
	}		
	@delta <- @lambda*@d # step for @x
	@x <-+ @delta
	@g <- @fun(@x, 1,@params)

	if (alltrue(@print > 0, @iter %% @print == 0)) {
		print(paste(@iter,": f(",paste(@x,sep:","),") = ",@f,sep:""))
		print(paste("gradient =",paste(@g,sep:",")))
	}
	if (alltrue(@record > 0, @iter %% @record == 0)) {
		@nrec <-+ 1
		@funvals[@nrec] <- @f
		@xvals[@nrec,] <- @x
		@gradients[@nrec,] <- @g
	}
	if (anymissing(vector(@f,@g))) {
		@critno <- -1
		break
	}
	@gamma <- @g - @oldg
	@hgam <- @h %*% @gamma
	if (@imeth == 1) { # bfs
		@denom1 <- sum(@delta*@gamma)
		if (@denom1 != 0) {
			@delhgam <- @delta * @hgam'
			@h <-+ ((@delta * @delta')*(1 + sum(@gamma*@hgam)/@denom1) -\
				(@delhgam + @delhgam'))/@denom1
		}
	} elseif (@imeth == 2) { # dfp
		@denom1 <- sum(@delta * @gamma)
		if (@denom1 != 0) {
			@h <-+ (@delta * @delta')/@denom1
		}
		@denom1 <- sum(@hgam * @gamma)
		if (@denom1 != 0) {
			@h <-- (@hgam * @hgam')/@denom1
		}
	} else { # broyden
		@delta <-- @hgam
		@denom1 <- sum(@delta * @gamma)
		if (@denom1 != 0){
			@h <-+ (@delta * @delta')/@denom1
		}
	}
	
	if (@iter >= @minit) {
		if (@relcoef > 0) {
			@xtmp <- vector(max(t(hconcat(rep(.5,@n),abs(@x)))))
			if (max(abs(@delta) - @relcoef*@xtmp) < 0){
				# relative accuracy for x accomplished
				@critno <- 1
				break
			}
		}
		if (@relcrit > 0) {
			if (abs(@f - @oldf) <= @relcrit*max(.5,abs(@f))){
				# relative accuracy for f accomplished
				@critno <- 2
				break;
			}
		}
		if (alltrue(@dg > 0,sqrt(sum(@g^2)) < @dg)) {
			# size of ||gradient|| accomplished
			@critno <- 3
			break;
		}
	}
}
if (alltrue(@record > 0, @nrec > 0)) {
	@recname <- if (@imeth == 1) {
		"BFSRECORD"
	} elseif(@imeth == 2) {
		"DFPRECORD"
	} else {
		"BROYDNRECORD"
	}
	
	<<@recname>> <- structure(xvals:padto(delete(@xvals,return:T),@nrec),\
					    funvals:padto(delete(@funvals,return:T),@nrec),\
						gradients:padto(delete(@gradients,return:T),@nrec))
	delete(@nrec,@recname)
}

if (@gold > 0) {
	delete(@x1,@x2,@left,@phi1,@phi2,@j,@l,@r)
}
delete(@n,@gold,@d,@gamma,@delta,@hgam,@delhgam,@denom1,\
	   @maxit,@minit,@print,@record,\
	   @lambda,@oldf,@oldg,@delhgam,silent:T)
structure(x:vector(delete(@x,return:T)), f:delete(@f,return:T),\
		gradient:vector(delete(@g,return:T)),h:delete(@h,return:T),\
		iterations:delete(@iter, return:T), status:delete(@critno,return:T))
%minimizer%

bfs  macro dollars
) Macro to use Broyden-Fletcher-Shanno method to minimize a function
) Based on Dahlquist and Bjorck, Numerical methods, Prentice Hall, 1974
) p. 443
) usage:
)   bfs(x0, fun [,params] [,h:invhes] [,goldsteps:m]\
)            [,criteria:vector(nsigx,nsigfun,dg)]\
)            [, maxit:maxit] [,minit:minit]\
)            [,printwhen:d1] [,recordwhen:d2]))
)    x0          REAL length k vector of starting values
)    fun         macro to evaluate function (f <- fun(x,0 [,params]))
)                and gradient (g <- fun(x,1 [,params]))
)    params      optional variable with information for fun
)    invhes      symmetric k by k REAL matrix, a starting approximation
)                to the inverse of the Hessian matrix [identity]
)    m           positive integer = number of cycles using golden mean
)                line search, default [5]
)    maxit > 0   integer = maximum number of iterations [30]
)    minit >= 0  integer = minimum number of iterations [0]
)    nsigx       integer; desired significant digits in coefficients [5]
)    nsigfun     integer; desired significant digits in minimum value [8]
)    dg          REAL convergence criterion for ||gradient|| [-1]
)    d1 >= 0     integer [0].  When d1 > 0, partial results are printed
)                at iterations d1, 2*d1, ...
)    d2 >= 0     integer [0].  When d2 > 0, at iterations d2, 2*d2, ...,
)                partial results are recorded in side effect variable
)                BFSRECORD = structure(xvals, funvals, gradient) where
)                xvals[i,], funvals[i], gradient[i,] are values of x,
)                f(x) and gradient(x) at iteration i*d2.
) Return value:
)   structure(x:minimum, f:minimizedValue,gradient:gradient,\
)      h:invhessian, iterations:number_of_interations,status:critno)
) C. Bingham (kb@stat.umn.edu)
)) 011112 Component for interations is now 'iterations' to match help
)) 011116 added keywords minit, printwhen, recordwhen
)) 011118 cleaned up golden section search; removed keyword 'epsilon';
))        added keyword 'criteria'
)) 011210 Changed to use macro minimizer()
)) Version 011210
#$S(start, fun [, params] [,h:invhes] [, goldsteps:m] \
#   [,criteria:vector(nsigx,nsigfun,dg)] [, maxit:maxit] [, minit:minit] \
#   [,printwhen:d1] [,recordwhen:d2])
if ($v < 2 || $v > 3) {
	error("usage:$S(start, fun [, params] [,h:invhes] \\
    [, goldsteps:m] [,criteria:vector(nsigx,nsigfun,dg)]  \\
    [, maxit:maxit] [, minit:minit] [,printwhen:d1] [,recordwhen:d2])",\
	macroname:F)
}
if (!ismacro(minimizer)) {
	getmacros(minimizer,silent:T)
}
minimizer($0,method:"$S")
%bfs%
		
dfp  macro dollars
) Macro to use Davidon-Fletcher-Powell method to minimize a function
) Based on Dahlquist and Bjorck, Numerical methods, Prentice Hall, 1974
) p. 442
) usage:
)   dfp(x0, fun [,params] [,h:invhes] [,goldsteps:m]\
)            [,criteria:vector(nsigx,nsigfun,dg)]\
)            [, maxit:maxit] [,minit:minit]\
)            [,printwhen:d1] [,recordwhen:d2]))
)    x0          REAL length k vector of starting values
)    fun         macro to evaluate function (f <- fun(x,0 [,params]))
)                and gradient (g <- fun(x,1 [,params]))
)    params      optional variable with information for fun
)    invhes      symmetric k by k REAL matrix, a starting approximation
)                to the inverse of the Hessian matrix [identity]
)    m           positive integer = number of cycles using golden mean
)                line search, default [5]
)    maxit > 0   integer = maximum number of iterations [30]
)    minit >= 0  integer = minimum number of iterations [0]
)    nsigx       integer; desired significant digits in coefficients [5]
)    nsigfun     integer; desired significant digits in minimum value [8]
)    dg          REAL convergence criterion for ||gradient|| [-1]
)    d1 >= 0     integer [0].  When d1 > 0, partial results are printed
)                at iterations d1, 2*d1, ...
)    d2 >= 0     integer [0].  When d2 > 0, at iterations d2, 2*d2, ...,
)                partial results are recorded in side effect variable
)                DFPRECORD = structure(xvals, funvals, gradient) where
)                xvals[i,], funvals[i], gradient[i,] are values of x,
)                F(x) and gradient(x) at iteration i*d2.
) Return value:
)   structure(x:minimum, f:minimizedValue,gradient:gradient,\
)      h:invhessian, iterations:number_of_interations,status:critno)
)) C. Bingham (kb@stat.umn.edu)
)) 011113 Component for interations is now 'iterations' to match help
)) 011116 added keywords minit, printwhen, recordwhen
)) 011118 cleaned up golden section search; removed keyword 'epsilon';
))        added keyword 'criteria'
)) 011210 Changed to use macro minimizer()
)) Version 011210
#$S(start, fun [, params] [,h:invhes] [, goldsteps:m] \
#   [,criteria:vector(nsigx,nsigfun,dg)] [, maxit:maxit] [, minit:minit] \
#   [,printwhen:d1] [,recordwhen:d2])
if ($v < 2 || $v > 3) {
	error("usage:$S(start, fun [, params] [,h:invhes] \\
    [, goldsteps:m] [,criteria:vector(nsigx,nsigfun,dg)]  \\
    [, maxit:maxit] [, minit:minit] [,printwhen:d1] [,recordwhen:d2])",\
	macroname:F)
}	
if (!ismacro(minimizer)) {
	getmacros(minimizer,silent:T)
}
minimizer($0,method:"$S")
%dfp%

broyden  macro dollars
) Macro to use Broyden's method to minimize a function.  No line
) search is used.  Some safeguards are said to be needed but are not
) implemented
) Based on Dahlquist and Bjorck, Numerical methods, Prentice Hall, 1974
) p. 443
) usage:
)   broyden(x0, fun [,params] [,h:invhes]\
)            [,criteria:vector(nsigx,nsigfun,dg)]\
)            [, maxit:maxit] [,minit:minit]\
)            [,printwhen:d1] [,recordwhen:d2]))
)    x0          REAL length k vector of starting values
)    fun         macro to evaluate function (f <- fun(x,0 [,params]))
)                and gradient (g <- fun(x,1 [,params]))
)    params      optional variable with information for fun
)    invhes      symmetric k by k REAL matrix, a starting approximation
)                to the inverse of the Hessian matrix [identity]
)    maxit > 0   integer = maximum number of iterations [30]
)    minit >= 0  integer = minimum number of iterations [0]
)    nsigx       integer; desired significant digits in coefficients [5]
)    nsigfun     integer; desired significant digits in minimum value [8]
)    dg          REAL convergence criterion for ||gradient|| [-1]
)    d1 >= 0     integer [0].  When d1 > 0, partial results are printed
)                at iterations d1, 2*d1, ...
)    d2 >= 0     integer [0].  When d2 > 0, at iterations d2, 2*d2, ...,
)                partial results are recorded in side effect variable
)                BROYDNRECORD = structure(xvals, funvals, gradient) where
)                xvals[i,], funvals[i], gradient[i,] are values of x,
)                F(x) and gradient(x) at iteration i*d2.
) Return value:
)   structure(x:minimum, f:minimizedValue,gradient:gradient,\
)      h:invhessian, iterations:number_of_interations, status:critno)
)) C. Bingham (kb@stat.umn.edu)
)) 011113 Component for interations is now 'iterations' to match help
)) 011116 added keywords minit, printwhen, recordwhen
)) 011118 removed keyword 'epsilon'; added keyword 'criteria'
)) 011210 Changed to use macro minimizer()
)) Version 011210
#$S(start, fun [, params] [,h:invhes] \
#   [,criteria:vector(nsigx,nsigfun,dg)] [, maxit:maxit] [, minit:minit] \
#   [,printwhen:d1] [,recordwhen:d2])
if ($v < 2 || $v > 3) {
	error("usage:$S(start, fun [, params] [,h:invhes] \\
    [,criteria:vector(nsigx,nsigfun,dg)] [, maxit:maxit] [, minit:minit] \\
    [,printwhen:d1] [,recordwhen:d2])",\
	macroname:F)
}	
if (!ismacro(minimizer)) {
	getmacros(minimizer,silent:T)
}
minimizer($0,method:"$S")
%broyden%
	
neldermead   MACRO  DOLLARS
) A macro implementing function minimization using the simplex method.
) The minimum found will often be a local,  not a global,  minimum.
)
) For details,  see Nelder & Mead,  The Computer Journal,  January 1965
)
) Based on a Fortran program with the following history
)   Programmed by D.E.Shaw, CSIRO,  Division of Mathematics & Statistics
)      P.O. Box 218,  Lindfield,  N.S.W. 2070
)    With amendments by R.W.M.Wedderburn,Rothamsted Experimental Station,
)      Harpenden,  Hertfordshire,  England
)    Further amended by Alan Miller, CSIRO, Division of Mathematics &
)      Statistics, Private Bag 10,  Clayton,  Vic. 3168
)
) MacAnova version by C. Bingham, kb@stat.umn.edu, June 2000 & Dec 2001
) based on the revision of 11 August 1991
)
) Usage:
)   neldermead(fun, start, steps [,data] [, maxeval:m] [, print:ip]\
)      [,stopcrit:crit] [,checkwhen:n] [,quad:F or T] [,simpcrit:s] \
)      [,tuning:vector(alpha,beta,gamma)])
)
)   fun     macro.  fun(x, data) should return the value F(x) of the
)           function to be minimized, where x is a REAL vector.  fun()
)           should ignore second argument if not used.
)   start   REAL vector of starting values, including elements of x
)           that remain fixed as controlled by steps
)   steps   REAL vector of initial step sizes for building the simplex.
)           step[i] = 0 means x[i] remains fixed.
)   data    optional variable of data or fixed parameters; may be structure
)   m       Integer > 0, the maximum no. of function evaluations allowed;
)           default = 1000
)   crit    Scalar > 0, stopping criterion on standard deviation of values
)           on a simplex; default = 1e-5
)   n       the stopping rule is applied after every n function
)           evaluations; default = 10
)   quad:T  Fit a quadratic surface for final polishing of minimum
)   quad:F  Don't fit a quadratic surface for final polishing of minimum
)           quad:F is default
)   s       Small scalar > 0, criterion for expanding the simplex to
)           overcome rounding errors before fitting the quadratic surface;
)           default = 1e-8.
)   ip      Integer controlling printing: ip < 0 suppresses printing;
)           ip = 0 prints values of x and fun(x) after initial evidence
)           of convergence; ip > 0 is like ip = 0 but also prints progress
)           reports after every ip evaluations, and prints the initial
)           simplex. Default = -1
)   alpha   Reflection coefficient, a positive scalar, default = 1
)   beta    Contraction coefficient, a positive scalar, default = .5
)   gamma   Expansion coefficient, a positive scalar, default = 3
)           These are Nelder-Mead suggested defaults
)
) The value returned is
)   structure(x:xmin, f:Fmin, invhessian:V, neval:Neval, status:Status)
)    xmin    REAL vector of location of minimum
)    Fmin    Attained minimum = value of F(xmin)
)    V       REAL matrix = inverse of Hessian matrix.  When
)            F(x) = -log(L(x)), L(x) a likelihood function, V is
)            the inverse of the observed information, and is an estimate
)            of the variance matrix of the MLE minimizing F(x)
)            When F(x) = sum(residuals^2), 2*MSE*V = estimated variance
)            matrix of the parameters, where MSE = F(xmin)/(n-np), n =
)            sample size and np = number of parameters
)    Neval   Number of function evaluations
)    Status  Termination status
)             = 0 for successful termination
)             = 1 if maximum no. of function evaluations exceeded
)             = 2 if information matrix is not + ve semi - definite
)   With ip >= 0, the returned structure is "invisible"; it can be
)   assigned but will not be printed automatically
)
) For advice on usage, see help subtopic neldermead:"advice_on_usage")
)
)) 000825 kb@stat.umn.edu
)) 011203 minor bugs corrected; optional argument data added; fossil
))        syminv removed
)) 011207 Added keyword 'tuning' for experimenting with reflection,
))        contraction and expansion coefficients.
))        Changed name of component 'covariance' to 'invhessian' and
))        removed factor of 0.5 in it's definition.
)) 011213 steps can now be scalar, equivalent to rep(steps,length(start))
)) 011216 cleaned up handling of near singular hessian; modified printed
))        output and made output invisible with print:ipm ip >= 0
))        changed orientation so that columns of @g are points, not rows
))        keyword 'nloop' is now 'checkwhen'
#$S(fun, bstart, steps, maxit:m, print:k, crit:crit, nloop:m,\
#   quad:T, stopcrit:crit)
# returns structure(x:xmin, f:Fmin, invhessian:V, neval:neval,status:s)
@fun <- argvalue($1,"argument 1","macro")
@p <- vector(argvalue($2, "starting values", "real vector nonmissing"))
@steps <- vector(argvalue($3, "step sizes", "nonneg vector"))
@data <- if ($v > 3) {
	argvalue($4, "data")
} else {
	NULL
}
@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(NoTaKeY:NULL)
}
@maxit <- keyvalue(@keys, "maxit*", "count", default:50)
@maxeval <- keyvalue(@keys, "maxeval*", "count", default:1000)
@print <- keyvalue(@keys, "print*", "integer scalar", default:-1)
@stopcr <- keyvalue(@keys, "stopcrit*", "positive number", default:1e-5)
@checkwhen <- keyvalue(@keys, "checkwhen", "positive count")
if (isnull(@checkwhen)) {
	# for backward compatibility
	@checkwhen <- keyvalue(@keys, "nloop", "positive count",default:10)
}
@quad <- keyvalue(@keys, "quad*", "TF", default:F)
@simpcrit <- keyvalue(@keys, "simpcrit*", "positive number", default:1e-8)
@tuning <- keyvalue(@keys, "tun*","positive real vector",\
					default:vector(1,.5,2)) #Nelder & Mead defaults
delete(@keys)

@nop <- length(@p)
if (isscalar(@steps)) {
	@steps <- rep(@steps,@nop)
} elseif (length(@steps) != @nop){
	error("Number of step sizes different from number of parameters")
}
@active <- run(@nop)[@steps != 0]
@nap <- length(@active)
if (@nap == 0){
	@f <- @fun(@p,delete(@data,return:T))
	delete(@fun,@steps,@maxit,@maxeval,@print,@stopcr,@checkwhen,@quad,\
		   @simpcrit,@nop,@nap)
	return(structure(x:vector(delete(@p,return:T)), f:vector(delete(@f,return:T)),\
				   invhessian:NULL, neval:1,status:0))
}
@np1 <- @nap + 1 # number of points in each simplex

 # reflection, contraction and expansion coefficients
@reflect  <- @tuning[1]  # alpha (default 1 suggested by Nelder-Mead)
@contract <- @tuning[2]  # beta (default .5)
@expand   <- @tuning[3]  # gamma (default 2)

@g <- matrix(rep(@p,@nap+1),@nop) # current simplex, @nop x @nap+1
@g[hconcat(@active,run(@nap)+1)] <- @p[@active] + @steps[@active]

@h <- rep(0, @np1) # function values at simplex vertices
for (@i,1,@np1){
	@pstar <- @g[,@i]
	@h[@i] <- @fun(@pstar,@data)
}
@neval <- @np1
if (@print > 0){
	print(@h,name:"Values on initial simplex")
}

@loop <- @iflag <- 0

@keepPstar <- F
while(T){
	# Start of main cycle
	@loop <-+ 1
	@tmp <- grade(@h)
	@imin <- @tmp[1]
	@imax <- @tmp[@np1]

	@hmin <- @h[@imin]
	@hmax <- @h[@imax]
	
	# find the centroid of the vertices other than p[imax]
	@pbar <- describe(@g[,-@imax]',mean:T)

	@pstar <- @pbar + @reflect*(@pbar-@g[,@imax])
	@hstar <- @fun(@pstar,@data)
	@neval <-+ 1
	if (alltrue(@print > 0, @neval %% @print == 0)){
		print(paste(iw:4,"neval =",@neval,", value =", @hstar))
		print(@pstar, name:"location", labels:F)
	}
	# if hstar < hmin,  reflect pbar through pstar,
	#   hstst = function value at pstst.
	
	if (@hstar < @hmin){
		@pstst <- @pbar + @expand * (@pstar - @pbar)
		@hstst <- @fun(@pstst,@data)
		@neval <-+ 1
		if (alltrue(@print > 0, @neval %% @print == 0)){
			print(paste(iw:4,"neval =",@neval,", value =", @hstst))
			print(@pstst, name:"location", labels:F)
		}

		if (@hstst < @hmin){
			@g[@active,@imax] <- @pstst[@active]
			@h[@imax] <- @hstst
		} else {
			@keepPstar <- T
		}
	} else {
		# hstar is not < hmin.
		# Test whether it is < function value at some point other than
		# p[imax]. If it is replace p[imax] by pstar & hmax by hstar.
		
		if (@hstar < max(@h[-@imax])) {
			@keepPstar <- T
		}
		
		# hstar > all function values except possibly hmax.
		# if hstar <= hmax,  replace p[imax] by pstar & hmax by hstar.
		if (!@keepPstar) {
			if (@hstar <= @hmax) {
				@g[@active,@imax] <- @pstar[@active]
				@hmax <- @h[@imax] <- @hstar
			}
			
			# contracted step to the point pstst,
			# hstst = function value at pstst.
			@pstst <- @pbar + @contract*(@g[,@imax] - @pbar)
			@neval <-+ 1

			if (alltrue(@print > 0, @neval %% @print == 0)){
				print(paste(iw:4,"neval =",@neval,", value =", @hstst))
				print(@pstst, name:"location", labels:F)
			}
			@hstst <- @fun(@pstst,@data)
			
			# if hstst < hmax replace p[imax] by pstst & hmax by hstst
			if (@hstst <= @hmax) {
				@g[@active,@imax] <- @pstst[@active]
				@h[@imax] <- @hstst
			} else {
				# hstst > hmax.
				# Shrink the simplex by replacing each point,  other than
				# the current minimum,  by a point mid - way between its
				# current position and the minimum.

				for (@i,1,@np1) {
					if (@i != @imin) {
						@g[@active,@i] <- (@g[@active,@i] + @g[@active,@imin])/2
						@p <- @g[,@i]
						@h[@i] <- @fun(@p,@data)
						@neval <-+ 1
						if (alltrue(@print > 0, @neval %% @print == 0)){
							print(paste(iw:4,"neval =",@neval,", value =", @h[@i]))
							print(@p, name:"location", labels:F)
						}
					}
				}
			}
		}
	}

	# replace maximum point by pstar & h[imax] by hstar

	if (@keepPstar){
		@keepPstar <- F
		@g[@active,@imax] <- @pstar[@active]

		@h[@imax] <- @hstar
	}
	# if loop >= checkwhen test for convergence,  otherwise repeat main cycle
	if (@loop >= @checkwhen) {
		# Calculate mean & standard deviation of function values for the
		# current simplex.
		@loop <- 0
		@hmean <- sum(@h)/@np1
		@hstd <- sqrt(sum((@h - @hmean)^2)/@np1)

		# If the rms > stopcr,  set iflag & loop to 0.0 and go to the
		# start of the main cycle again.
		
		if (@hstd > @stopcr && @neval <= @maxeval) {
			# not yet converged or hit limit
			@iflag <- 0
		} else {
			# find the centroid of the current simplex
			# and the function value there
			@p[@active] <- describe(@g[@active,]',mean:T)
			@func <- @fun(@p,@data)
			@neval <-+ 1
			if (alltrue(@print > 0, @neval %% @print == 0)){
				print(paste(iw:4,"neval =",@neval,", F(x) =", @func))
				print(x:@p, labels:F)
			}
			# Test whether the number of function values allowed,
			# maxeval,  has been overrun;if so, exit with ifault = 1.
			
			if (@neval > @maxeval) {
				break
			}
			# Convergence criterion satisfied.
			# If iflag = 0,  set iflag & save hmean.
			# If iflag = 1 & change in hmean <= stopcr then
			# search is complete.

			if (@print >= 0) {
				print("")
				print("Evidence of convergence")
				print(@p, name:"Centroid of last simplex", labels:F)
				print(paste("Function value at centroid =",@func))
			}
			if (@iflag <= 0) {
				@iflag <- 1
			} elseif (abs(@savemn - @hmean) < @stopcr) {
				break
			}
			@savemn <- @hmean
		}
	}
}

if (@neval > @maxeval) {
	if (@print >= 0) {
		print("")
		print(paste("WARNING: Number of function function evaluations >",\
			@maxeval))
		print(paste("RMS of function values of last simplex =",@hstd))
		print(@p, labels:F, name:"Centroid of last simplex")
		print(paste("Function value at centroid =",@func))
	}
	return(structure(x:vector(delete(@p,return:T)), f:vector(delete(@func,return:T)),\
				   invhessian:NULL,neval:delete(@neval,return:T),status:1))
}

if (@print > 0 || !@quad && @print == 0) {
	print(paste("Minimum found after",@neval,"function evaluations"))
	print(vector(@p),labels:F,name:"Location of minimum")
	print(paste("Function value at minimum =",@func))
}

if (@quad) {
	if (@print >= 0) {
		print("")
		print("Fitting quadratic surface near supposed minimum")
	}
	# If necessary, expand the final simplex to overcome rounding errors.
	
	@neval1 <- 0

	for (@i, 1, @np1) {
		while(abs(@h[@i] - @func) < @simpcrit) {
			@g[@active,@i] <- @g[@active,@i] + (@g[@active,@i] - @p)
			@pstst <- @g[,@i]
			@h[@i] <- @fun(@pstst,@data)
			@neval1 <-+ 1
		}
	}
	# @g has simplex; @h has function values at vertices
	# Function values are calculated and saved at an additional nap points.
	@aval <- rep(0,@nap)
	for (@i, 1, @nap) {
		@pstar <- (@g[,1] + @g[,@i+1])/2
		@aval[@i] <- @fun(@pstar,@data)
	}
	# @aval has nap values half-way down the simplex edges from @g[,1]
	# Calculate and save in bmat the matrix if 2nd derivatives as
	# estimated from quadratic surface fitting the simplex and all the
	# points bisecting the edges.
	@a0 <- @h[1]

	# Set diagonal of bmat and part of off-diagonal
	@bmat <- dmat(2*@h[-1]) + 2*@a0 - 2*(@aval + @aval')
	# Compute remaining part of off-diagonal
	if (@nap > 1) {
		for (@i, 2, @nap) {
			@pstst <- (@g[,run(2,@i)] + @g[,@i + 1])/2
			for (@j, 1, @i-1) {
				@pststj <- @pstst[,@j]
				@hstst <- @fun(@pststj,@data)
				@bmat[@i,@j] <- @bmat[@j,@i] <- @bmat[@j,@i] + 2*@hstst
			}
		}
	}
	@neval1 <-+ @nap*(@nap + 1)/2
	
	# Save the vector of estimated first derivatives in @aval.
	@aval <- 2*@aval - (@h[-1] + 3*@a0)/2

	# @aval is .5*(1st derivative of function computed by @fun()), length(@nap)
	# @bmat is @nap x @nap .5*(matrix of 2nd derivatives)
	# Derivatives are in coordinate system where g[1,] <-> (0,0,...0)
	# and g[j,] <-> ej = (0,...,0,1,0,...0), parallel the jth coordinate axis
	NMHESSIAN <- 2*@bmat
	@R <- cholesky(@bmat, pivot:T, nonposok:T)
	if (isnull(@R)) { # not positive definite
		@eigs <- eigen(@bmat)
		@vals <- @eigs$values
		@vecs <- @eigs$vectors
		@nonpos <- abs(min(@vals)) >= 1e-8*max(@vals)
		if (!@nonpos) { #might be close enough to nonnegative definite
			@J <-  abs(@vals) < 1e-8 * abs(max(@vals))
			if (sum(@J) > 0) {
				# try decomposing @bmat slightly adjusted to make PD
				@bmat1 <- @bmat + 1e-8*@vals[1] * @vecs[,@J] %*% @vecs[,@J]'
				@R <- cholesky(@bmat1, pivot:T, nonposok:T)
				@nonpos <- isnull(@R)
			}
		}
		if (@nonpos) {
			if (@print >= 0) {
				print("WARNING: Matrix of estimated second derivatives not positive definite")
				print("WARNING: Minimum probably not found")
			}
			@neval <- @neval + delete(@neval1,return:T)
			return(structure(x:vector(delete(@p,return:T)), f:vector(delete(@func,return:T)),\
						   invhessian:NULL, neval:delete(@neval,return:T),status:2))
		}
	}
	@pivot <- @R$pivot
	@R <- @R$r
	@invpivot <- grade(@pivot)
	# @bmat[@pivot,@pivot] = @R' %*% @R
	# @bmat = @R[@invpivot,@invpivot]' %*% @R[@invpivot,@invpivot]

	@J2 <- diag(@R)/@R[1,1] > 5e-5
	@bmatinv <- dmat(@nap,0)
	@R <- solve(@R[@J2,@J2])
	@bmatinv[@J2,@J2] <- @R %C% @R
	@bmatinv <- @bmatinv[@invpivot,@invpivot] # rearrange
	@h <- @bmatinv %*% @aval

	# Find the position,  pmin,  & value,  ymin,  of the
	# minimum of the quadratic.
	@ymin <- @a0 - sum(@h * @aval) # @a0 - @aval' %*% solve(@bmat, @aval)

	# The matrix Q of Nelder & Mead is calculated and stored in q.
	@q <- @g[,-1] - @g[,1]
	@pmin <- @g[,1] - @q %*% @h
	delete(@g)
	if (@print >= 0) {
		print(paste("Minimum of quadratic surface =",@ymin))
		print(vector(@pmin), name:"Location of minimum")
		print("If this differs by much from the minimum found by the simplex search")
		print("the minimum may be false or the inverse Hessian matrix may be inaccurate")
	}
	# Calculate function value at the minimum of the quadratic
	@hstar <- @fun(@pmin, @data)
	@neval1 <-+ 1
	
	# if hstar < func,  replace search minimum with quadratic minimum
	if (@hstar < @func) {
		@func <- @hstar
		@p <- @pmin
		if (@print >= 0) {
			print(paste("Function value at minimum of quadratic =",@func))
		}
	}

	# q*bmatinv*q'/2 is calculated & is stored in vc

	@vc <- @q %*% (@bmatinv %C% @q)/2

	# @vc is inverse of matrix of 2nd derivatives of F(x) computed by @fun()
	# If F(x) = -log(L(x)), @vc = inverse of the observed information matrix,
	# that is, @vc = estimated variance matrix of the parameters
	# If F(x) = sum(residuals^2), 2*MSE*@vc = estimated variance matrix
	# of the parameters, where MSE = F(xmin)/(n-np), n = sample size
	# and np = number of parameters
	
	if (@print >= 0) {
		print("")
		print(@vc ,name:"V = Generalized inverse of the Hessian matrix")
		print("If the function minimized is -log(likelihood),")
		print("  V = the variance/covariance matrix of the parameters.");
		print("If the function was a sum of squared residuals,");
		print("  2*MSE*V = the variance/covariance matrix of the parameters");

		# Compute and print Hessian matrix in actual coordiates
		@q <- solve(@q[@active,])
		NMHESSIAN <- @tmp <- 2 * @q %c% @bmat %*% @q
		@imat <- dmat(@nop, 0)
		@imat[@active,@active] <- delete(@tmp,return:T)
		print("")
		print(delete(@imat,return:T), name:"Hessian matrix (information matrix)")

		# Compute and print standard deviations of parameters
		@sd <- rep(0,@nop)
		@sd[@active] <- @d <- sqrt(diag(@vc)[@active])
		print("")
		print(delete(@sd,return:T),name:"Standard deviations of parameters")

		# Compute and print correlation matrix of parameters
		@cor <- dmat(@nop,0)
		@cor[@active,@active] <- @vc[@active,@active]/(@d*@d')
		print("")
		print(@cor, name:"Correlation matrix of parameters")
		delete(@d,@cor)
	}
	if (@print > 0) {
		print(paste(@neval1,"more function evaluations were used"))
	}
	delete(@q)
	@neval <-+ delete(@neval1,return:T)
} else {
	@vc <- NULL
}
delete(@hstar,@h,@g, @keepPstar, @hstd, @hmean, @stopcr, @loop,\
	@checkwhen, @savemn, silent:T)

@result <- structure(x:vector(delete(@p,return:T)), f:vector(delete(@func,return:T)),\
	neval:delete(@neval,return:T), invhessian:delete(@vc,return:T),\
	status:0)
delete(@result,return:T,invis:delete(@print,return:T) >= 0)
%neldermead%

===> levmar <===
levmar        MACRO DOLLARS
) Macro for fitting by minimizing sums of squares of residuals.
) Usage:
) levmar(b, x, y [, f], param [,resid:residmac] [,deriv:deriv]\
)        [,active:active] [,crit:crvec, maxit:itmax, minit:itmin, print:T])
) b        REAL vector of starting values for coefficients
) x        REAL variable, usually a vector or matrix with nrows(x) = nrows(y)
) y        REAL vector of data to be fit
) f        macro called as fit <- f(b,x,param); required without
)          resid:residmac
) param    NULL or a vector or structure of additional parameters for f
) residmac A macro called as residmac(b, x, y, param) to compute a vector
)          of residuals of length nrows(y).  When f is an argument,
)          'resid' should not be used and residmac(b,x,y,param is
)          essentially y - f(b,x,param)
) crvect   vector(numsig, nsigsq, delta), 3 criteria for convergence
) deriv    optional macro; deriv(b,x,y,param,j) computes derivative of
)          f(b,x,param) or of -resmac(b,x,y,param) with respect to b[j],
)          returning a vector of length nrows(y)
) active   LOGICAL vector the same length as b; active[i] = F means
)          b[i] remains constant
) itmax    maximum number of iterations permitted (default 30)
) itmin    minimum number of iterations performed (default 1)
) print    When T, partial results are printed on each iteration
)
) levmar() returns structure(coefs,hessian,jacobian,gradient,rss,residuals,
)        nobs, iter, iconv)
)
)) Other macros used (should be loaded automatically if available)
))  _cgrad()
))  _lmout() (only with print:T)
)) Written by C. Bingham, December 1998 based on a Fortran program of
)) Ken Brown.  See
))  Brown,K.,M. and Dennis,J.,E., Derivative free analogues
))  of the Levenberg-Marquardt and Gauss algorithms for
))  nonlinear least squares approximation.  Numerische
))  Mathematik, Vol. 18, pp. 289-297 (1972)
))
))  Brown,K.,M., Computer oriented methods for fitting tabular data
))  in the linear and nonlinear least squares sense,
))  Technical Report No. 72-13, University of Minnesota
))  Department of Computer and Information Sciences
))
)) 010724 modified so that hessian returned is always computed from
))        the jacobian returned
) Version 010724
# $S(b, x, y [, f], param [,resid:residmac] [,deriv:deriv]\
#        [,active:active] [,crit:crvec, maxit:itmax, minit:itmin, print:T])
# b        REAL vector of starting values for coefficients
# x        REAL variable, usually a vector or matrix with nrows(x) = nrows(y)
# y        REAL vector of data to be fit
# f        macro called as fit <- f(b,x,param); required without
#          resid:residmac
# param    NULL or a vector or structure of additional parameters for f
# residmac A macro called as residmac(b, x, y, param) to compute a vector
#          of residuals of length nrows(y).  When f is an argument,
#          'resid' should not be used and residmac(b,x,y,param) is
#          essentially y - f(b,x,param)
# crvect   vector(numsig, nsigsq, delta), 3 criteria for convergence
# deriv    optional macro; deriv(b,x,y,param,j) computes derivative of
#          f(b,x,param) or of -resmac(b,x,y,param) with respect to b[j],
#          returning a vector of length nrows(y)
# active   LOGICAL vector the same length as b; active[i] = F means
#          b[i] remains constant
# itmax    maximum number of iterations permitted (default 30)
# itmin    minimum number of iterations performed (default 1)
# print    When T, partial results are printed on each iteration
#
# levmar() returns structure(coefs,hessian,jacobian,gradient,rss,residuals,
#        nobs, iter, iconv)
@b <- argvalue($1,"argument 1","nonmissing real vector")
@x <- argvalue($2,"argument 2","nonmissing real matrix")
@y <- argvalue($3,"argument 3","nonmissing real vector")

@npar <- length(@b)

@crit <- keyvalue($K,"crit*","nonmissing real vector",\
		default:vector(5, 8, -1))
if (max(@crit) <= 0){
	error("You must have at least one positive criterion value")
}

if (isscalar(@crit)){
	@crit <- vector(@crit,8, -1)
} elseif (length(@crit) == 2){
	@crit <- vector(@crit, -1)
}

@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(notakey:NULL)
}
_resid <- keyvalue(@keys, "resid", "macro")
@deriv <- keyvalue(@keys, "deriv*", "macro")

@active <- keyvalue(@keys,"active","nonmissing logic vector",\
	default:rep(T,@npar))
if (length(@active) != @npar){
	error("length of value 'active' must be length of b")
}

@maxit <- keyvalue(@keys,"maxit*", "count", default:30)
@minit <- keyvalue(@keys,"minit*", "count", default:0)
if(@maxit < @minit){
	error("value for maxit < value for minit")
}

@print <- keyvalue(@keys,"print", "TF", default:F)
delete(@keys)

if (isnull(_resid)){
	@_FUNC <- argvalue($4,"argument 4","macro")
	_resid <- macro(paste("\$3 - ",nameof(@_FUNC),"(\$1,\$2,\$04)",sep:""))
	@param <- if($v > 4){
		$5
	}else{
		NULL
	}
} else {
	@_FUNC <- NULL
	@param <- if($v > 3){
		$4
	}else{
		NULL
	}
}

@relcon <- @crit[1]
@relssq <- @crit[2]
if (@relcon != floor(@relcon) || @relssq != floor(@relssq)){
	error("crit[1] and crit[1] must be integers")
}

@relcon <- if(@relcon > 0){
	10^-@relcon
}else{
	0
}

@relssq <- if(@relssq > 0){
	10^-@relssq
}else{
	0
}

@delta <- @crit[3]

@nobs <- nrows(@y) # not necessarily = nrows(@x)

@iconv <- -5
@ibad <- -99
@dtst <- sqrt(5e-11)
@deps <- sqrt(1e-8)
@isw <- 1

@r <- _resid(@b, @x, @y, @param)

@ssq <- sum(@r^2)

if (!ismacro(_cgrad)){
	getmacros(_cgrad,silent:T)
}
@stuff <- _cgrad(@isw, @b, @x, @y, @param, @r, active:@active, deriv:@deriv)

@gradient <- @stuff$gradient
@jacobian <- @stuff$jacobian

delete(@stuff)
@erl2 <- sqrt(sum(@gradient^2))
if (@print){
	if (!ismacro(_lmout)){
		getmacros(_lmout, silent:T)
	}
	_lmout(0, @b, @ssq, @gradient, @erl2, @iconv)
}

if (@maxit > 0){
	for(@iter,1,@maxit){
		@oldssq <- @ssq

		if (@erl2 <= @delta && @iter >= @minit){
			@iconv <- 3
			break
		}

		@hes <- @jacobian %c% @jacobian
		@dnorm <- @erl2/sqrt(sum(diag(@hes)^2))
		while(T){# as yet ibad is not ever set other than to 0
			@d <- @dnorm * sqrt(diag(@hes))
			@hes1 <- @hes + dmat(@d)
			if (@ibad > 0){
				@hes1 <-+ dmat(.5*diag(@hes1) + @deps)
			}
			@temp <- vector(solve(delete(@hes1,return:T), @gradient))
			break
		}
		@b[@active] <- @b[@active] - @temp
		if (@iter >= @minit && @relcon > 0){
			@btmp <- vector(max(t(hconcat(rep(.5,@npar),abs(@b)))))
			if (max(abs(@temp) - @relcon*@btmp[@active]) < 0){
				@iconv <- 1
				break
			}
		}
		for (@ldiv, 1, 11){
			@r <- _resid(@b,@x,@y,@param)
			@ssq <- sum(@r^2)
			if (@iter >= @minit && @ldiv == 1 && @relssq > 0){
				if (abs(@ssq - @oldssq) <= @relssq*max(.5,@ssq)){
					@iconv <- 2
					break 2;
				}
			}
			if (@ssq < @oldssq){
				break
			}
			@temp <-/ 2
			@b[@active] <- @b[@active] + @temp
		}
		@stuff <- _cgrad(@isw, @b, @x, @y, @param, @r,\
			active:@active, deriv:@deriv)

		@gradient <- @stuff$gradient
		@jacobian <- @stuff$jacobian
		@erl2 <- sqrt(sum(@gradient^2))
		delete(@stuff)

		if (@print && @iter < @maxit){
			_lmout(@iter,@b,@ssq,@gradient,@erl2,@iconv)
		}
		if (@ssq >= @oldssq){
			@iconv <- 4
			break
		}
	}# for(@iter,1,@maxit)
}else{# if (@maxit > 0)
	@iter <- 0
}
@hes <- @jacobian %c% @jacobian # ensure hessian matches jacobian

@iconv <- max(@iconv, 0)
if (@iter > 0 && @print){
	_lmout(@iter,@b,@ssq,@gradient,@erl2,@iconv)
	if (@iconv == 0){
		print(paste("No convergence in",@maxit,"iterations"))
	}elseif (@iconv < 4){
		print(paste("Convergence by criterion", @iconv))
	}else{
		print(paste("WARNING: Can't reduce RSS after halving step 10 times on iteration",@iter))
	}
}# if (@iter > 0 && @print)

if (!isnull(@_FUNC)){
	delete(_resid)
}# clean up
delete(@x, @y, @_FUNC, @crit, @maxit, @minit, @print, @relcon, \
	@relssq, @delta, @npar, @ibad, @dtst, @deps, @oldssq, \
	@erl2, @dnorm, @d, @btmp, silent:T)

structure(coefs:delete(@b, return:T),\
	hessian:delete(@hes, return:T),\
	jacobian:-delete(@jacobian, return:T),\
	gradient:delete(@gradient, return:T),\
	rss:delete(@ssq, return:T),\
	residuals:delete(@r, return:T),\
	nobs:delete(@nobs, return:T),\
	iter:delete(@iter, return:T),\
	iconv:delete(@iconv, return:T))
%levmar%

===> _cgrad <===
_cgrad        MACRO DOLLARS
) Macro used by levmar() to compute the Jacobian and gradient =
) jacobian %c% residuals.  It computes residuals using macro
) _resid() which must be set by levmar() or the macro calling levmar()
) It computes derivatives numerically
) Usage:
)   stuff <- _cgrad(isw, b, x, y, params [,resids],\
)            [, active: active] [, deriv:deriv])
) isw     1 use (f(b+h)-f(b))/h with small h
)         2 use (f(b_h)-f(b-h)/(2*h) with larger h
) b       REAL vector of coefficients
) x       REAL vector or matrix used as argument to _resid
) y       REAL vector being fitted
) resids  optional REAL vector of residuals computed with b
) params  vector or structure of parameters used as argument to _resid or deriv
) active  LOGICAL vector the same length as b
) deriv   macro for computing derivatives called by
)         deriv(b, x, y, param, j), where 1 <= j <= length(b) which
)         needs to work for any j with active[j] true
)
) _cgrad returns structure(gradient, jacobian)
)
) Unless deriv:deriv is an argument, _resid is called as _resid(b,x,y,params)
)
) It does very little argument checking, assuming that the calling
) macro is doing it right.
)
)) Other macros used
))  _resid (not checked for, must be set by calling macro)
)) Written by C. Bingham, December 1998, based on a Fortran subroutine of
)) Ken Brown.  See macro levmar() for references.
) Version 000909
#  stuff <- _cgrad(isw, b, x, y, params [,resids],\
#           [, active: active] [, deriv:deriv])
# returns structure(gradient, jacobian)
@isw <- $1
@b <- $2
@x <- $3
@y <- $4
@npar <- length(@b)
@active <- keyvalue($K, "active", default:rep(T,@npar))
@deriv <- keyvalue($K,"deriv*")
@dud <- isnull(@deriv)

@nactive <- sum(@active)
@nobs <- nrows(@y)
@factor <- if(@isw == 2){
	1e3
}else{
	1
}

if (@isw == 1){
	@rdec <- if ($v < 6){
		_resid(@b,@x,@y,$05)
	}else{ # residuals supplied as arg 6
		$6
	}
}
@k <- 0
for(@j,1,@npar){
	if (@active[@j]){
		@k <-+ 1
		if (@dud){ # compute derivatives numerically
			@bhold <- @b[@j]
			@hh <- max(1e-8 * @factor * abs(@bhold), 5e-11)
			@b[@j] <- @b[@j] + @hh
			@rinc <- _resid(@b, @x, @y, $05)
			if (@isw == 2){
				@b[@j] <- @bhold - @hh
				@rdec <- _resid(@b, @x, @y, $05)
				@hh <- @hh + @hh
			}
			@jac <- (@rinc - @rdec)/@hh
			@b[@j] <- @bhold
		}else{
			@jac <- -@deriv(@b, @x, @y, $05, @j)
		}
		if (!isdefined(@jacobian)){
			@nobs <- nrows(@jac)
			@jacobian <- matrix(rep(0,@nactive*@nobs),@nobs)
		}
		@jacobian[, @k] <- @jac
	}# if (@active[@j])
}# for(@j,1,@npar)

if (@isw == 2){
	@rdec <- if($v < 6){
		_resid(@b,@x,@y,$05)
	}else{
		$6
	}
}
delete(@isw,@bhold,@b,@x,@y,@npar,@nobs,@rinc,@factor,\
	@j,@k,@active,@jac,@deriv,silent:T)

structure(gradient:vector(@jacobian %c% delete(@rdec,return:T)),\
	jacobian:delete(@jacobian,return:T))
%_cgrad%

===> _lmout <===
_lmout        MACRO DOLLARS
) Macro used by levmar() to print out information on each iteration
) _lmout(iter,b, ssq, gradient, erl2, iconv))
) iter     iteration number
) b        current value of parameters
) ssq      current RSS
) gradient current gradient vector
) erl2     norm of gradient (sqrt(sum(gradient^2)))
) iconv    convergence indicator (< 0 means not converged)
)
)) Other macros used
))  None
)) Written by C. Bingham, December 1998
) Version 000909
#$S(iter,theta, ssq, gradient, erl2, iconv)
print(paste(sep:"","Iteration ", $1,":",sep:" ","Coefs =",$2))
print(paste("  RSS = ",$3,", ||Gradient|| = ", $5,sep:""))
if (($6) < 0){
	print(paste("  Gradient =",$4))
}
%_lmout%

binom               MACRO DOLLARS
) Macro to compute binomial coefficients
) usage:
)  bc <- binom(n,k)
)   n and k   non-negative REAL variables.  Elements of n and k
)             need not be integers.  When n and k are both not
)             scalars, they must be vectors or matrices of the same
)             shape.  Corresponding elements of n and k must
)             satisfy n[i] >= k[i]
)   bc        REAL variable the same size and shape as n and k if
)             they have the same dimensions, or the same size and
)             shape as the non-scalar argument otherwise
)             bc[i] = gamma(n[i]+1)/(gamma(k[i]+1)*gamma(n[i]-k[i]+1))
))000220 stripped $$ and fixed bug when n or k non integer
#usage: $S(n,k) computes binomial coefficients
if($v != 2){
	error("usage: $S(n, k), n >= k >= 0", macroname:F)
}
@n <- argvalue($1, "n", "nonnegative")
@k <- argvalue($2, "k", "nonnegative")
if (!isscalar(@n) && !isscalar(@k) &&\
	anytrue(ndims(@n) != ndims(@k), sum(dim(@n) != dim(@k)) != 0)){
	error("n and k different shapes and neither are scalars")
}

if(min(vector(@n-@k)) < 0){
	error("value for k > value for n")
}
@ans <- exp(lgamma(@n+1) - lgamma(@k+1) - lgamma(@n-@k+1))
if (sum(vector(@n,@k) != floor(vector(@n,@k))) == 0) {
	@ans <- round(@ans) # ensure the result is integral
}
delete(@n,@k)
delete(@ans, return:T)
%binom%

factorial       MACRO DOLLARS
) Macro to compute factorial x!, where x > 0.
) Usage:
)  y <- factorial(x)
)   x         REAL variable
)   y         REAL variable the same size and shape as x
)             y[i] = MISSING when x[i] <= -1 or is MISSING
)             y[i] = gamma(x[i]+1) otherwise
)  Computation is as exp(lgamma(x+1)) except that a check
)  is made so that when x[i] is an integer <= 20 the result is exact
)) 000830 C. Bingham, kb2stat.umn.edu
# $S(x), x non-negative REAL
@x <- argvalue($1,"argument","real")
@J <- ismissing(vector(@x))
if (sum(@J) > 0) {
	print("WARNING: $S of MISSING set to MISSING")
}
if (sum(!@J) > 0) {
	@x1 <- vector(@x)[!@J]
	if (min(@x1) <= -1) {
		print("WARNING: element of x <= -1; result is MISSING",\
			macroname:T)
	}
	if (max(@x1) >= 171.6243769563026973) {
		print("WARNING: element of x too large; result is MISSING",\
			macroname:T)
	}
	@warnopt <- getoptions(warnings:T)
	setoptions(warnings:F)
	@J1 <- @x1 == floor(@x1)
	@J2 <- @J1 && @x1 > 14.5
	@J1 <- @J1 && @x1 <= 14.5
	@x1 <- exp(lgamma(@x1+1))
	if (sum(@J1) > 0){
		@x1[@J1] <- round(@x1[@J1])
	}
	if (sum(@J2) > 0){
		@x1[@J2] <- round(@x1[@J2], -3)
	}
	delete(@J1, @J2)
	setoptions(warnings:delete(@warnopt, return:T))
	@x[!@J]  <- @x1
	delete(@x1)
}
delete(@J)
delete(@x,return:T)
%factorial%

continfrac      MACRO DOLLARS
) Macro to compute continued fraction b0+a1/(b1+a2/(b2+a3/(b3+a4/... )))
) usage:
)  result <- contfrac(a,b)
)   a and b         REAL vectors or matrices with nrows(b) = nrows(a)
)                   or nrows(b) = nrows(a) + 1 and ncols(a) = ncols(b)
)                   unless ncols(a) = 1 or ncols(b) = 1
)                   a[1,] contains a1, a[2,] contains a2, etc.
)                   When nrows(b) = nrows(a), b0 is assumed 0,
)                   b[1,] contains b1, b[2,] contains b2, etc.
)                   When nrows(b) = nrows(a) + 1, b[1,] contains b0,
)                   b[2,] contains b1, etc.
)   result          REAL vector of length max(ncols(a),ncols(b))
) Written by C. Bingham (kb@stat.umn.edu)
) Version: 990227
# $S(a,b), a, b vectors or matrices, nrows(a) = nrows(b) or nrows(b)+1
@a <- matrix(argvalue($1,"numerator coefficients","real matrix nonmissing"))
@b <- matrix(argvalue($2,"denominator coefficients","real matrix nonmissing"))
@n <- nrows(@a)

if (@n == nrows(@b)){
	@b <- vconcat(0*@b[1,],@b)
} elseif (@n != nrows(@b)-1){
	error("nrows(arg1) > nrows(arg2) or nrows(arg1) < nrows(arg2) - 1")
}

if(min(ncols(@a),ncols(@b)) != 1 && ncols(@a) != ncols(@b)){
	error("ncols(arg1)>1 and ncols(arg1)>1 with ncols(arg1) != ncols(arg2)")
}
@A1 <- @B2 <- 1
@B1 <- 0
@A2 <- @b[1,]
for(@i,1,@n){
	@A0 <- @A1
	@A1 <- @A2
	@B0 <- @B1
	@B1 <- @B2
	@A2 <- @b[@i+1,]*@A1 + @a[@i,]*@A0
	@B2 <- @b[@i+1,]*@B1 + @a[@i,]*@B0
}
delete(@i,@n,@a,@b,@A0,@A1,@B0,@B1)
vector(delete(@A2,return:T)/delete(@B2,return:T))
%continfrac%

i0              MACRO DOLLARS
) Compute modified Bessel function of the first kind I0(x) using
) approximations from Abramowitz & Stegun.
) Usage:
)  y <- i0(x)
)    x            REAL scalar, vector, matrix, array or structure
)                 with REAL components
)    y            REAL scalar, vector, matrix, array or structure
)                 the same size and shape as x
))000220 Stripped $$, use argvalue(), minor other changes
# usage: y <- $S(x)   Modified Bessel function of first kind
@x <- argvalue($1, "argument")
if(isstruc(@x)){
	@ans <- @x
	for(@i,run(ncomps(@x))){
		@arg <- @x[@i]
		@tmp <- $S(@x[@i])
		@ans[@i] <- if (!isstruc(@arg)){
			@tmp
		}elseif (ncomps(@arg) > 1){
			@tmp
		}else{
			structure(@tmp, compnames:compnames(@arg))
		}
	}
	delete(@i,@arg,@tmp) # prevents symbol buildup
}else{
	if(!isreal(@x)){error("Argument 1 is not REAL")}
	@ans <- @x # construct an object of same shape as @x
	if(anymissing(@x)){
		@J <- ismissing(vector(@x))
		@x[@J] <- 0 # ensure @x has no MISSING
	}
	@z <- abs(@x)
	@I <- vector(@z < 3.75)
	if(sum(@I) != 0){# computations for values < 3.75
		if(!isdefined(__i0_b_)){ #A&S 9.8.1
			__i0_b_ <- vector(1,3.5156229,3.0899424,1.2067492,.2659732,.0360768,\
				.0045813)
		}
		@ans[@I] <- rational((@z[@I]/3.75)^2,__i0_b_);;
	}
	@I <- !@I
	if(sum(@I) != 0){# computations for values >= 3.75
		if(!isdefined(__i0_c_)){ #A&S 9.8.2
			__i0_c_ <- vector(.39894228,.01328592,.00225319,-.00157565,.00916281,\
				-.02057706,.02635537,-.01647633,.00392377)
		}
		@z <- @z[@I]
		@ans[@I] <- exp(@z)*rational(3.75/@z,__i0_c_)/sqrt(@z)
	}
	if(isdefined(@J)){@ans[@J] <- ?; delete(@J)}
	delete(@I,@z)
}
delete(@x)
delete(@ans, return:T) # return value for i0
%i0%

i1              MACRO DOLLARS
) Compute modified Bessel function of the first kind I1(x) using
) approximations from Abramowitz & Stegun.
) Usage:
)  y <- i1(x)
)    x            REAL scalar, vector, matrix, array or structure
)                 with REAL components
)    y            REAL scalar, vector, matrix, array or structure
)                 the same size and shape as x
))000220 Stripped $$, use argvalue(), minor other changes
# usage: y <- $S(x)   Modified Bessel function of first kind
@x <- argvalue($1, "argument")
if(isstruc(@x)){
	@ans <- @x
	for(@i,run(ncomps(@x))){
		@arg <- @x[@i]
		@tmp <- $S(@x[@i])
		@ans[@i] <- if (!isstruc(@arg)){
			@tmp
		}elseif (ncomps(@arg) > 1){
			@tmp
		}else{
			structure(@tmp, compnames:compnames(@arg))
		}
	}
	delete(@i,@arg,@tmp) # prevents symbol buildup
}else{
	if(!isreal(@x)){error("1st argument to $S not REAL")}
	@ans <- @x # construct an object of same shape as @x
	if(anymissing(@x)){
		@J <- ismissing(vector(@x))
		@x[@J] <- 0 # ensure @x has no MISSING
	}

	@I <- vector(abs(@x) < 3.75)
	if(sum(@I) > 0){# computations for values < 3.75
		if(!isdefined(__i1_b_)){ #A&S 9.8.3
			__i1_b_ <- vector(.5,.87890594,.51498869,.15084934,.02658733,\
				.00301532,.00032411)
		}
		@z <- @x[@I]
		@ans[@I] <- @z*rational((@z/3.75)^2,__i1_b_)
	}
	@I <- !@I
	if(sum(@I) > 0){# computations for values >= 3.75
		if(!isdefined(__i1_c_)){ #A&S 9.8.4
			__i1_c_ <- vector(.39894228,-.03988024,-.00362018,.00163801,\
				-.01031555,	.02282967,-.02895312,.01787654,-.00420059)
		}
		@z <- abs(@x[@I])
		@ans[@I] <- sqrt(@z)*exp(@z)*\
			rational(3.75/@z,__i1_c_)/@x[@I]
	}
	if(isdefined(@J)){@ans[@J] <- ?; delete(@J)} #set MISSING
	delete(@I,@z)
}
delete(@x)
delete(@ans,return:T) # return value for i1
%i1%

orthopoly  MACRO DOLLARS
) Macro to compute orthogonal polynomials
)  Usage: orthopoly(x,n [,polycode [,parameters]])
)    x a REAL vector with no MISSING  values
)    n >= 1 an integer
)    polycode an unquoted letter, default 'p'.
)    parameters a REAL scalar or vector
)
)    Generates P-sub-j(x),0<=j<=n
)      Family                polycode  Parameter
)       Legendre                p       none
)       Jacobi                  j       vector(alpha,beta)
)       Gegenbauer              g       alpha
)       Chebyshev 1             t       none
)       Shifted Chebyshev 1     T       none
)       Chebyshev 2             u       none
)       Shifted Chebyshev 2     U       none
)       Laguerre                l       alpha
)       Hermite                 h       none
)       Discrete                d       vector of weights
)
) All but the discrete follow the normalizations used in Abramowitz
) and Stegun and are computed from the recurrence relations given there.
) The discrete polynomials satisfy sum(weights*Poly^2)/sum(weights) = 1
)
)) Written by C. Bingham 000603, kb@stat.umn.edu
)) 000825
# $S(x,n, polycode, [,parameters])
if ($v < 2) {
	error("usage: $S(x,n [,code [,parameters]]), code one of p,j,g,t,u,l,h or d")
}
@x <- vector(argvalue($1,"x", "real nonmissing vector"))
@degree <- argvalue($2,"max degree", "positive count")

@polycode <- if ($v > 2) {
	"$3"
} else {
	"p" # default is Legendre
}

@params <- if ($v > 3){
	argvalue($4, "parameters", "nonmissing real vector")
}else{NULL}

@typeno <- match(@polycode,\
	vector("p","j","g","t","T","u","U","l","h","H","d"),0)
if (@typeno == 0) {
	error(paste("unknown polynomial code '",@polycode,"'",sep:""))
}
@polyname <- vector("Legendre","Jacobi","Gegenbauer","Chebyshev 1",\
	"Shifted Chebyshev 1","Chebyshev 2","Shifted Chebyshev 2","Laguerre",\
	"Hermite","Hermite_e","Discrete")[delete(@typeno,return:T)]

@N <- length(@x)
@poly <- matrix(rep(1,(@degree + 1)*@N),@N)
delete(@N)
@poly[,2] <- @x
@n <- run(@degree)

@discrete <- @polycode == "d"

if (@polycode == "j") {#Jacobi, Abramowitz & Segun 22.7.1
	if (anytrue(length(@params) != 2, min(@params) <= -1)) {
		error(paste("length(params) != 2 or min(params) <= -1 with polycode \"",\
			@polyname,"\"",sep:""))
	}
	@alpha <- @params[1]
	@beta <- @params[2]
	@tmp <- 2*@n + @alpha + @beta
	@a1 <- 2*(@n + 1)*@tmp*(@tmp - @n + 1)
	@a2 <- (@tmp + 1)*(@alpha - @beta)*(@alpha + @beta)/@a1
	@a3 <- @tmp*(@tmp + 1)*(@tmp + 2)/@a1
	@a4 <- 2*(@n + @alpha)*(@n + @beta)*(@tmp + 2)/delete(@a1,return:T)
	delete(@tmp)
	@poly[,2] <- (@alpha - @beta + (@alpha + @beta + 2)*@x)/2
}elseif(@polycode == "p"){#Legendre, Abramowitz & Segun 22.7.10
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a1 <- @n + 1
	@a2 <- rep(0,@degree)
	@a3 <- (2*@n + 1)/@a1
	@a4 <- @n/delete(@a1,return:T)
} elseif(@polycode == "g"){#Gegenbauer, Abramowitz & Segun 22.7.3
	if (!alltrue(isscalar(@params), @params[1] != 0)) {
		error(paste("params not a non - zero scalar with polycode \"",\
			@polyname,"\"",sep:""))
	}
	@alpha <- @params
	@a1 <- @n + 1
	@a2 <- rep(0,@degree)
	@a3 <- 2*(@n + @alpha)/@a1
	@a4 <- (@n + 2*@alpha - 1)/delete(@a1,return:T)
	@poly[,2] <- 2*@alpha*@x
} elseif(@polycode == "t"){#Chebyshev 1, Abramowitz & Segun 22.7.4
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(0,@degree)
	@a3 <- @a2 + 2 #rep(2,@degree)
	@a4 <- @a2 + 1 #rep(1,@degree)
} elseif(@polycode == "T"){#Shifted Chebyshev 1, Abramowitz & Segun 22.7.8
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(-2,@degree)
	@a3 <- @a2 + 6 #rep(4,@degree)
	@a4 <- @a2 + 3 #rep(1,@degree)
	@poly[,2] <- 2*@x - 1
} elseif(@polycode == "u"){#Chebyshev 2, Abramowitz & Segun 22.7.5
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(0,@degree)
	@a3 <- @a2 + 2 #rep(2,@degree)
	@a4 <- @a2 + 1 #rep(1,@degree)
	@poly[,2] <- 2*@x
} elseif(@polycode == "U"){#Shifted Chebyshev 2, Abramowitz & Segun 22.7.9
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(-2,@degree)
	@a3 <- @a2 + 6 #rep(4,@degree)
	@a4 <- @a2 + 3 #rep(1,@degree)
	@poly[,2] <- 4*@x - 2
} elseif(@polycode == "l"){#Laguerre, Abramowitz & Segun 22.7.12
	if (!alltrue(isscalar(@params), @params[1] > -1)) {
		error(paste("params not a scalar > -1 with polycode \"",\
			@polyname,"\"",sep:""))
	}
	@alpha <- @params
	@a1 <- @n + 1
	@a2 <- (2*@n + @alpha + 1)/@a1
	@a3 <- -1/@a1
	@a4 <- (@n + @alpha)/delete(@a1,return:T)
	@poly[,2] <- -@x + @alpha + 1
} elseif(@polycode == "h"){#Hermite, Abramowitz & Segun 22.7.13
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(0,@degree)
	@a3 <- @a2 + 2
	@a4 <- 2*@n
	@poly[,2] <- 2*@x
} elseif(@polycode == "H"){#HermiteE, Abramowitz & Segun 22.7.14
	if (!isnull(@params)){
		error(paste("parameters not allowed with polycode",@polyname))
	}
	@a2 <- rep(0,@degree)
	@a3 <- @a2 + 1
	@a4 <- @n
	@poly[,2] <- @x
} elseif(@discrete){#discrete
	@w <- if(isnull(@params)) {
		  0*@x + 1
	} else {
		if (length(@params) != length(@x) || min(@params) <= 0) {
			error("wrong number of weights or not all positive")
		}
		@params
	}
	if (@degree >= length(unique(@x))) {
		error("degree >= number of unique values in x")
	}
	@w <-/ sum(@w)
	@x <-- sum(@w*@x) #center x
	# polynomials satisfy sum(weights*poly^2) = sum(weights)
	@poly[,2] <- @x/sqrt(sum(@w*@x^2))
}else{
	error(paste(@polyname,"polynomials not yet implemented"))
}
delete(@polycode,@polyname,@params)
if (@degree > 1){
	if (!@discrete) {
		for (@i,2,@degree){#@i is order
			@poly[,@i + 1] <-\
				(@a2[@i-1] + @a3[@i-1]*@x)*@poly[,@i] - @a4[@i-1]*@poly[,@i - 1]
		}
	} else { # normalized so that sum(weights*poly^2) = sum(weights)
		@tmp2 <- @poly[,2];
		for (@i,2,@degree){#@i is order
			@tmp1 <- @poly[,@i-1];
			@a2 <- -sum(@w*@x*@tmp2^2)
			@a4 <- sum(@w*@x*@tmp2*@tmp1)
			@tmp2 <- (@a2 + @x)*@tmp2 - @a4*@tmp1
			@tmp2 <-/ sqrt(sum(@w*@tmp2^2))
			@poly[,@i + 1] <- @tmp2
		}
		delete(@w,@tmp1,@tmp2)
		@a3 <- NULL
	}
}
delete(@a2,@a3,@a4,@i,@n, silent:T)
delete(@poly,return:T)
%orthopoly%

economize    macro  dollars
) Macro to "economize" a power series expansion on [-1, 1]
) Usage:
)   b <- economize(vector(a0,a1,...,an),m)
) where fn(x) = a0 +a1*x + a2*x^2 + ... + an*x^n is the n-th partial
) sum of a power series providing a good approximation on an
) interval I contained in (-1,1).
) b = vector(b0,b1,...,bm) are the coefficients of the powers of
) x in c0 + c1*T1(x) + c2*T2(x) + ... + cm*Tm(x), where
)    c0 + c1*T1(x) + c2*T2(x) + ... + cn*Tn(x) = fn(x)
) and T1(x), ..., Tn(x) are Chebyshev polynomials.
) C. Bingham, kb@stat.umn.edu
))000607
@a <- argvalue($1,"coefficients","nonmissing real vector")
@m <- argvalue($2,"degree","positive count")
@n <- length(@a)-1
if (@m > @n) {
	error(paste("order > highest power"))
}
if (@m < @n) {
	@c <- vector(0, -1)
	@b <- @a[@n+1]
	for(@i,2,@n+1){
		@b <- movavg(@c,vector(@b, @a[@n+2-@i]))
		@b[run(@i-1)] <- @b[run(@i-1)]/2
	}
	@b <- @b[-run(@n - @m)]
	@a <- padto(@a, @m+1)
	@k <- 1
	for (@i, @m+1,1) {
		@b0 <- @b[@i]
		@b <- autoreg(@c,2*@b)
		@a[@k] <- (@b[@i] - @b0)/2
		@b <- @b[-@i]
		@k <-+ 1
	}
	delete(@b,@c,@i,@k,@b0)
	@a <-* 2
}
delete(@m,@n)
delete(@a,return:T)
%economize%

chebcoefs  macro dollars
) Macro to find the coefficients of the expansion in Chebysev polynomials
) of a polynomial.
) Usage:
)   b <- chebcoefs(a), where a = vector(a0,a1,...,an)
)     are coefficients of the polynomial Pn(x)=a0+x1*x+...+an*x^n
)   b is vector(b0,b1,...,bn), where Pn(x) = b0+b1*T1(x)+ ... + bn*Tn(x)
)   is the unique representation in Chebyshev polynomial Tj(x)
) C. Bingham, kb@stat.umn.edu
)) 000607
# usage: b <- $S(a), REAL vector of polynomial coefficients
@a <- argvalue($1,"arg 1","nonmissing real vector")
@n <- length(@a)-1;

if (@n <= 1) {
	@b <- @a
} else {
	@c <- vector(0, -1)
	@b <- @a[@n+1]
	for(@i,2,@n+1){
		@b <- movavg(@c,vector(@b, @a[@n+2-@i]))
		@b[run(@i-1)] <- @b[run(@i-1)]/2
	}
	@b <- reverse(@b)
	@b[-1] <- 2*@b[-1]
	delete(@c,@i)
}
delete(@a,@n)
delete(@b,return:T)
%chebcoefs%

invchebcoefs macro dollars
) Macro to find the coefficients of the expansion in Chebysev polynomials
) of a polynomial.
) Usage:
)   b <- invchebcoefs(a), where a = vector(a0,a1,...,an)
)     are coefficients of the polynomial Pn(x)=a0+T1(x)*x+...+an*Tn(x)
)   b is vector(b0,b1,...,bn), where Pn(x) = b0+b1*x+ ... + bn*x^n
)   Tj(x) is the jth Chebyshev polynomial.
) C. Bingham, kb@stat.umn.edu
)) 000607
@b <- argvalue($1,"arg 1","nonmissing real vector")
@n <- length(@b) - 1
if (@n > 1) {
	@b[-1] <- @b[-1]/2
	@b <- reverse(@b)
	@d <- vector(0, -1)
	@a <- rep(0,@n+1)
	@k <- 1
	for (@i, @n+1,1) {#print(paste("@i =",@i,"@k =",@k,"@n =",@n))
		@b0 <- @b[@i]
		@b <- autoreg(@d,2*@b)
		@a[@k] <- (@b[@i] - @b0)/2
		@b <- @b[-@i]
		@k <-+ 1
	}
	@b <- 2*@a
	delete(@d,@k,@i,@a)
}
delete(@b,return:T)
%invchebcoefs%

invertseries MACRO DOLLARS
) Macro to compute the a power series of the solution
) x = B(y) = sum(b[j]*y^j,j=1,.) where y satisifies
) A(x) =  = sum(a[j]*x^j,j=1,.)
) usage:
)  b <- invertseries(a), REAL vector a
)
)) Written by C. Bingham, kb@stat.umn.edu
)) Version 000901
@a <- argvalue($1,"coefficients","nonmissing real vector")
@a1 <- @a[1]
if (@a1 == 0) {
	error("a[1] = 0")
}

@n <- length(@a)
@a <- if (@n == 1) {
	1
} else {
	@a <-/ @a1
	@A <- @a * padto(1,@n)'
	for(@j,2,@n){
		@I <- run(@j-1 ,@n) # length @n - @j + 2
		@A[@I[-1],@j] <- movavg(-@a[-1], @A[@I,@j-1])[@I[-1] - @j + 1]
	}
	delete(@j,@I)
 #	@b <- padto(1,@n) # alternate method, not using solve
 #	for(@i,2,@n){
 #		@b[@i] <- -@A[@i,] %*% @b
 #	}
 #	delete(@b,return:T)
 	vector(solve(delete(@A,return:T),padto(1,@n)))
}/@a1^run(@n)
delete(@n,@a1)
delete(@a,return:T)
%invertseries%

factors         MACRO DOLLARS
) Macro to factor positive integers.
) Usage:
)  m <- factors(n)
)    n       vector of positive integers
)    m       vector of prime factors of scalar n or structure with
)            m[i] = vector of prime factors of n[i].  The component
)            names are either 'composite' or 'prime' when n is not
)            a scalar
)
) This is almost identical to function primefactors(), except, when n is
) a vector, it can handle n[i] >= 1000000000000.  Also the component
) names are either 'prime' or 'composite' instead of the value of n[i]
)
))000830 no longer uses changestr but assigns directly to components
))011216 now uses primefactors() to get factors
# usage: $S(n)  prime factors of elements of vector n of positive integers
@N <- vector(argvalue($1,"argument to $S","positive integer vector"))
@length <- length(@N)
if(@length > 1){
	@which <- rep(0,@length)
	@ans <- split(@which')

	for(@i,run(@length)){
		@tmp <- primefactors(@N[@i])
		@which[@i] <- isscalar(@tmp) + 1
		@ans[@i] <- @tmp
	}
	@ans <- strconcat(@ans,compnames:vector("composite","prime")[@which])
	delete(@tmp,@i,@which)
}else{
	@ans <- primefactors(@N)
}
delete(@length, @N)
delete(@ans,return:T)
%factors%

factors         MACRO DOLLARS
) Macro to factor positive integers.  factors(n) returns a vector of
) factors
))000830 no longer uses changestr but assigns directly to components
# usage: $S(n)  prime factors of elements of vector n of positive integers
@N <- vector(argvalue($1,"argument to $S","positive integer vector"))
@length <- length(@N)
if(@length > 1){
	@which <- rep(0,@length)
	@ans <- split(@which')

	for(@i,run(@length)){
		@tmp <- $S(@N[@i])
		@which[@i] <- isscalar(@tmp) + 1
		@ans[@i] <- @tmp
	}
	@ans <- strconcat(@ans,compnames:vector("composite","prime")[@which])
	delete(@tmp,@i,@which)
}else{
	@ans <- @N
	for(@Prime,vector(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,\
		67,71,73,79,83,89,97)){
		while(@N %% @Prime == 0){
			@ans <- vector(@ans, @Prime)
			@N <-/ @Prime
		}
	}
	@Prime <- 101
	@LIM <- sqrt(@N)
	while(@N > 1 && @Prime <= @LIM){
		for(@i,1,100){
			while(@N %% @Prime == 0){
				@ans <- vector(@ans, @Prime)
				@N <-/ @Prime
				@LIM <-/ sqrt(@Prime)
			}
			@Prime <-+ 2
			if(@Prime > @LIM || @N <= 1){
				break;
			}
		}
	}
	if(length(@ans) > 1){
		@ans <- @ans[-1]
		if(@N > 1){
			@ans <- vector(@ans,@N)
		}
	}
	delete(@Prime, @LIM)
}
delete(@length, @N)
delete(@ans,return:T)
%factors%

printfactors         MACRO DOLLARS
) macro to print out factors of all elements of an integer vector
) Usage:
)   printfactors(vector(n1,n2,...))
)   where n1, n2, ... are positive integers
)) version 990901
))000219 stripped $$
# usage: $S(n), prints factorization of elements of integer vector n
@vec <- argvalue($1,"argument 1","positive integer vector")
for(@n, @vec){
	@tmp <- primefactors(@n)
	if(isscalar(@tmp)){
		print(paste(@n,"is prime"))
	}else{
		print(paste(@n,"=",@tmp[1]+1e-1,sep:"*",@tmp[-1]+1e-1,format:".0f"))
	}
}
delete(@n,@vec,@tmp)
%printfactors%

partitions           MACRO DOLLARS
) Usage:
)   parts <- partitions(n)
)   allparts <- partitions(n, all:T)
)    n            positive integer
)    parts        matrix of non-negative integers with n columns with
)                 whose rows are the partitions of n (padded with 0),
)                 that is, each row sums to n
)    allparts     structure with n components, with allparts[m] the same
)                 as partitions(m), m = 1, 2, ..., n
)
)  The number N of rows of partitions(n) grows rapidly
)
)  n     N     n     N     n     N     n     N
)  1     1     7    15    13   101    19   490
)  2     2     8    22    14   135    20   627
)  3     3     9    30    15   176    21   792
)  4     5    10    42    16   231    22  1002
)  5     7    11    56    17   297    23  1255
)  6    11    12    77    18   385    24  1575
)
)) Version of 030512; written by C. Bingham, kb@umn.edu
#usage: parts <- $S(n [,all:T])
@N <- argvalue($1,"n","positive count")
@all <- keyvalue($K,"all","TF",default:F)
if (@N == 1){
	@R <- structure(Parts_of_1:1)
}else{
	@R <- split(rep(1,@N)',\
		compnames:getlabels(vector(rep(0,@N),labels:"Parts_of_")))
	for (@n, 2, @N){
		@X <- padto(@n,@n)
		for(@i,1,@n-1){
			@m <- @n - @i
			@r <- @R[@i]
			for(@j,run(ncols(@r))){
				if (@r[1,@j] <= @m){
					@X <- hconcat(@X,padto(vector(@m,@r[,@j]),@n))
				}
			}
		}
		@R[@n] <- @X
	}
	delete(@X,@n,@i,@j,@r,@m)
}
if (!delete(@all,return:T)) {
	@R <- @R[@N]
}
delete(@N)
delete(@R,return:T)'
%partitions%

===> cmatmultc <===
cmatmultc   macro  dollars
) Macro to compute matrix product A %*% B of complex matrices A and B
) Usage:
)  c <- cmatmultc(a, b [, "%*%"])
)  c <- cmatmultc(a, b, "%c%")
)  c <- cmatmultc(a, b, "%C%")
)    a, b              REAL matrices with no MISSING elements representing
)                      complex matrices A and B in fully complex form
)                      (Re(A) and Re(B) in odd columns of a and b, Im(A)
)                      and Im(B) in even columns, possibly after appending
)                      a column of zeros).  It is an error if
)                      floor((ncols(a) + 1)/2) != nrows(b)
)    c                 A REAL matrix representing the complex matrix product
)                      A %*% B, A %c% B or A %C% B in fully complex form
))030317 written by C. Bingham, kb@umn.edu
# c <- $S(a,b [,op]), a, b fully complex matrices
# op = "%*%" (default), "%c%", or "%C%"
@a <- argvalue($1,"left operand","nonmissing real matrix")
@b <- argvalue($2,"right operand","nonmissing real matrix")
@re_a <- creal(@a)
@re_b <- creal(@b)
if ($v < 3) {
	@op <- "%*%"
	@n1 <- ncols(@re_a)
	@n2 <- nrows(@re_b)
} else {
	@op <- argvalue($3,"operation", "string")
	if (alltrue(@op != "%*%",@op != "%c%",@op != "%C%")){
		error(paste("'",@op,"' not one of '%*%', '%c%' or '%C'",sep:""))
	}
	@n1 <- if (@op == "%c%") {
		nrows(@re_a)
	} else {
		ncols(@re_a)
	}
	@n2 <- if (@op == "%C%") {
		ncols(@re_b)
	} else {
		nrows(@re_b)
	}
}
if (delete(@n1,return:T) != delete(@n2,return:T)){
	error(paste("Dimension mismatch:",nrows(@re_a),"by",ncols(@re_a),\
		  @op,nrows(@re_b),"by",ncols(@re_b)))
}
@im_a <- cimag(delete(@a,return:T))
@im_b <- cimag(delete(@b,return:T))
if(!ismacro(@matmultit)){
  @matmultit <- macro(paste("cmplx(",nameof(@re_a),@op,nameof(@re_b),"-",\
	  nameof(@im_a),@op,nameof(@im_b),",",nameof(@re_a),@op,nameof(@im_b),"+",\
	  nameof(@im_a),@op,nameof(@re_b),")",sep:""))
}
# (re_a + i*im_a) OP (re_b + i*im_b)
#  = re_a OP re_b - im_a OP im_b + i*(re_a OP im_b + im_a OP re_b)
@result <- @matmultit()
delete(@re_a,@im_a,@re_b,@im_b)
delete(@result,return:T)
%cmatmultc%

===> ctranspose <===
ctranspose   macro dollars
) Macro to compute the transpose of a complex matrix
)  Usage:
)   b <- ctranspose(a)
)    a           REAL matrix representing a complex matrix A in fully
)                complex form (reals in odd columns, imaginaries in
)                even columns, possibly after appending a column of 0's
)    b           REAL matrix representing the complex matrix A' in fully
)                complex form
)) 030317 Written by C. BIngham, kb@umn.edu
# b <- $S(a), complex matrix a in fully complex form
@a <- argvalue($1, "argument", "real matrix")
@re <- creal(@a)'
cmplx(delete(@re,return:T),cimag(delete(@a,return:T))')
%ctranspose%

===> cjtranspose <===
cjtranspose   macro dollars
) Macro to compute the transpose of the complex conjugate of a
) complex matrix
)  Usage:
)   b <- cjtranspose(a)
)    a           REAL matrix representing a complex matrix A in fully
)                complex form (reals in odd columns, imaginaries in
)                even columns, possibly after appending a column of 0's
)    b           REAL matrix representing the complex matrix conj(A') in
)                fully complex form
)) 030317 Written by C. BIngham, kb@umn.edu
# b <- $S(a), complex matrix a in fully complex form
if (!ismatrix(ctranspose)){
	getmacros(ctranspose,silent:T)
}
ctranspose(cconj(argvalue($1, "argument", "real matrix")))
%cjtranspose%

===> ceigen <===
ceigen     macro dollars
) Macro to compute the real eigenvalues and complex eigenvectors of a
) a Hermitian symmetric matrix (ctranspose(cconj(a)) = a)
)  Usage:
)   result <- ceigen(a)
)    a                REAL matrix representing a square fully complex
)                     Hermitian matrix A (A' = conj(A)).  ncols(a)
)                     must be 2*n or 2*n-1, where n = nrows(a)
)    result           structure(values:vals, vectors:vecs)
)                     vals, length n vector of eigenvalues
)                     vecs, n by 2*n matrix, with columns 2*i-1 and
)                     2*i containing the real and imaginary parts of
)                     complex eigenvectors of A
)   Caution: with equal eigenvalues, some vectors may be linearly
)   dependent
))030317 written by C. Bingham, kb@umn.edu
# str <- $S(a)
@a <- matrix(argvalue($1,"argument","nonmissing real matrix"))
@n <- nrows(@a)
if(@n != floor((ncols(@a) + 1)/2)){
	error("Argument is not square complex matrix")
}
@R <- creal(@a)
@fuzz <- sum(vector(abs(@a)))/prod(dim(@a))*1e-7
if (@fuzz > 0){
	@symmetric <- max(abs(vector(@R-@R'))) <= @fuzz
	if (@symmetric){
		@I <- cimag(@a)
		@symmetric <- max(abs(vector(@I + @I'))) <= @fuzz
	}
	if (!delete(@symmetric,return:T)){
		error("Argument is not Hermitian")
	}
	@A <- vconcat(hconcat(@R,-@I),hconcat(@I,@R))
	delete(@R,@I)
} else {
	@A <- matrix(rep(0,4*@n^2),2*@n)
	delete(@R)
}
delete(@a,@fuzz)
@eigs <- eigen(delete(@A,return:T))
@J <- run(1,2*@n,2)
@eigs <- structure(values:@eigs$values[@J],\
			vectors:matrix(@eigs$vectors[,@J],delete(@n,return:T)))
delete(@J)
delete(@eigs,return:T)
%ceigen%

ctrace    macro    dollars
) Macro to find the trace of a complex matrix in fully complex form
) Usage:
)  b <- ctrace(a)
)   a              REAL matrix interpreted as a square complex matrix
)                  A in fully complex form.  This imples that
)                  floor((ncols(a)+1)/2) = nrows(a) is required.
)   b              REAL 1 by 2 matrix representing the complex scalar
)                  tr(A) = trace(creal(a)) + i*trace(cimag(a))
))030317 written by C. Bingham, kb@umn.edu
#b <- $S(a), REAL matrix interpreted as a square complex matrix in fully
#            complex form
@a <- argvalue($1,"argument","nonmissing real matrix")
@n <- nrows(@a)
if (@n != floor((ncols(@a)+1)/2)){
	error("argument not square complex matrix in fully complex form")
}
@s <- cmplx(trace(creal(@a)),trace(cimag(@a)))
delete(@n,@a)
delete(@s,return:T)
%ctrace%

===> csolve <===
csolve     macro  dollars
) Macro to invert a complex matrix in fully complex form
) Usage:
)  ainv <- csolve(a)
)   a           REAL matrix with no MISSING elements representing
)               a square complex matrix A in fully complex form.
)   ainv        REAL matrix representing a square complex matrix
)               Ainv in fully complex form satisfying Ainv %*% A
)               is the identity matrix.
) It is mathematically possible but unlikely that a will be considered
) singular when it is not
) This may require function cprdc() which was introduced with
) release 2 of version 4.13 (March 2003)
))030317 written by C. Bingham, kb@umn.edu
# ainv <- csolve(a)
@a <- argvalue($1,"argument","nonmissing real matrix")
@n <- nrows(@a)
if (@n != floor((ncols(@a)+1)/2)){
	error("argument not square complex matrix in fully complex form")
}
@re <- creal(@a)
@im <- cimag(@a)
@q <- solve(@re,@im,singok:T)
@usez <- isnull(@q)
if (@usez){
	# @z is arbitrary, in hopes Re(cprdc(@z,@a)) non-singular
	if (!isfunction(cprdc)){
		error("$S() won't work in this version of MacAnova", macroname:F)
	}
	@z <- vector(0.5460400862681714, 3.775912244316442)'
	@a <- cprdc(@z,@a)
	@re <- creal(@a)
	@im <- cimag(@a)
	@q <- solve(@re,@im,singok:T)
	if (isnull(@q)){
		error("Argument apparently singular")
	}
}
@re <- solve(@re + @im %*% @q)
@im <- delete(@q,return:T) %*% @re
@result <- cmplx(delete(@re,return:T),-delete(@im,return:T))
if (delete(@usez,return:T)){
	@result <- cprdc(delete(@z,return:T), @result)
}
delete(@result,return:T)
%csolve%

csubscr    macro   dollars
) Macro to simulate subscripted complex vector or matrix
)  Usage:
)   cy <- csubscr(cx,i,j)
)    i, j     REAL or LOGICAL vectors that would be legal subscripts
)             for m by n matrix
)    cy       REAL matrix containing complex elements X[i,j].
)
)   cy <- csubscr(cx,i)
)    i        REAL vector or matrix or LOGICAL vector that would be
)             legal subscript for X[i].
)    cy       REAL matrix containing complex elements X[i].
)
)   cy <- csubscr(cx)
)    cy       REAL matrix cmplx(creal(cx),cimag(cx))
)
)   An empty or NULL i or j is treated as an empty subscript.
))030317 written by C. Bingham, kb@umn.edu
@a <- argvalue($1,"argument","real matrix")
@real <- creal(@a)
@imag <- cimag(@a)
@ndims <- if (ncols(@a) <= 2){
	1
} else {
	2
}

@result <- if ($v == 1){
	cmplx(@real,@imag)
} else {
	@sub1 <- argvalue($02, "subscript")
	if (@ndims == 1){
		if (isvector(@sub1)) {
			cmplx(@real[@sub1],@imag[@sub1])
		}else {
			error("subscript not vector")
		}
	} elseif (@ndims == 2) {
		if ($v == 2){
			if (isvector(@sub1)){
				cmplx(@real[@sub1,],@imag[@sub1,])
			} elseif (ismatrix(@sub1)) {
				cmplx(@real[@sub1],@imag[@sub1])
			} else {
				error("single subscript not vector or matrix")
			}
		} else {
			@sub2 <- argvalue($03, "subscript 2")
			if (isnull(@sub1)){
				if (isnull(@sub2)){
					cmplx(@real,@imag)
				} elseif (isvector(@sub2)) {
					cmplx(@real[,@sub2],@imag[,@sub2])
				} else {
					error("non-vector subscript 2")
				}
			} elseif (isvector(@sub1)){
				if (isnull(@sub2)){
					cmplx(@real[@sub1,],@imag[@sub1,])
				} elseif (isvector(@sub2)) {
					cmplx(@real[@sub1,@sub2],@imag[@sub1,@sub2])
				} else {
					error("non-vector subscript 2")
				}
			} else {
				error("non-vector subscript 1")
			}
		}
	}
}
delete(@real,@imag,@ndims,@sub1,@sub2,silent:T)
delete(@result,return:T)
%csubscr%

cdiag    macro    dollars
) Macro to extract the diagonal of a square complex matrix stored in
) fully complex form
) Usage:
)  d <- cdiag(a)
)   a         REAL matrix representing square complex matrix A in fully
)             complex form (ncols(a) = 2*nrows(a) or 2*nrows(a) - 1
)   d         REAL nrows(a) by 2 matrix representing complex vector
)             diag(A) in fully complex form
)) 030318 written by C. Bingham (kb@umn.edu)
# d <- cdiag(a), REAL matrix a representing square complex matrix A
@a <- argvalue($1,"argument", "real matrix")
@m <- nrows(@a)
if (@m != floor((ncols(@a) + 1)/2)){
	error("argument is not square complex matrix")
}
delete(@m)
@a <- cmplx(diag(creal(@a)),diag(cimag(@a)))
delete(@a,return:T)
%cdiag%

mathhelp   MACRO
) Macro to get help on macros in math.mac
# usage $S(topic1 [, topic2 ...] [help keywords])
if(!ismacro(_gethelp)){
	getmacros(_gethelp,silent:T)
}
___HELPFILE_ <- "math.mac"
___INDEXTOP_ <- "math_index"
___MACRO_ <- "$S"
_gethelp($0)
%mathhelp%

_E_N_D_O_F_M_A_C_R_O_S_

Help file for math.mac for MacAnova
(C) 2003 by Gary W. Oehlert and Christopher Bingham
Updated 030814 CB

!!!! Starting marker for message of the day
!!!! Ending marker for message of the day

???? Starting marker for list of up to 32 comma/newline separated keys
Bessel functions
Binomial coefficients
Complex matrices
Continued fractions
Direct search
Expansions
General
Generalized inverse
Integers
Matrices
Minimize
Nonlinear least squares
Orthogonal polynomials
Power series
Prime factors
Quasi-newton
QR decomposition
Special functions
Variable metric
???? Ending marker for keys

011007 Minor changes to help
011118 Update help to match changes to bfs(), dfp() and broyden()
011216 Update help to match changes to neldermead()
       Added help keys

====bfs()#minimize,quasi-newton,variable metric
%%%%
bfs(x0, fun [,params:params]  [, goldsteps:ngold] [, maxit:maxiter]\
  [,minit:miniter] [,criteria:vector(nsigx,nsigfun,dgrad)]\
  [printwhen:d1] [,recordwhen:d2]), REAL vector x0, macro
  fun(x,i [,params]), integers ngold > 0, maxiter >= 0, miniter > 0,
  nsigx, nsigfun, d1 >= 0, d2 >= 0, dgrad REAL scalar
%%%%
@@@@introduction#Introduction
Macro bfs() uses the Broyden-Fletcher-Shanno variable metric algorithm
to minimize a function iteratively.  A golden mean line search is made
at each step.  See Dahlquist and Bjorck, Numerical methods, Prentice
Hall, 1974, p. 443.

bfs() is a "front-end" to macro minimizer() which it calls with
all the arguments to bfs() plus argument 'method:"bfs"'.

@@@@usage#Usage
result <- bfs(x0, fun [, params] [,optional keywords]) computes the
minimum of a real function F(x1,x2,...,xk) starting the Broyden-
Fletcher-Shanno iteration at x = x0 = vector(x01,x02,...,x0k),
a REAL vector with no MISSING elements.

See minimizer() for details on the arguments, keywords and the value.

@@@@see_also#Cross references
See also mnimizer(), dfp(), broyden(), and neldermead()
@@@@______

====binom()#binomial coefficients,integers
%%%%
binom(n,k), REAL n and k with non negative elements.  If both are non-
  scalars, they must have the same dimensions
%%%%
@@@@usage#Usage
binom(n,k), where n >= 0 and k >= 0 are REAL scalars with n >= k returns
a binomial coefficient.  n and k need not be integers, but when they are
binom(n,k) returns n!/(k!*(n-k)!) as an exact integer.  Otherwise it
returns gamma(n+1)/(gamma(k+1)*gamma(n-k+1)), where gamma(x) is the
Gamma function.

When just one of n and k is a scalar, binom(n,k) returns a REAL vector,
matrix or array consisting of binomial coefficients computed from the
scalar and each of the elements of the other argument.

When neither n or k is a scalar, both must have exactly the same
dimensions and the result is an array with the same dimesions
consisting of of binomial coefficients computed from corresponding
elements of n and k.

@@@@examples
Examples:
  Cmd> binom(4,run(0,4)) # vector(binom(4,0),...,binom(4,4))
  (1)          1          4          6          4          1

  Cmd> binom(run(3,7),3) # vector(binom(3,3),...,binom(7,3))
  (1)          1          4         10         20         35

  Cmd> binom(run(3,7),run(0,4)) # vector(binom(3,0),...,binom(7,4))
  (1)          1          4         10         20         35

@@@@see_also#Cross reference
See also lgamma().
@@@@______

====blockdmat()#matrices
%%%%
blockdmat(A1,A2,...,Ak), A1, A2, ... matrices, all of the same type,
  REAL, LOGICAL or CHARACTER
%%%%
@@@@usage#Usage
B <- blockdmat(A1,A2,...,Ak), where A1, A2, ..., Ak are matrices,
creates a block diagonal matrix B with diagonal blocks A1, ..., Ak.  B
will have the same type as A1, ..., Ak which must all have the same
type.  Elements of B outside the blocks are 0 , F or "", depending on
the type.

If Aj is mj by nj, j = 1,...,k, then B is m1 + m2 + ... + mk by n1 + n2
+ ... + nk.

@@@@example
Example:                                                 [1 1 1 0 0]
           [1 1 1]       [2 2]                           [1 1 1 0 0]
When  A1 = [1 1 1], A2 = [2,2], then blockdmat(A1, A2) = [0 0 0 2 2]
                         [2,2]                           [0 0 0 2 2]
                                                         [0 0 0 2 2]

@@@@see_also#Cross references
See also dmat(), diag().
@@@@______

====broyden()#minimize,quasi-newton,variable metric
%%%%
broyden(x0, fun [,params:params] [, maxit:maxiter] [,minit:miniter]\
  [,criteria:vector(nsigx,nsigfun,dgrad)] [printwhen:d1]\
  [,recordwhen:d2]), REAL vector x0, macro fun(x,i [,params]), integers
  ngold > 0, maxiter >= 0, miniter > 0, nsigx, nsigfun, d1 >= 0, d2 >=
  0, dgrad REAL scalar
%%%%
@@@@introduction#Introduction
Macro broyden() minimizes a function iteratively using a variable metric
algorithm due to Broyden.  It has no linear search step.  See Dahlquist
and Bjorck, Numerical methods, Prentice Hall, 1974, p. 443.

broyden() is a "front-end" to macro minimizer() which it calls with
all the arguments to minimizer() plus argument 'method:"broyden"'.

@@@@usage#Usage
result <- broyden(x0, fun [, params] [,optional keywords]) computes the
minimum of a real function F(x1,x2,...,xk) starting the Broyden
iteration at x = x0 = vector(x01,x02,...,x0k), a REAL vector with no
MISSING elements.

See minimizer() for details on the arguments, keywords and the value.
Keyword 'golden' is ignored by broyden().

@@@@see_also#Cross references
See also minimizer() bfs(), dfp(), and neldermead().
@@@@______

====cdiag()#complex matrices,matrices
%%%%
d <- cdiag(a), REAL matrix a representing the fully complex form of a
  square complex matrix A
%%%%
@@@@usage#Usage
d <- cdiag(a), where a is a REAL matrix representing a square complex
matrix A in fully complex form.

d is a nrows(a) by 2 REAL matrix representing diag(A) in fully complex
form.

It is an error for A to not be square, that is for ncols(a) != 2*nrows(a)
and ncols(a) != 2*nrows(a) - 1.

@@@@see_also#Cross references
See also csubscr(), diag(), 'complex'.

====ceigen()#complex matrices,matrices
%%%%
eigs <- ceigen(a), REAL matrix a representing the fully complex form of
  a complex Hermitian matrix A (A' = conj(A)).  Result is
  structure(values:vals, vectors:vecs)
%%%%
@@@@usage#Usage
result <- ceigen(a) computes the real eigenvalues and complex eigen-
vectors of a complex matrix A with Hermitian symmetry (A' = conj(A)),
coded in fully complex form in REAL matrix a with no MISSING elements.

result is structure(values:V, vectors:U).  V is a length n vector of
eigen- values where n = nrows(a).  U is a n by 2*n REAL matrix; columns
2*i-1 and 2*i contain the real and imaginary parts of the i-th complex
eigenvector of A.

@@@@caution#Caution
When A has duplicate eigenvalues, some of the eigenvectors computed may
be linearly dependent.

@@@@example
Example:
  Cmd> a <- cmplx(matrix(vector(8,2,2,1),2),matrix(vector(0,-3,3,0),2))

  Cmd> a # 2 by 2 Hermitian symmetric complex matrix
  (1,1)            8            0            2            3
  (2,1)            2           -3            1            0

  Cmd> eigs <- ceigen(a)

  Cmd> eigs
  component: values
  (1)       9.5249     -0.52494
  component: vectors
  (1,1)     -0.70598     -0.59149     -0.31138     -0.23405
  (2,1)     -0.37378      0.10967      0.86883     -0.30561

  Cmd> cdivc(cmatmultc(a,eigs$vectors),eigs$vectors)
  (1,1)       9.5249   6.4749e-16     -0.52494  -2.4739e-16
  (2,1)       9.5249  -4.0029e-16     -0.52494   2.2808e-16

@@@@see_also#Cross references
See also cdivc(), cmatmultc(), eigen(), 'complex'.

====chebcoefs()#expansions,power series,orthogonal polynomials
%%%%
chebcoefs(vector(a0,a1,...,an)), a0, a1, ..., an non MISSING REAL
  scalars
%%%%
@@@@usage#Usage
Suppose Pn(x) = a0+a1*x+ ... + an*x^n is a polynomial in x of degree
n.  Then Pn(x) has a unique expansion, Pn(x) = b0 + b1*T1(x) + ... +
bn*Tn(x), where Tj(x) = the Chebyshev polynomial of degree j defined
for -1 <= x <= 1 as Tj(x) = cos(j*acos(x)).

chebcoefs(vector(a0,a1,...,an)), where a0, a1, ..., an are non MISSING
REAL scalars returns vector(b0,b1,...,bn), where the b's are the
coefficients in the Chebyshev expansion.

@@@@see_also#Cross reference
See also topic invchebcoefs().
@@@@______

====cjtranspose()#complex matrices,matrices
%%%%
b <- cjtranspose(a), REAL matrix a representing a complex matrix A in
  fully complex form.  b is the transpose conj(A)' in fully complex
  form.
%%%%
@@@@usage#Usage
b <- cjtranspose(a) is equivalent to b <- ctranspose(cconj(a)) and
computes the transpose of the complex conjugate of complex matrix a in
fully complex form.

@@@@see_also#Cross references
See also ctranspose(), transpose(), cconj(), 'complex'.

====cmatmultc()#complex matrices,matrices
%%%%
c <- cmatmultc(a, b), a and b REAL matrices representing complex
  matrices A and B in fully complex form, with nrows(b) =
  floor((ncols(a) + 1)/2)
c <- cmatmultc(a, b, "%c%"), requires nrows(a) = nrows(b)
c <- cmatmultc(a, b, "%C%"), requires floor(ncols(a) + 1)/2) =
  floor(ncols(b) + 1)/2
%%%%
@@@@usage#Usage
c <- cmatmultc(a, b), computes the matrix product of REAL matrices a and
b interpreted as complex matrices A and B in fully complex form (real
parts in odd columns, imaginary parts in even).

The result c is a REAL matrix interpreted as the complex matrix A %*% B,
in fully complex form.

c <- matmultc(a, b, "%*%") does the same.

It is required that nrows(b) = floor((ncols(a) + 1)/2).

c <- cmatmultc(a, b, "%c%") does the same except the matrix product
is A %c% B = A' %*% B and nrows(a) = nrows(b) is required.

c <- cmatmultc(a, b, "%C%") does the same except the matrix product
is A %C% A = A %*% B' and floor((ncols(a)+1)/2) = floor((ncols(b)+1)/2)
is required.

@@@@see_also#Cross references
See also 'matrices', 'complex'.

====continfrac()#continued fractions,special functions
%%%%
continfrac(a, b), a and b REAL vectors or matrices with nrows(b) =
  nrows(a) or nrows(b) = nrows(a) + 1
%%%%
@@@@usage#Usage
continfrac(a,b), where a and b are REAL vectors with no MISSING
elements, and nrows(a) = nrows(b) = m, evaluates the continued
continued fraction a[1]/(b[1] + a[2]/(b[2] + a[3]/(b[3] + a[4]/... +
a[m]/b[m])).

continfrac(a,vector(b0, b)) returns b0 + continfrac(a,b), when b0 is a
non MISSING real scalar.

a and b can also be matrices with the same shape.  In that case,
continfrac(a,b) returns hconcat(continfrac(a[,1], b[,1]),...,
continfrac(a[,m], b[,m]).

If either a or b is a vectir and the other has more than 1 column, the
vector is used in each column of the result.

continfrac(a,vconcat(b0',b)) returns b0 + continfact(a,b), where b0 is
is a vector of length ncols(b).
@@@@______

====csolve()#complex matrices,matrices
%%%%
ainv <- csolve(a), REAL matrix a interpreted as a square complex matrix
  in fully complex form
%%%%
@@@@usage#Usage
ainv <- csolve(a) computes the complex inverse of REAL matrix a,
interpreted as a square complex matrix A in fully coplex form.

It is an error if A is singular.

@@@@caution#Caution
It is possible but unlikely that csolve() will report that A is singular
when that is not the case.

@@@@example#Example
  Cmd> areal <- matrix(vector(0.57,-0.24,-0.33,-0.55),2)

  Cmd> aimag <- matrix(vector(0.43,-0.08,-0.16,0.26),2)

  Cmd> a <- cmplx(areal,aimag)

  Cmd> ainv <- csolve(a)

  Cmd> prd <- cmatmultc(ainv,a)

  Cmd> creal(prd)
  (1,1)            1  -2.7756e-17
  (2,1)   3.4694e-17            1

  Cmd> cimag(prd)
  (1,1)            0            0
  (2,1)   6.9389e-17   5.5511e-17

@@@@see_also#Cross references
See also cmplx(), cmatmultc(), creal(), cimag(), solve(), 'matrices',
'complex'.

====csubscr()#complex matrices,matrices
%%%%
y <- csubscr(x,I), REAL matrix x containing complex matrix X in fully
  complex form, legal REAL or LOGICAL subscript or NULL I
y <- csubscr(x,I,J), legal REAL or LOGICAL subscripts or NULL, I and J
%%%%
@@@@usage#Usage
csubscr() simulates subscript extraction from a complex vector X or m by
n complex matrix X stored in REAL matrix x in fully complex form (real
parts in odd columns, imaginary parts in even columns).  The result y
contains a complex vector or matrix Y in fully complex form.

y <- csubscr(x,I) simulates Y <- X[I,] when I is a vector or Y <- X[I]
when I is a matrix with two columns.  When I is empty or NULL, it
simulates Y <- X[].

y <- csubscr(x,I,J) simulates Y <- X[I,J].  When I is empty or NULL it
simulates Y <- X[,J]; when J is empty or NULL it simulates Y <- X[I,].

y <- csubscr(x) is equivalent to cmplx(creal(x),cimag(y)).

@@@@caution#Caution
Unlike real subscripts you cannot assign to csubscr(x, I, J).

@@@@example#Example
Example with 2 by 4 REAL a representing 2 by 2 complex A:
  Cmd> creal(a) # real part
  (1,1)        0.57       -0.33
  (2,1)       -0.24       -0.55

  Cmd> cimag(a) # imaginary part
  (1,1)        0.43       -0.16
  (2,1)       -0.08        0.26

  Cmd> csubscr(a,1,2) # simulates A[1,2]
  (1,1)       -0.33       -0.16

  Cmd> csubscr(a,hconcat(run(2),run(2))) # simulates cdiag(a)
  (1,1)        0.57        0.43
  (2,1)       -0.55        0.26

  Cmd> csubscr(a,,-1) # or csubscr(a,NULL,-1), simulates A[,-1]
  (1,1)       -0.33       -0.16
  (2,1)       -0.55        0.26

@@@@see_also#Cross references
See also cdiag(), 'complex', 'subscripts', creal(), cimag().

====ctrace()#complex matrices,matrices
%%%%
b <- ctrace(a), REAL matrix a interpreted as a square complex matrix in
  fully complex form.
%%%%
@@@@usage#Usage
c <- ctrace(a) computes the complex trace of REAL matrix a, interpreted
as a square complex matrix A in fully complex form.  c has value
cmplx(trace(creal(a)),trace(cimag(a))).  For A to be squre, nrows(a)
must be floor((ncols(a)+1)/2).'

@@@@example
Example:
  Cmd> areal <- matrix(vector(0.57,-0.24,-0.33,-0.55),2)

  Cmd> aimag <- matrix(vector(0.43,-0.08,-0.16,0.26),2)

  Cmd> ctrace(cmplx(areal,aimag))
  (1,1)        0.02        0.69

  Cmd> cmplx(trace(areal),trace(aimag)) # check
  (1,1)        0.02        0.69

@@@@see_also#Cross references
See also trace(), cmplx(), 'complex'.

====ctranspose()#complex matrices,matrices
%%%%
b <- ctranspose(a), REAL matrix a representing a complex matrix in fully
  complex form
%%%%
b <- ctranspose(a) computes the complex transpose of the complex matrix
a in fully complex form.  That is creal(b) = creal(a)' and cimag(b) =
cimag(a)'.

@@@@example#Example
  Cmd> areal <- matrix(vector(0.57,-0.24,-0.33,-0.55),2)

  Cmd> aimag <- matrix(vector(0.43,-0.08,-0.16,0.26),2)

  Cmd> a <- cmplx(areal,aimag)

  Cmd> atrans <- ctranspose(a)

  Cmd> creal(atrans) # transpose of areal
  (1,1)        0.57       -0.24
  (2,1)       -0.33       -0.55

  Cmd> cimag(atrans) # transpose of aimag
  (1,1)        0.43       -0.08
  (2,1)       -0.16        0.26

@@@@see_also#Cross references
See also cjtranspose(), 'complex', 'transpose'.

====dfp()#minimize,quasi-newton,variable metric
%%%%
dfp(x0, fun [,params:params]  [, goldsteps:ngold] [, maxit:maxiter]\
  [,minit:miniter] [,criteria:vector(nsigx,nsigfun,dgrad)]\
  [printwhen:d1] [,recordwhen:d2]), REAL vector x0, macro
  fun(x,i [,params]), integers ngold > 0, maxiter >= 0, miniter > 0,
  nsigx, nsigfun, d1 >= 0, d2 >= 0, dgrad REAL scalar
%%%%
@@@@introduction#Introduction
Macro dfp() uses the Davidon-Fletcher-Powell variable metric algorithm
to minimize a function iteratively.  A golden mean line search is made
at each step.  See Dahlquist and Bjorck, Numerical methods, Prentice
Hall, 1974, p. 442.

dfp() is a "front-end" to macro minimizer() which it calls with
all the arguments to bfs() plus argument 'method:"dfp"'.

@@@@usage#Usage
result <- dfp(x0, fun [, params] [,optional keywords]) computes the
minimum of a real function F(x1,x2,...,xk) starting the Davidon-
Fletcher-Powell iteration at x = x0 = vector(x01,x02,...,x0k),
a REAL vector with no MISSING elements.

See minimizer() for details on the arguments, keywords and the value.

@@@@see_also#Cross references
See also minimizer(), bfs(), broyden(), and neldermead()
@@@@______

====economize()#expansions,power series
%%%%
economize(vector(a0,a1,...,an),m), a0, ..., an REAL scalars, m > 0 an
  integer scalar
%%%%
@@@@introduction#Introduction
Macro to "economize" a power series expansion on [-1, 1]

Suppose Fn(x) = a0 +a1*x + a2*x^2 + ... + an*x^n is the n-th partial
sum of the Taylor series expansion of a function F(x) defined on an
interval I contained in (-1, 1).  Then Fn(x) has a unique expansion
         Fn(x) = c0 + c1*T1(x) + c2*T2(x) + ... + cn*Tn(x)
where Tj(x) is the j-th Chebyshev polynomial defined as Tj(x) =
cos(acos(j*x)) for -1 <= x <= 1.

@@@@usage#Usage
b <- economize(vector(a0,a1,...,an),m) computes b = vector(b0,b1,...,bm)
where
  Bm(x) = b0 + b1*x + b2*x^2 + ... + bm*x^m =
                c0 + c1*T1(x) + c2*T2(x) + ... + cm*Tm(x)
is the Chebyshev series truncated at Tm.

When Fn(x) is a good approximation to F(x) on I, and, as is the case
with many functions, the coefficients c0, c1, ... converge to 0 much
faster than do a0, a1, ..., Bm(x) with m << n, may be a good
polynomial approximation for F(x) on I.
@@@@______

====factorial()#special functions,integers
%%%%
factorial(x), x REAL
%%%%
@@@@usage#Usage
factorial(x) computes x! (x factorial), where x is a REAL scalar,
vector, matrix or array.  The result is the same size and shape as
as x.  If any element is MISSING, <= -1 or such that x! is too large
to be computed, the corresponding element of the result is MISSING.

Elements of x need not be integers.  x! is computed as exp(lgamma(x+1)),
except that if x is an integer <= 20 the value should be exact.

@@@@see_also#Cross references
See also binom() and lgamma().
@@@@______

====factors()#prime factors
%%%%
factors(vector(n1 [, n2, ...]), positive integer n1, n2, ...
%%%%
@@@@usage#Usage
factors(n), where n > 0 is an integer returns a REAL scalar or vector
containing the prime factors of n.

result <- factors(N), where N = vector(n1,n2,...), integers n1 > 0, n2 >
0, ..., sets result to structure(factor(n1), factor(n2), ...), that is
result[j] contains the factors of N[j].  The name of component j is
'prime' if N[j] is prime and 'composite' otherwise.

Note: factors() has been superceded by function primefactors().

@@@@example
Example:
  Cmd> write(factors(2^2^run(5) + 1),format:"10.0f")
  STRUCTURE:
  component: prime
  (1)          5
  component: prime
  (1)         17
  component: prime
  (1)        257
  component: prime
  (1)      65537
  component: composite
  (1)        641    6700417

@@@@see_also#Cross references
See also printfactors() and primefactors().
@@@@______

====i0()#bessel functions,special functions
%%%%
i0(x), x a REAL scalar, vector, matrix, array or structure of REAL
  components.
%%%%
@@@@usage#Usage
i0(x) computes values of I-sub-0(x), the modified Bessel function of the
first kind, where x is a REAL scalar, vector, matrix or array.  The
result is REAL of the same shape as x.

When x1, x2, ... are REAL, i0(structure(x1,x2,...)) returns
structure(i0(x1),i0(x2), ...)

i0 is based on equations 9.8.1 and 9.8.2 in Abramowitz & Stegun,
Handbook of Mathematical functions.

@@@@see_also#Cross reference
See also i1().
@@@@______

====i1()#bessel functions,special functions
%%%%
i1(x), x a REAL scalar, vector, matrix, array or structure of REAL
  components.
%%%%
@@@@usage#Usage
i1(x) computes values of I-sub-1(x), the modified Bessel function of the
first kind, where x is a REAL scalar, vector, matrix or array.  The
result is REAL of the same shape as x.

When x1, x2, ... are REAL, i1(structure(x1,x2,...)) returns
structure(i1(x1),i1(x2), ...)

i1 is based on equations 9.8.3 and 9.8.4 in Abramowitz & Stegun,
Handbook of Mathematical functions.

@@@@see_also#Cross reference
See also i0().
@@@@______

====invchebcoefs()#expansions,power series
%%%%
invchebcoefs(vector(b0,b1,...,bn)), b0, b1, ... non MISSING REAL
  scalars
%%%%
@@@@usage#Usage
Suppose Pn(x) = b0+b1*T1(x)+...+bn*Tn(x), where Tj(x) is the Chebyshev
polynomial of degree j, defined as Tj(x) = cos(j*acos(x)) for -1 <= x
<= 1.  Then Pn(x) is a polynomial and can be expressed as Pn(x) =
a0+a1*x+ ... + an*x^n.

invchebcoefs(vector(b0,b1,...,bn)), where b0, ..., bn are non MISSING
real scalars, returns vector(a0,a1,...,an), the coefficients of the
powers of x of Pn(x).

@@@@see_also#Cross reference
See also topic chebcoefs().
@@@@______

====invertseries()#expansions,power series
%%%%
invertseries(a), REAL vector a with non-MISSING elements
%%%%
@@@@usage#Usage
Suppose A(x) is a function with A(0) = 0 and with Taylor series a1*x +
a2*x^2 + a3^x^3 + ..., where a1 != 0.  Then in a neighborhood of 0
there is an inverse function B(y) is defined such that B(A(x)) = x for
y near 0.  Let b1*y + b2*y^2 + b3*y^3 ... be the Taylor series for
B(y).  The coefficients bj are unique functions of a1, a2, ..., aj; in
particular, b1 = 1/a1.

b <- invertseries(a), where a = vector(a1,a2,...,an) is a REAL vector,
computes the coefficients b = vector(b1,b2,...,bn) of the Taylor series
for B(y).

Because of numerical instability, higher order coefficients in b may
not be accurate.

@@@@example
Example:
Let A(x) = log(1-x) = -x - x^2/2 - x^3/3 - ... . Then B(y) = 1 - exp(y)
= -y - y^2/2 - y^3/6 - y^4/24 - y^5/120 - ... satisfies B(A(x)) = 1 -
exp(log(1 - x)) = 1 - (1 - x) = x.

  Cmd> a <- -1/run(10); a # coefficients of -log(1-x)
  (1)          -1        -0.5    -0.33333       -0.25        -0.2
  (6)    -0.16667    -0.14286      -0.125    -0.11111        -0.1

  Cmd> b <- invertseries(a); b # coefficients of 1 - exp(y)
  (1)          -1        -0.5    -0.16667   -0.041667  -0.0083333
  (6)  -0.0013889 -0.00019841 -2.4802e-05 -2.7557e-06 -2.7557e-07
@@@@______

====kronecker()#matrices
%%%%
kronecker(A,B), A, B REAL matrices with no MISSING values
%%%%
@@@@usage#Usage
C <- kronecker(A, B), where A and B are REAL matrices with no missing
values, computes the Kronecker product of A and B.

C is a nrows(A)*nrows(B) by ncols(A)*ncols(B) made of nrows(a)*ncols(a)
blocks of the form a[i,j]*B.

@@@@see_also#Cross reference
See also topic 'matrices'.
@@@@______

====mathhelp()#general
%%%%
mathhelp(topic1 [, topic2 ...] [,usage:T] [,scrollback:T])
mathhelp(topic, subtopic:Subtopics), CHARACTER scalar or vector
  Subtopics
mathhelp(topic1:Subtopics1 [,topic2:Subtopics2 ...])
mathhelp(key:Key), CHARACTER scalar Key
mathhelp(index:T [,scrollback:T])
%%%%
@@@@usage#Usage
mathhelp(Topic1 [, Topic2, ...]) prints help on topics Topic1, Topic2,
... related to macros in file math.mac.  The help is taken from file
math.mac.

mathhelp(Topic1 [, Topic2, ...] , usage:T) prints usage information
related to these macros.

mathhelp(index:T) or simply mathhelp() prints an index of the topics
available using mathhelp().

mathhelp(Topic, subtopic:Subtopic), where Subtopic is a CHARACTER scalar
or vector, prints subtopics of topic Topic.  With subtopic:"?", a list
of subtopics is printed.

mathhelp(Topic1:Subtopics1 [,Topic2:Subtopics2], ...), where Suptopics1
and Subtopics2 are CHARACTER scalars or vectors, prints the specified
subtopics.  You can't use any other keywords with this usage.

In all the first 4 of these usages, you can also include help() keyword
phrase 'scrollback:T' as an argument to mathhelp().  In windowed
versions, this directs the output/command window will be automatically
scrolled back to the start of the help output.

@@@@key_keyword#Keyword key
mathhelp(key:key) where key is a quoted string or CHARACTER scalar lists
all topics cross referenced under Key.  mathhelp(key:"?") prints a list
of available cross reference keys for topics in the file.
@@@@______

mathhelp() is implemented as a predefined macro.

@@@@see_also#Cross reference
See help() for information on direct use of help() to retrieve
information from math.mac.
@@@@______

====math_index*#general
%%%%
Topics in this file:
 Macros related to matrices
   blockdmat(), kronecker(), matsqrt(), moorepenrose(), qrdcomp()
 Macros related to complex matrices
   cdiag(), ceigen(), cmatmultc(), cjtranspose(), csolve(), csubscr(),
   ctrace(), ctranspose()
 Macros related to optimization
   dfp(), bfs(), broyden(), levmar(), _cgrad(), _lmout(), neldermead()
 Macros related to polynomials and series in powers of x
   chebcoefs(), economize(), invchebcoefs(), invertseries(), orthopoly()
 Macros related to special sequences and functions
   binom(), continfrac(), factorial(), i0(), i1()
 Other math related macros
   factors(), printfactors(), partitions()
%%%%
Macros related to matrices
  blockdiag    Construct block diagonal matrix bith given blocks
  kronecker    Compute the Kronecker product of two matrices
  matsqrt      Compute upper or lower triangular or symmetric matrix B
               such that B' %*% B = A, for given positive definite A
  moorepenrose Compute the Moore-Penrose inverse of a matrix
  qrdcomp      Compute the QR decomposition of a matrix

Macros for working with fully complex forms of complex matrices A and B
  cmatmultc    matrix product A %*% B, A %c% B or A %C% B
  ctranspose   A'
  cjtranspose  conj(A)'
  cdiag        diag(A)
  ctrace       trace(A)
  ceigen       eigenvalues and eigenvectors of Hermitian A
  csolve       inverse of A
  csubscr      simulated A[], A[i], A[i,], A[,j] or A[i,j]

Macros related to optimization
  dfp()        Optimization by Davidon-Fletcher-Powell method using
               golden section search
  bfs()        Optimization by Broyden-Fletcher-Shanno method using
               golden section search
  broyden()    Optimization by a method due to Broyden with no linear
               search
  neldermead() Optimization by Nelder Mead simplex method with optional
               quadratic polish
  levmar()     Macro implementing Levenberg-Marquart non-linear least
               squares based on a Fortran program of K. M. Brown
  _cgrad()     Macro used by levmar() to compute gradient, Jacobian,
               Hessian
  _lmout()     Macro used by levmar() to print out info on each
               iteration

Macros related to polynomials and series in powers of x
  chebcoefs()  Compute the coefficients of the expansion of a polynomial
               in Chebysev polynomials
  economize()  "Economize" a power series expansion on [-1, 1]
  invchebcoefs() Compute the coefficients powers of x in a finite
               series involving Chebysev polynomials
  invertseries() Macro to find coefficients of y^j of solution to
               y = sum(a[i]*x^i,i=1,n)
  orthopoly()  Compute standard orthogonal polynomials by recursion

Macros related to special sequences and functions
  binom()      Compute binomial coefficients
  continfrac() Compute continued fraction
               b0+a1/(b1+a2/(b2+a3/(b3+a4/... )))
  i0()         Compute modified Bessel function of the first kind I0(x)
  i1()         Compute modified Bessel function of the first kind I1(x)

Other math related macros
  factors()    Macro to compute prime factors of each element in a
               vector of positive integers
  printfactors() Macro to print output from factors
  partitions() Macro to compute partitions of integers n

====matsqrt()#matrices
%%%%
matsqrt(A [, symmetric:T or [lower:T]), square positive semi-definite
   matrix A with no MISSING values
%%%%
@@@@usage#Usage
B <- matsqrt(A) computes a matrix square root B of a positive
semi-definite REAL matrix A with no MISSING values.  B satisfies
B' %*% B = A.  B can also be computed by cholesky(A).

B <- matsqrt(A, lower:T) returns the lower triangular matrix square
root of A.

B <- matsqrt(A, symmetric:T) returns the symmetric matrix square
root of A.  B satisfies B %*% B = A

@@@@see_also#Cross reference
See also cholesky().
@@@@______

====minimizer()#minimize
%%%%
minimizer(x0, fun [,params:params] [, method:M] [, goldsteps:ngold] \
  [, maxit:maxiter] [,minit:miniter] [,criteria:vector(nsigx, \
  nsigfun,dgrad)] [printwhen:d1] [,recordwhen:d2]), REAL vector x0,
  macro fun(x,i [,params]), CHARACTER scalar M (one of "bfs", "dfp",
  "broyden", integers ngold > 0, maxiter >= 0, miniter > 0, nsigx,
  nsigfun, d1 >= 0, d2 >= 0, dgrad REAL scalar
%%%%
@@@@introduction#Introduction
Macro minimizer() uses a quasi-Newton variable metric algorithm to
minimize a function iteratively.  There is a choice of three methods,
Broyden-Fletcher-Shanno (method:"bfs"), Davidon-Fletcher-Powell
(method:"dfp") and Broyden's (method:"broyden").  For the first two a
golden mean line search is made at each step.  See Dahlquist and
Bjorck, Numerical methods, Prentice Hall, 1974, p. 441-444.

@@@@usage#Usage
result <- minimizer(x0, fun [, params] [,optional keywords]) computes
the minimum of a real function F(x1,x2,...,xk) starting iteration at x
= x0 = vector(x01,x02,...,x0k), a REAL vector with no MISSING elements.
The Broyden-Fletcher-Shanno updating is used by default.

result <- minimizer(x0, fun [, params] , method:M [,optional
keywords]), where M is one of "bfs", "dfp", or "broyden", does the
same, using the Broyden-Fletcher-Shanno, Davidon-Fletcher-Powell or
Broyden update methods.

Optional argument params is a variable, possibly a structure, with
additional constant information used to compute F(x) and its gradient
vector.

fun is a macro such that fun(x,0 [,params]) returns F(x[1],x[2],
...,x[k]) and fun(x,1 [,params]) returns a length k REAL gradient
vector (vector of partial derivatives of F(x) with respect to the
elements of x).  fun() should ignore its third argument if not needed.

fun(x,-1 [,params]) should carry out any initialization needed.  If
starting values or elements of param are not appropriate, fun(x, -1,
[params]) should return MISSING.  Otherwise, any non-MISSING value
should be returned.

fun() can use either a formula for derivatives, when one is known, or
compute them by numerical differentiation.

@@@@result#Result
result is structure(x:xmin, f:minVal, gradient:gradient, h:invhessian,
iterations:niter, status:N) where xmin is a REAL vector such that
minVal = F(xmin) is a local minimum, gradient is the length k gradient
vector at xmin, invhessian is a k by k approximation to the inverse of
the Hessian matrix (matrix of second order derivatives of F), niter is
the number of iterations taken and N is an integer indicating
convergence status; see below.

@@@@convergence_control#Convergence control
Optional keyword phrase criterion:vector(nsigx, nsigfun, dgrad) allows
control over how convergence is determined.  nsigx and nsigfun must be
integers and dgrad a small REAL scalar.  At least one of nsigx,
nsigfun and dgrad must be positive.  A value <= 0 is ignored.

Iteration ends when any of the following occurs
  1.  nsigx > 0 and all elements of x have relative change <=
  10^-nsigx.  For any x[i] with abs(x[i]) < .5, change is relative to
  .5.

  2.  nsigfun > 0 and relative change in F(x) <= 10^-nsigfun.  When
  abs(f(x)) < .5, change is relative to .5.

  3.  dgrad > 0 and ||gradient|| < dgrad

The default value for 'criterion' is vector(8,5,-1).

Keyword phrases minit:miniter and maxit:maxiter specify no check for
convergence is made until iteration miniter and no more than maxiter
iterations will be done.  miniter >= 0 (default 0) and maxiter > 0
(default 30) are integers.

@@@@convergence_status#Convergence status
When N = 0 (value of component 'status' of the result), no convergence
criterion was satisfied.

When N < 0, iteration was terminated because an illegal value was
encountered.

When N = 1 iteration ended by test using nsigx.

When N = 2 iteration ended by test using nsigfun.

When N = 3, iteration ended by test using dgrad.

After
  Cmd> result <- minimizer(x0, fun ..., maxit:n)

You can restart the iteration where it stopped by
  Cmd> result <- minimizer(result$x, fun, ..., h:result$h)

@@@@other_keywords#Other keywords
There other optional keyword phrases that can be arguments.

  h:invhessian      invhessian is a k by k REAL symmetric matrix
                    (default dmat(k,1)) used as starting approximation
                    to inverse Hessian matrix
  goldsteps:m       integer m >= 0 (default 5), the number of cycles
                    to be used in the golden mean linear search; with
                    m == 0 no linear search is done; ignored with
                    method:"broyden"
  printwhen:d1      Integer d1 >= 0.  When d1 > 0, current values of x,
                    F(x), and the gradient are printed on iterations
                    d1, 2*d1, 3*d1, ...
  recordwhen:d2     Integer d2 >= 0.  When d2 > 0, current values of
                    x, F(x) and the gradient on saved in components
                    'xvalx', 'funvals' and 'gradients' of side-
                    effect structure BFSRECORD, DJPRECORD or
                    BROYDNRECORD, depending on the method used.

@@@@see_also#Cross references
See also bfs(), dfp(), broyden(), and neldermead()
@@@@______

====moorepenrose()#matrices,generalized inverse
%%%%
moorepenrose(A), A a real matrix with no MISSING values
%%%%
@@@@usage#Usage
B <- moorepenrose(A), where A is a real matrix with no MISSING values
computes the Moore-Penrose inverse of a A.

B satisfies A %*% B %*% A = A and B %*% A %*% B = B.

When A is square and non-singular B = solve(A).

When A is m by n, B is n by m.  When m >= n and a is full rank b is
solve(a %c% a, a').

The result is computed from the singular value decomposition of A.

@@@@see_also#Cross references
See also solve(), svd().
@@@@______

====levmar()#nonlinear least squares,minimize
%%%%
levmar(b,x,y,f,param [,deriv:deriv,crit:crvec,active:active,\
  maxit:itmax,minit:itmin,print:T])
or
levmar(b,x,y,param ,resid:res [,deriv:deriv,crit:crvec,active:active,\
  maxit:itmax,minit:itmin,print:T])
  b        REAL vector of starting values for coefficients
  x        REAL variable.   Without 'resid:res'  a vector or matrix
           with nrows(x) = nrows(y); with 'resid:res' nrows(x) =
           norows(y) is not required
  y        REAL vector of data to be fit
  f        Macro:  fit <- f(b,x,param) returns a vector of fitted values
           with length nrows(x) = nrows(y); not allowed with 'resid:res'
  param    vector or structure of additional parameters for f or NULL
  res      Macro: res(b, x, y, param) computes a vector of residuals of
           length nrows(y).  'resid:res' is required when argument f is
           omitted and is not allowed when f is an argument.
  crvect   vector(numsig, nsigsq, delta), 3 criteria for convergence
  deriv    optional macro: deriv(b,x,y,param,j) computes derivative of
           f(b,x,param) or -res(b,x,y,param)  with respect to b[j]
  active   LOGICAL vector the same length as b
  itmax    integer >= 0, maximum number of iterations permitted (default
           = 30)
  itmin    0 <= minimum <= itmax = number of iterations performed (default
           = 1)
  print    if T, partial results are printed on each iteration

  Returned value is structure(coefs:b_hat,hessian:hes,jacobian:jac,
    gradient:g,rss:Rssmin,residuals:resids,nobs:n,iter:niter,
    iconv:convflag)
%%%%
@@@@introduction#Introduction
levmar() uses an analogue of the Levenberg-Marquardt and Gauss algorithms
to minimize a sum of squares.  Specifically, the criterion minimized is
sum(resids^2) where resids is computed as
    resids <- y - f(b, x, param)
or as
    resids <- res(b, x, y, param)
where f() or res() is a macro provided by the user

You can supply a macro to compute derivatives with respect to the elements
of b analytically or rely on difference-based numerical differentiation.

levmar() is intended for primarily use in higher level macros such as
nlreg() (non-linear regression) and arima() (estimation of ARIMA time
series models).

levmar() is based on a Fortran program for nonlinear least squares of Ken
Brown.  See below for references.

@@@@usage_and_arguments#Usage and arguments
levmar(b,x,y,f,param) uses an iterative algorithm to find a REAL vector
b_hat which minimizes Rss = sum((y - f(b_hat,x,param))^2), a sum of
squared residuals.

levmar(b,x,y,param, resid:res) does the same except the quantity minimized
is Rss = sum(res(b, x, y, param)^2).

In these usages, derivatives are computed by differences.

levmar(b,x,y,f,param, deriv:der) and levmar(b,x,y,param,resid:res,
deriv:der) do the same, except derivatives are computed by macro der()
instead of by differencing.

The arguments to levmar() are as follows
  b        REAL vector of starting values for b_hat
  x        REAL variable.  When f is an argument, x must be a vector or
           matrix with nrows(x) = nrows(y).  When 'resid:res' is an
           argument, nrows(x) is not required; if x is not used, it should
           be 0
  y        REAL vector of data to be fit
  param    vector, matrix or structure of additional fixed parameters or
           NULL
  f        macro called as fit <- f(b,x,param).  fit is a vector of
           containing nrows(x) = nrows(y) function values
  res      macro called as r <- res(b, x, y, param).  r is a vector
           of nrows(y) residuals
  der      macro called as derivs_j <- der(x,y,param,j).  derivs_j is
           vector containing the nrows(y) derivatives with respect to
           b[j] of the elements of f(b,x,param) or -res(b,x,y,param).

@@@@return_value#Return value
levmar() returns structure(coefs:b_hat,hessian:hes,jacobian:jac,
gradient:g,rss:Rssmin,residuals:resids,nobs:n,iter:niter,iconv:convflag)
where component values are as follows:

  b_hat     REAL vector which minimizes Rss
  hes       jac' %*% jac, an approximation to the Hessian matrix H,
            where H[i,j] = 2nd order partial derivative of Rss/2
            with respect to b_hat[i] and. b_hat[j]
  jac       the nrows(y) by nrows(b) Jacobian matrix; -jac[,j] = vector
            of partial derivatives of the elements of the residual
            vector with respect to b_hat[j].
  g         the nrows(b) REAL gradient vector with g[j] = partial
            derivative of Rss/2 with respect to b_hat[j].  g should be
            close to a vector of zeros.
  Rssmin    the minimized value of Rss
  resids    vector of residuals of length nrows(y)
  n         positive integer = nrows(y) = nrows(resids) = nrows(jac)
  niter     positive integer = the number of iterations
  conflag   Convergence status flag; 0 = not converged, 1 = met relative
            change in b_hat criterion, 2 = met relative change in Rss
            criterion, 3 = met norm of gradient vector criterion, 4 =
            failed to reduce Rss on a step of the iteration.

When any parameters are inactive as specified by keyword 'active' (see
below), jac is nrows(y) by p, hes is p by p and g has length p where p =
number of active parameters.

@@@@keyword_arguments#Keyword phrase arguments
There are several of keywords that can be used to control the iteration
and hold parameters at fixed values.

  Keyword Phrase   Value and Explanation
  crit:crvec       vector(numsig, nsigsq, delta), 3 criteria for
                   convergence (default = vector(8,5,-1)
                   numsig = desired number of significant digits in
                   the elements of b_hat (conflag = 1 when met)
                   nsigsq = desired number of significant digits in the
                   minimized Rss (conflag = 2 when met)
                   delta is a threshold for ||g||.  Iteration is stopped
                   when ||g|| <= delta (conflag = 3 when met)
                   A negative criterion is not used
  active:act       LOGICAL vector the same length as b.  b[j]
                   "participates" in the iteration only if act[j] is
                   True (default = rep(T,length(b))).  When act[j] is
                   False, b_hat[j] remains at the starting value
  maxit:itmax      Non-negative integer specifying the maximum number of
                   iterations (default = 30).  When itmax = 0, no iterations
                   are done and the quantities returned are computed at
                   the starting value b
  minit:itmin      Non-negative integer < itmax specifies the mininum
                   number of iterations
  print:T          When T, partial results are printed on each iteration

Iteration stops when (i) itmax is exceeded, (ii) any of the three
convergence criteria are satisfied, or (iii) when there has been no
reduction in Rss on an iteration after 10 halvings of the initial step
size.

@@@@convergence_criteria#Convergence criteria
There are 3 possible convergence criteria, at least one of which must be
enabled.  levmar() terminates iteration when at least itmin iterations
have been completed and any convergence criterion is satisfied or when
itmax iterations have been completed.

The criteria are specified by optional argument crit:vec, where vec is
vector(numsig [, nsigsq [, delta]]) with length <= 3.  The default values
for numsig, nsiqsq and delta are 8, 5 and -1, respectively.

A negative value for a criterion means it is not active.

  numsig   Desired number of significant digits in every active
           coefficient.  Specifically, the criterion is satisfied if,
           for every active coefficient b[j], the change d_j satisfies
           abs(d_j) < 10^-numsig*max(.5,abs(b_j)) where b_j is the
           updated value of b[j].  Component 'iconv' of the return value
           is 1 when satisfied.

  nsigsq   Desired number of significant digits in Rss = the residual
           sum of squares.  Specifically, the criterion is satisfied
           when abs(Rss_new - Rss_old) < 10^-nsigsq*max(.5,Rss_new).
           Component 'iconv' of the return value is 2 when satisfied.

  delta    Desired maximum norm ||g|| of the gradient vector g.
           Specifically, the criterion is satisfied when
           sqrt(sum(g^2)) <= delta.  Component 'iconv' of the return
           value is 3 when satisfied.

For numsig and nsigsq, the returned coefficients are the values updated
on that iteration.  For delta, the returned coefficients are the values
found on the previous iteration.

Criteria are checked in the order delta, numsig and nsigsq.

@@@@examples
Example:
Fit the function b1 + b2*b3^x to data from Snedecor and Cochran with
starting values b1 = b2 = 40 and b3 = 1.

  Cmd> x <- vector(0,1,2,3,4,5)

  Cmd> y <- vector(57.5, 45.7, 38.7, 35.3, 33.1, 32.2)

  Cmd> func <- macro("@b <- $1; @b[1]+@b[2]*@b[3]^($2)", dollars:T)

  Cmd> startVal <- vector(40,40,1) # starting values for iteration

  Cmd> stuff <- levmar(startVal,x,y,func,NULL); stuff
  component: coefs
  (1)      30.724      26.821     0.55184
  component: hessian
  (1,1)           6      2.1683      111.39
  (2,1)      2.1683      1.4367      30.242
  (3,1)      111.39      30.242      2675.8
  component: jacobian
  (1,1)           1           1           0
  (2,1)           1     0.55184      26.821
  (3,1)           1     0.30453      29.602
  (4,1)           1     0.16805      24.503
  (5,1)           1    0.092736      18.029
  (6,1)           1    0.051176      12.436
  component: gradient
  (1)  5.5896e-08 -8.4145e-09  2.4503e-05
  component: rss
  (1)    0.097248
  component: residuals
  (1)   -0.044919     0.17523    -0.19159    0.068869    -0.11115
  (6)     0.10356
  component: nobs
  (1)           6
  component: iter
  (1)           7
  component: iconv
  (1)           1

  Cmd> # compute MSE and approximate standard errors

  Cmd> mse <- stuff$rss/(stuff$nobs - length(stuff$coefs)); mse
  (1)     0.032416

  Cmd> sqrt(mse*diag(solve(stuff$hessian)))
  (1)      0.23099       0.2577     0.008448

  Cmd> # don't iterate over coefficient 1

  Cmd> startVal1 <- vector(30,40,1)

  Cmd> levmar(startVal1,x,y,func,NULL,active:vector(F,T,T))
  component: coefs
  (1)          30      27.418     0.57447
  component: hessian
  (1,1)      1.4907      34.492
  (2,1)      34.492      3136.2
  component: jacobian
  (1,1)           1           0
  (2,1)     0.57447      27.418
  (3,1)     0.33002      31.502
  (4,1)     0.18959      27.145
  (5,1)     0.10891      20.792
  (6,1)    0.062567      14.931
  component: gradient
  (1)  3.5552e-09  0.00067866
  component: rss
  (1)     0.38885
  component: residuals
  (1)     0.08214    -0.05082    -0.34842     0.10193     0.11385
  (6)     0.48454
  component: nobs
  (1)           6
  component: iter
  (1)           7
  component: iconv
  (1)           1

@@@@references#References
For information on the algorithm, see K M Brown and J E Dennis,
Derivative free analogues of the Levenberg-Marquardt and Gauss
algorithms for nonlinear least squares approximation, Numerische
Mathematik, Vol. 18, pp. 289-297 (1972), and K M Brown, Computer
oriented methods for fitting tabular data in the linear and nonlinear
least squares sense, Technical Report No. 72-13, University of
Minnesota Department of Computer and Information Sciences.
@@@@______

====neldermead()#minimize,direct search
%%%%
neldermead(fun, xstart, steps [,data] [, maxeval:maxeval, print:ip,\
  stopcrit:eps, nloop:nloop, quad:F or T, simpcrit:s]),
  fun a macro, xstart and steps REAL vectors with no MISSING elements,
  data a variable as required by fun(), integers maxeval > 0, nloop > 0
  and ip, REAL scalars eps > 0 and s > 0.  fun() specifies a
  function to be minimized and is called as fun(x) or fun(x,data),
  where x is a REAL vector the same length as xstart.
%%%%
@@@@introduction#Introduction
neldermead() is a macro implementing function minimization by direct
search using the simplex method.  The minimum found may be a local,
not a global, minimum.  For details, see Nelder & Mead, The Computer
Journal, January 1965.

@@@@usage#Usage
result <- neldermead(fun, xstart, steps) attempts to find a minimum of
the function F(x) computed by macro fun().

xstart, a REAL vector with no MISSING values, contains starting values
for x.

Macro fun() will be called as f <- fun(x, NULL) to evaluate
F(x), with argument 2 to be ignored.

steps is a REAL vector of non-negative numbers the same length as
xstart.  When steps[i] = 0, x[i] will remain fixed at xstart[i].  When
steps[i] > 0, it specifies an initial step size for building the
simplex.

result <- neldermead(fun, xstart, steps, data), where data is an
arbitrary variable, possibly a structure, does the same, except that
fun() will be called as f <- fun(x, data).  Argument data can contain
data or other constant information such as a fixed parameter.

result <- neldermead(fun, xstart, steps [,data], quad:T) does the same,
except a quadratic approximation is fitted to F(x) near the minimum
found by searching and then the minimum of that approximation is
found.

The value returned has the form structure(x:xmin, f:minval,
invhessian:v, neval:n, status:s).

@@@@returned_value#Returned value
The components of the structure returned are as follows:
  x:xmin        REAL vector contains the location of the minimum found
  f:minval      REAL scalar = F(xmin), the minimum attained
  invhessian:V  NULL or REAL square matrix describing the quadratic surface
                fitted with quad:T.
  neval:n       Integer > 0, the number of function evaluations required
  status:s      Integer >= 0 specifying the termination status. s = 0
                means successful termination at the apparent minimum; s
                = 1 means termination because of excessive function
                evaluations; s = 2 means v is not positive semidefinite
                (only with quad:T), indicating stationary point is not a
                minimum.

When F(x) = -log L(x), where x is a vector of parameters and L(x) is a
likelihood function, V is the inverse of the observed information matrix
and can be used as a variance-covariance matrix of the maximum
likelihood estimates.

When F(x) = sum((y - yhat(x))^2) is a sum of squared residuals, 2*MSE*V
is an estimate of the variance-covariance matrix of the least squares
estimates, where MSE = minval/edf, edf = error degrees of freedom.

@@@@keywords#Keyword phrase arguments
neldermead recognizes a number of additional keyword phrases.

  maxeval:m   Integer m > 0, the maximum number of function evaluations
              allowed; default is m = 1000
  crit:eps    REAL scalar eps > 0.  Search will terminate when SD <= eps,
              where SD = standard deviation of values of F(x) on a
              simplex.  The default is eps = 1e-5
  checkwhen:n Integer n > 0; the stopping rule is applied after every
              n function evaluations. The default is n = 10.
  simp:s      Small REAL scalar s > 0, a criterion for expanding the
              final simplex to overcome rounding errors before fitting
              the quadratic surface with quad:T.  The default is s =
              1e-8.  See below for a guideline.
  print:ip    Integer scalar ip controlling printing.  With ip < 0 (the
              default), nothing will be printed unless an error is
              found.  With ip = 0, the found minimum value and its
              location will be printed as well as warning messages.
              With ip > 0, the current value of x and F(x) will be
              printed every ip function evaluations.  When ip >= 0
              the returned structure is "invisible"; it can be
              assigned but won't be printed automatically.

@@@@advice_on_usage
                          Advice on usage
When the function minimized can be expected to be smooth in the vicinity
of the minimum, you are are strongly urged to use 'quad:T' to specify
the quadratic-surface fitting option.  This is the only satisfactory
way of testing that the minimum has been found.  When the fitted
quadratic surface is not positive definite (value of status = 2), it
probably means that the search terminated prematurely and you have not
found the minimum.

You should use simp:s, where s >= 10000 * E, where E = the rounding
error in calculating F(x). For example, the default simp:1e-8 is
appropriate when the rounding error in calculating F(x) may be of the
order of 1e-12.  If, say, the F(x) is computed by numerical integration
with accuracy on the order of 1e-05, simp:0.1 would be appropriate..

This advice is derived from the comments in the Fortran source on which
neldermead() was modeled.

@@@@history_and_references#History and references
neldermead() is based on a Fortran program with the following history
  Programmed by D.E.Shaw, CSIRO,  Division of Mathematics & Statistics
     P.O. Box 218,  Lindfield,  N.S.W. 2070
  With amendments by R.W.M.Wedderburn,Rothamsted Experimental Station,
     Harpenden,  Hertfordshire,  England
  Further amended by Alan Miller, CSIRO, Division of Mathematics &
     Statistics, Private Bag 10,  Clayton,  Vic. 3168

@@@@see_also#Cross references
For other macros to minimize functions, see minimizer(), bfs(), dfp()
and broyden().

@@@@______

====orthopoly()#orthogonal polynomials
%%%%
orthopoly(x, n, [,polycode [,parameters]]), x a REAL vector with no
  MISSING values, n >= 0 an integer, polycode one of p, j, g, t, u,l, h
  and d, parameters a REAL scalar or vector
%%%%
@@@@usage#Usage
orthopoly() is a macro which computes values of several of the standard
orthogonal polynomials.

P <- orthopoly(x, n), where x is a REAL vector with no MISSING values
and n > 0 is an integer, computes the matrix
  P = hconcat(P0(x), P1(x), ..., Pn(x))
where Pj(x) = vector(Pj(x[1], Pj(x[2]), ...) and Pj is the j-th
Legendre polynomial.  When x is a scalar, P is vector(P0(x),..., Pn(x)).

P <- orthopoly(x, n, polycode [,parameters]) does the same except
the type of polynomials is determined by polycode, an unquoted letter
which must be one of p, j, g, t, T, u, U, l, h and d.

parameters is required when polycode is j, G, or l, is optional when
polycode is d and should not be an argument otherwise.

@@@@available_polynomials#Available polynomials
The following table summarizes the options
  Polynomial type      polycode  parameters
  Legendre                p       none
  Jacobi                  j       vector(alpha,beta)
  Gegenbauer              g       alpha
  Chebyshev 1             t       none
  Shifted Chebyshev 1     T       none
  Chebyshev 2             u       none
  Shifted Chebyshev 2     U       none
  Laguerre                l       alpha
  Hermite                 h       none
  Discrete                d       vector w of with w[i] > 0; default
                                  is w = rep(1,length(x)).

@@@@discrete_polynomials#Discrete polynomials
The discrete polynomials are defined by both vector x and the vector of
weights.  They are orthogonal on the discrete set x[1], ..., x[m],
where m = length(x).  That is, they satisfy
  w[1]*Pj(x[1])*Pk(x[1]) + w[2]*Pj(x[2])*Pk(x[2]) + . . . +
     w[m]*Pj(x[m])*Pk(x[m]) = 0, j != k

@@@@recurrence_relations#Recurrence relations
All but the discrete polynomials are computed from the recurrence
relations in Table 22.7 of Handbook of Mathematical Functions by
Abramowitz and Stegun and satisfy the normalizations in Table 22.4.
The discrete polynomials are also computed by recursion and are
standardized so that sum(wj*poly(x[j])^2)/sum(wj) = 1
@@@@______

====partitions()#integers
%%%%
partitions(n [,all:T])
%%%%
@@@@usage#Usage
partitions(n), where n > 0 is an integer, computes a matrix with n
columns whose rows, possibly padded with 0, are the partitions of n.  A
partition of n is a non-increasing set of integers summing to n.

partitions(n,all:T) returns a structure R with n components, such that
R[j] is a matrix with j columns whose rows are the partitions of j
(padded with 0).

The number N of partitions of n grows fairly rapidly:
  n     N     n     N     n     N     n     N
  1     1     7    15    13   101    19   490
  2     2     8    22    14   135    20   627
  3     3     9    30    15   176    21   792
  4     5    10    42    16   231    22  1002
  5     7    11    56    17   297    23  1255
  6    11    12    77    18   385    24  1575

@@@@examples
Examples:
  Cmd> partitions(4) # all 5 partitions of 4; rows add to 4
  (1,1)          4          0          0          0
  (2,1)          3          1          0          0
  (3,1)          2          2          0          0
  (4,1)          2          1          1          0
  (5,1)          1          1          1          1

  Cmd> partitions(3,all:T) # partitions of 1, 2 and 3
  component: Parts_of_1
  (1,1)          1
  component: Parts_of_2
  (1,1)          2          0
  (2,1)          1          1
  component: Parts_of_3
  (1,1)          3          0          0
  (2,1)          2          1          0
  (3,1)          1          1          1
@@@@______

====printfactors()#prime factors
%%%%
printfactors(vector(n1 [n2, ...,]), integers n1 > 0, n2 > 0, ...
%%%%
@@@@usage#Usage
printfactors(vector(n), where n > 0 is an integer prints n and its prime
factors, if any.  It does not return any values.

printfactors(N), where N = vector(n1,n2,...), integers n1 > 0, n2 > 0,
..., prints each element of N together with its prime factors, if any.

@@@@example
Example:
  Cmd> printfactors(2^2^run(5) + 1)
  5 is prime
  17 is prime
  257 is prime
  65537 is prime
  4294967297 = 641*6700417

@@@@see_also#Cross references
See also factors() and primefactors().
@@@@______

====qrdcomp()#matrices,qr decomposition
%%%%
qrdcomp(x) or qrdcomp(qr(x)), where x is a REAL matrix with no MISSING
  values
%%%%
@@@@usage#Usage
qrdcomp(X) computes a QR decomposition without pivoting of M by N REAL
matrix X with no MISSING values.  It returns structure(q:Q,r:R), where
Q and R are REAL matrices satisfying X = Q %*% R.

When M >= N, Q is M by N orthonormal (Q' %*% Q = N by N identity
matrix) and R is N by N upper triangular.  When M < N, Q is M by M
orthonormal and R is M by N with the first M columns an upper
triangular matrix.

qrdcomp() also works with the output from function qr().  That is
qrdcomp(qr(X)) is the same as qrdcomp(X).

Note:  Computation of the QR decomposition without pivoting is not
reliable when X is not of full rank.

@@@@pivot_keyword#Keyword pivot
qrdcomp(X, pivot:T) or qrdcomp(qr(X, pivot:T)) a QR decomposition or X
with pivoting and return a structure(q:Q,r:R,pivot:J), where J is a
permutation of run(N), such that X[,J] = Q %*% R, where Q and R are as
before.  This usage is accurate even when X is not of full rank.

@@@@see_also#Cross reference
See also qr()
@@@@______

====_cgrad()*#
%%%%
stuff <- _cgrad(isw, b, x, y, params [,resids], [, active: active]\
  [, deriv:deriv]), integer isw = 1 or 2, REAL vector b, REAL vector or
  matrix x, REAL vector y, params, vector or structure, REAL vector
  resids, LOGICAL vector active.
%%%%
@@@@usage#Usage
_cgrad() computes the gradient and jacobian used in non-linear
estimation.

Usage is
  stuff <- _cgrad(isw, b, x, y, params [,resids], [, active: active]\
   [, deriv:deriv])
where
  isw     1 use (f(b+h)-f(b))/h with small h to compute derivatives
          2 use (f(b_h)-f(b-h)/(2*h) with larger h
  b       REAL vector of coefficients
  x       REAL vector or matrix
  y       REAL vector being fitted
  resids  optional REAL vector of residuals computed with b
  params  vector or structure of parameters used as argument to _resid()
          or deriv()
  active  LOGICAL vector the same length as b; a derivative with respect
          to b[j] is computed only if active[j] is True; default is
          rep(T,nrows(b))
  deriv   macro for computing derivatives called by
          deriv(b, x, y, param, j), where 1 <= j <= length(b) which
          needs to work for any j with active[j] True

The value returned is of the form structure(gradient:gradient,
jacobian:jacobian).

_cgrad() is used by levmar() and is unlikely to be useful otherwise.
@@@@______

====_lmout()*#
%%%%
_lmout(iter,theta, ssq, gradient, erl2, iconv)
%%%%
Macro used by levmar() to print partial results at each iteration

_E_O_F_#This should be the last line, an internal End Of File marker
