* The riemann theta function is defined by its mulit-dimensional Fourier series *
*
* θ(z|B) = ∑n∈Zg exp( ½(Bn,n)+(z,n) )
* ,
*
z∈Cg,B∈Cgxg.
* B is symmetric and has strictly negative definite real part.
*
* The exponential growth in z of the Riemann theta function can be factored
* out and the remaining oscillatory part can be approximated with a prescribed error
* tolerance. Thus we have:
*
*
* θ(z|B) = exp( factor(z|B) ) · thetaSum(z|B)
* ,
*
* Separating the exponential factor is crucial when dealing with theta functions. In practice,
* they often appear as ratios, in which case the exponential growth can be canceled.
* For example,
* let a,b,c,d ∈ Cg with
* a+b = c+d mod 2πiZg. Then
*
*
*
*
*
*
* θ(z+a|B)·θ(z+b|B)
* f(z) =
* θ(z+c|B)·θ(z+d|B)
* defines an Abelian function.
* The following code fragment is taken from {@link AbelianFunction},
* it shows how you evaluate f(z) using Theta:
*
*
*
*
*
*
* Theta theat = new Theta( B, tol );
*
* theta.theta( z.plus(a), factorOfZPlusA, thetaSumOfZPlusA );
* theta.theta( z.plus(b), factorOfZPlusB, thetaSumOfZPlusB );
* theta.theta( z.plus(c), factorOfZPlusC, thetaSumOfZPlusC );
* theta.theta( z.plus(d), factorOfZPlusD, thetaSumOfZPlusD );
*
* factor.assign( factorOfZPlusA );
* factor.assignPlus( factorOfZPlusB );
* factor.assignMinus( factorOfZPlusC );
* factor.assignMinus( factorOfZPlusD );
*
* fOfZ.assignExp( factor );
* fOfZ.assignTimes( thetaSumOfZPlusA );
* fOfZ.assignTimes( thetaSumOfZPlusB );
* fOfZ.assignDivide( thetaSumOfZPlusC );
* fOfZ.assignDivide( thetaSumOfZPlusD );
* Evaluating the theta function may easily lead to overflows or numerical
* instabilities, which should be avoid. ( The code above avoids the creation of
* temporarily used objects, which can be very beneficial for the performance. )
* The first line of the code creates a theta function for a period matrix
* B and a tolerance of tol for the thetaSum(z|B).
* This results in an absolute error of
*
*
*
*
*
* thetaSumZPlusA+thetaSumZPlusB
· exp(factor) · tol thetaSumZPlusC·thetaSumZPlusD
* for fOfZ. This error has be be controled if you want to be really shure what
* is going on. In practice it suffies to avoid the poles of f. This is true because
* the estimates, which lead to the approximation of the oscillating part of the Riemann
* theta function, are rather weak, and therefore the error for the actual values often falls low
* of the tolerance. This is also the reason why we have never experienced
* inaccuracies when evaluating first and second order derivatives of the Riemann
* theta function. Hence, we have not incorporated error estimates for these.
*
* To create a theta function object you must provide the period matrix and you may provide * an error tolerance. On top of that, you may switch off two important features, which generically * are used. *
* The first is about the fill factor error. The fill factor error is a heuristic to improve
* the accuracy of the error prediction. As mentioned earlier the estimates which lead to
* the approximation error are rather coarse, which is improved by this heuristic for
* the price of loosing the certainty of the 100% error. The smaller
* the fill factor, the bigger its influence on the approximation, which
* anyway is usually irrelevant for small genus (g<4).
*
* The second feature is Siegel`s reduction. {@link SiegelReduction} uses the modular
* transformation property of the Riemann theta function and can reduce the cost of
* computing theta functions drasticly. Sometimes you already know that your period
* matrix cannot be reduced. In that case you would switch of this feature to save the effort
* of checking for a possible reduction.
* The {@link ModularPropertySupport} is implemented without controling
* a nasty 8th-root of unity,
* which leads to an ambiguity. It is, however, usually neglectable because one only deals
* with ratios of theta functions. In case you are interested in theta functions themselves
* you still can take advantage of Siegel`s reduction: just evaluate theta at zero
* (or any other point) with and
* without Siegel`s reduction and determine that 8th-root of unity.
*
* Further reading: *
| - | B.Deconinck, M.Heil, A.I.Bobenko, M.van Hoeij, and M.Schmies. * Computing Riemann Theta Functions. Is to appear. |
| - | M.Heil. * Numerical Tools for the study of finite gap solutions of integrable systems. * PhD thesis, Technische Universität Berlin, 1995. |
| - | E.D.Belokolos, A.I.Bobenko, V.Z.Enolskii, A.R.Its, and V.B.Matveev. * Algebro-geometric appoach to nonlinear integrable probems. * Springer-Verlag Berlin, 1994. |
| - | C.L.Siegel. * Topics in complex function theory. Vol.III. * John Wiley & Sons, Inc., New York, 1989 |
periodMatrix.
* The error tolerance is by default e-7, the fill factor error
* and Siegel`s reduction is used.
* @param periodMatrix symmetric complex matrix with negative definite real part
*/
public Theta(final ComplexMatrix periodMatrix) {
this(periodMatrix, 1e-7, true, true);
}
/**
* Creates a Riemann theta function with prescribed periodMatrix
* and error tolerance tol.
* Both, fill factor error and Siegel`s reduction is used.
* @param periodMatrix symmetric complex matrix with negative definite real part
* @param tol positive number
*/
public Theta(final ComplexMatrix periodMatrix, final double tol) {
this(periodMatrix, tol, true, true);
}
/**
* Creates a Riemann theta function with prescribed periodMatrix
* and error tolerance tol.
* It further enables the user to switch off the Siegel`s Reduction.
* @param periodMatrix symmetric complex matrix with negative definite real part
* @param tol positive number
* @param performSiegelReduction controles whether or not Siegel`s reduction is performed
*/
public Theta(final ComplexMatrix periodMatrix, final double tol, final boolean performSiegelReduction) {
this(periodMatrix, tol, performSiegelReduction, true);
}
/**
* Creates a Riemann theta function
* with prescribed periodMatrix
* and error tolerance tol.
* It further enables the user to switch off the Siegel`s Reduction and the
* heuristic of the fill factor error.
* @param periodMatrix symmetric complex matrix with negative definite real part
* @param tol positive number
* @param performSiegelReduction controles whether or not Siegel`s reduction is performed
* @param useFillFactorError controles whether the fill factor or the 100% error is used.
*/
public Theta(final ComplexMatrix periodMatrix, final double tol, final boolean performSiegelReduction, final boolean useFillFactorError) {
this.performSiegelReduction = performSiegelReduction;
this.tol = tol;
setPeriodMatrix(periodMatrix);
}
/**
* Returns the fill factor.
*/
public final double getFillFactor() {
return fillFactor;
}
/**
* Returns whether the heuristic of the fill factor error is used or not.
*/
public final boolean getUseFillFactorError() {
return useFillFactorError;
}
/**
* Switches the heuristic fill factor error on or off.
* @param useFillFactorError controles whether the fill factor or the 100% error is used.
*/
public final void setUseFillFactorError(final boolean useFillFactorError) {
if (useFillFactorError == this.useFillFactorError)
return;
this.useFillFactorError = useFillFactorError;
computeRadius();
computeLatticePoints();
}
/**
* Returns the error tolerance for the oscillatory part of the theta function.
*/
public final double getTol() {
return tol;
}
/**
* Sets the error tolerance for the oscillatory part of the theta function.
* @param tol positive number
*/
public final void setTol(final double tol) {
if (this.tol == tol)
return;
this.tol = tol;
computeRadius();
computeLatticePoints();
}
/** Returns the dimension of the complex vector space on which the Riemann theta
functions operates.
*/
public final int getDim() {
return dim;
}
void setDim(final int v) {
if (v != dim) {
dim = v;
transform = new TransformPropertySupport(dim);
siegel = new SiegelReduction(dim);
modular = new ModularPropertySupport(dim);
gammaOfHalfDim = Gamma.gamma(dim / 2.0);
two2PowOfDim = Math.pow(2, dim);
sqrtOfPI2PowOfDim = Math.pow(Math.sqrt(Math.PI), dim);
}
}
/**
* Returns number of lattice points which are used for the uniform approximation
* of the oscillatory part of the Riemann theta function.
*/
public final int getNumOfLatticePoints() {
return numOfLatticePoints;
}
/**
* Returns period matrix on which the Riemann theta function operates.
*/
public final ComplexMatrix getPeriodMatrix() {
return new ComplexMatrix(periodMatrix);
}
public final ComplexMatrix getB() {
return new ComplexMatrix(B);
}
/**
* Sets periodMatrix on which the Riemann theta function operates.
* @param periodMatrix symmetric complex matrix with negative definite real part
*/
public final void setPeriodMatrix(final ComplexMatrix periodMatrix) {
if (this.periodMatrix != null && periodMatrix.equals(this.periodMatrix))
return;
if (!periodMatrix.isSquared())
throw new IllegalArgumentException("matrix is not squared");
if (!periodMatrix.isSymmetric())
throw new IllegalArgumentException("matrix is not symmetric");
setDim(periodMatrix.getNumRows());
if (this.periodMatrix == null)
this.periodMatrix = new ComplexMatrix(periodMatrix);
else
this.periodMatrix.assign(periodMatrix);
lastChangeOfPeriodMatrix = System.currentTimeMillis();
compute();
}
final void compute() {
if (performSiegelReduction && dim >= 1) {
siegel.setPeriodMatrix(periodMatrix);
modular.setModularTransformation(siegel.modular);
B = siegel.tPM;
T.assignTranspose(siegel.tL); // -re(B) = tL * T;
// the lattice: sqrt(1/2) T * Z^g
detOfLattice = T.get(0, 0);
for (int i = 1; i < dim; i++)
detOfLattice *= T.get(i, i);
detOfLattice /= Math.sqrt(two2PowOfDim);
lSLV = T.get(0, 0) / Math.sqrt(2);
} else {
modular.assignId();
B = periodMatrix;
B.getRe(T);
T.assignTimes(-1);
Cholesky.decompose(T.re, T.re);
T.assignTranspose();
// the lattice: sqrt(1/2) T * Z^g
detOfLattice = T.get(0, 0);
for (int i = 1; i < dim; i++)
detOfLattice *= T.get(i, i);
detOfLattice /= Math.sqrt(two2PowOfDim);
T.assignTranspose();
LLL.reduce(T.re, (int[][]) null);
lSLV = 0;
for (int i = 0; i < dim; i++)
lSLV += T.get(0, i) * T.get(0, i);
lSLV = Math.sqrt(lSLV /= 2);
}
lSLV2PowOfDim = Math.pow(lSLV, dim);
fillFactor = 2 * sqrtOfPI2PowOfDim * lSLV2PowOfDim / two2PowOfDim / dim / gammaOfHalfDim / detOfLattice;
B.getRe(reB);
B.getIm(imB);
modularIsId = modular.isId();
if (!modularIsId) {
modular.setPeriodMatrix(periodMatrix);
}
transform.setPeriodMatrix(B);
// the lattice: sqrt(1/2) T * Z^g
if (latticePointsforUniformApproximation == null) {
//latticePointsforUniformApproximation = new LatticePointsForUniformApproximation(reB.times(-0.5), true);
latticePointsforUniformApproximation = new LatticePointsForUniformApproximation(reB.times(-0.5));
} else {
latticePointsforUniformApproximation.B.assignTimes(reB, -0.5);
latticePointsforUniformApproximation.setB( latticePointsforUniformApproximation.B );
latticePointsforUniformApproximation.uptodate = false;
}
computeRadius();
computeLatticePoints();
}
/**
* Returns radius of ellipsoid which is used to determine the lattice
* points which are used for the uniform approximation
* of the oscillatory part of the Riemann theta function.
*/
public final double getRadius() {
return radius;
}
void computeLatticePoints() {
latticePointsforUniformApproximation.setRadius(radius);
latticePointsforUniformApproximation.update();
latticePoints = latticePointsforUniformApproximation.latticePoints;
numOfLatticePoints = latticePointsforUniformApproximation.numOfLatticePoints / 2 + 1;
// the thetaSum procedures expect the zero at the beginning
double [] tmp = latticePoints.re[numOfLatticePoints-1];
latticePoints.re[numOfLatticePoints-1] = latticePoints.re[0];;
latticePoints.re[0] = tmp;
if (expOfHalfBnn == null || expOfHalfBnn.size() < numOfLatticePoints) {
expOfHalfBnn = new ComplexVector(numOfLatticePoints);
}
final double[][] BRe = B.re;
final double[][] BIm = B.im;
for (int i = 0; i < numOfLatticePoints; i++) {
final double[] nRe = latticePoints.re[i];
double re = 0;
double im = 0;
for (int j = 0; j < dim; j++) {
final double[] rowBRe = BRe[j];
final double[] rowBIm = BIm[j];
re += rowBRe[j] * nRe[j] * nRe[j] / 2;
im += rowBIm[j] * nRe[j] * nRe[j] / 2;
for (int k = j + 1; k < dim; k++) {
re += rowBRe[k] * nRe[j] * nRe[k];
im += rowBIm[k] * nRe[j] * nRe[k];
}
}
exponent.assignExp(re, im);
expOfHalfBnn.set(i, exponent);
}
}
final void error0(double x, double[] value) {
double f = dim * two2PowOfDim / lSLV2PowOfDim / 2.0;
if (useFillFactorError) {
f *= fillFactor;
}
value[0] = x;
value[1] = f * Gamma.gamma(dim / 2.0, x, gammaOfHalfDim);
value[2] = -1 * f * Math.exp(Math.log(x) * (dim / 2.0 - 1) - x);
}
final class ErrorRadiusSolverFunction implements NewtonRaphson.RealFunctionWithDerivative, Serializable {
double error;
int order = 0;
double[] rootValues = new double[3];
public void eval(double x, double[] f, int offsetF, double[] df, int offsetDF) {
switch (order) {
case 1 :
break;
case 2 :
break;
default :
double factor = dim * two2PowOfDim / lSLV2PowOfDim / 2.0;
if (useFillFactorError) {
factor *= fillFactor;
}
f[offsetF] = factor * Gamma.gamma(dim / 2.0, x, gammaOfHalfDim);
df[offsetDF] = -1 * factor * Math.exp(Math.log(x) * (dim / 2.0 - 1) - x);
}
f[offsetF] -= error;
}
double getRadius(double anError, int anOrder) {
error = anError;
order = anOrder;
NewtonRaphson.search(this, dim / 2.0, rootValues);
return lSLV / 2 + Math.sqrt(rootValues[0]);
}
}
final void computeRadius() {
double harmonicThreshRadius = (Math.sqrt(dim + 2 * highestOrder + Math.sqrt(dim * dim + 8 * highestOrder)) + lSLV) / 2;
try {
radius = Math.max(ersf.getRadius(tol, 0), harmonicThreshRadius);
} catch (RuntimeException e) {
radius = harmonicThreshRadius;
}
}
void thetaSum(final ComplexVector Z, final Complex thetaSumZ) {
final double[] zRe = Z.re, zIm = Z.im;
final double[][] intLatticePointsRe = latticePoints.re;
final double[] expOfHalfBnnRe = expOfHalfBnn.re;
final double[] expOfHalfBnnIm = expOfHalfBnn.im;
thetaSumZ.assign(1, 0);
for (int i = 1; i < numOfLatticePoints; i++) {
final double[] nRe = intLatticePointsRe[i];
double nZRe = 0, nZIm = 0;
for (int j = 0; j < dim; j++) {
nZRe += zRe[j] * nRe[j];
nZIm += zIm[j] * nRe[j];
}
expOfNZ.assignExp(nZRe, nZIm);
invOfExpOfNZ.assignInvert(expOfNZ);
// compute thetaSumZ
term.assignPlus(expOfNZ, invOfExpOfNZ);
term.assignTimes(expOfHalfBnn.re[i], expOfHalfBnn.im[i]);
thetaSumZ.assignPlus(term);
}
}
final void dThetaSum(final ComplexVector Z, final ComplexVector X, final Complex thetaSumZ, final Complex thetaSumX) {
final double[] zRe = Z.re, zIm = Z.im;
final double[] uRe = X.re, uIm = X.im;
final double[][] intLatticePointsRe = latticePoints.re;
final double[] expOfHalfBnnRe = expOfHalfBnn.re, expOfHalfBnnIm = expOfHalfBnn.im;
thetaSumZ.assign(1, 0);
thetaSumX.assign(0, 0);
for (int i = 1; i < numOfLatticePoints; i++) {
final double[] nRe = intLatticePointsRe[i];
double nZRe = 0, nZIm = 0;
double nXRe = 0, nXIm = 0;
for (int j = 0; j < dim; j++) {
final double n = nRe[j];
nZRe += zRe[j] * n;
nZIm += zIm[j] * n;
nXRe += uRe[j] * n;
nXIm += uIm[j] * n;
}
expOfNZ.assignExp(nZRe, nZIm);
invOfExpOfNZ.assignInvert(expOfNZ);
// compute thetaSumZ
term.assignPlus(expOfNZ, invOfExpOfNZ);
term.assignTimes(expOfHalfBnn.re[i], expOfHalfBnn.im[i]);
thetaSumZ.assignPlus(term);
// compute thetaSumX
term.assignMinus(expOfNZ, invOfExpOfNZ);
term.assignTimes(expOfHalfBnn.re[i], expOfHalfBnn.im[i]);
term.assignTimes(nXRe, nXIm);
thetaSumX.assignPlus(term);
}
}
final void ddThetaSum(
final ComplexVector Z,
final ComplexVector X,
final ComplexVector Y,
final Complex thetaSumZ,
final Complex thetaSumX,
final Complex thetaSumY,
final Complex thetaSumXY) {
final double[] zRe = Z.re, zIm = Z.im;
final double[] uRe = X.re, uIm = X.im;
final double[] vRe = Y.re, vIm = Y.im;
final double[][] intLatticePointsRe = latticePoints.re;
final double[] expOfHalfBnnRe = expOfHalfBnn.re, expOfHalfBnnIm = expOfHalfBnn.im;
thetaSumZ.assign(1, 0);
thetaSumX.assign(0, 0);
thetaSumY.assign(0, 0);
thetaSumXY.assign(0, 0);
for (int i = 1; i < numOfLatticePoints; i++) {
final double[] nRe = intLatticePointsRe[i];
double nZRe = 0, nZIm = 0;
double nXRe = 0, nXIm = 0;
double nYRe = 0, nYIm = 0;
for (int j = 0; j < dim; j++) {
final double n = nRe[j];
nZRe += zRe[j] * n;
nZIm += zIm[j] * n;
nXRe += uRe[j] * n;
nXIm += uIm[j] * n;
nYRe += vRe[j] * n;
nYIm += vIm[j] * n;
}
expOfNZ.assignExp(nZRe, nZIm);
invOfExpOfNZ.assignInvert(expOfNZ);
// compute thetaSumZ
term.assignPlus(expOfNZ, invOfExpOfNZ);
term.assignTimes(expOfHalfBnn.re[i], expOfHalfBnn.im[i]);
thetaSumZ.assignPlus(term);
// compute thetaSumXY
term.assignTimes(nXRe, nXIm);
term.assignTimes(nYRe, nYIm);
thetaSumXY.assignPlus(term);
// compute thetaSumX and thetaSumY
term.assignMinus(expOfNZ, invOfExpOfNZ);
term.assignTimes(expOfHalfBnn.re[i], expOfHalfBnn.im[i]);
tmp.assign(term);
term.assignTimes(nXRe, nXIm);
thetaSumX.assignPlus(term);
term.assign(tmp);
term.assignTimes(nYRe, nYIm);
thetaSumY.assignPlus(term);
}
}
private final void getCombinedFactor(Complex f) {
f.assignPlus(modular.factor, transform.factor);
}
private final void getDerivativeOfCombinedFactor(ComplexVector X, Complex dFByDX) {
tmpVector.assignTimes(transform.M, modular.H);
tmpVector.assignPlus(modular.R);
tmpVector.assignMinus(modular.ATrZ);
tmpVector.assignMinus(modular.AZ);
ComplexVector.dotBilinear(tmpVector, X, dFByDX);
}
private final void getDerivativeOfCombinedFactor(ComplexVector X, ComplexVector Y, Complex ddFByDXDY) {
tmpVector.assignTimes(modular.A, X);
ComplexVector.dotBilinear(tmpVector, Y, ddFByDXDY);
tmpVector.assignTimes(modular.A, Y);
ComplexVector.dotBilinear(tmpVector, X, tmp);
ddFByDXDY.assignPlus(tmp);
ddFByDXDY.assignNeg();
}
private final void getDervativeOfCombinedZTransformation(ComplexVector X, ComplexVector dTByDX) {
dTByDX.assignTimes(modular.H, X);
}
private final void getContinousFactor(final Complex factor, final Complex continousCorrection) {
double continousFactor = -0.5 * transform.zRe.dot(transform.reBInvReZ);
continousCorrection.assign(transform.factor);
continousCorrection.re -= continousFactor;
factor.assignMinus(continousCorrection);
continousCorrection.assignExp();
}
/**
* Evaluates the Riemann theta function at Z.
*
*
* θ(z|B) = exp( factor ) · thetaSumZ
*
*
X direction at Z.
*
*
*
*
*
*
* θ(z|B) = exp( factor ) · thetaSumZ
* DXθ(z|B) = exp( factor ) · thetaSumX
* Evaluating the Riemann theta funciton and its first derivative simultanesly
* can be performed with almost no extra cost.
* @param Z argument vector
* @param X direction of derivative
* @param factor exponential part
* @param thetaSumZ oscillatory part
* @param thetaSumX oscillatory part of first derivative in X direction
*/
public final void dTheta(final ComplexVector Z, final ComplexVector X, final Complex factor, final Complex thetaSumZ, final Complex thetaSumX) {
if (modularIsId) {
transform.setZ(Z);
factor.assign(transform.factor);
dThetaSum(transform.transfromedZ, X, thetaSumZ, thetaSumX);
ComplexVector.dotBilinear(X, transform.M, MX);
tmp.assignTimes(MX, thetaSumZ);
thetaSumX.assignPlus(tmp);
} else {
modular.setZ(Z);
transform.setZ(modular.tZ);
getCombinedFactor(factor);
getDerivativeOfCombinedFactor(X, dFactorByDX);
getDervativeOfCombinedZTransformation(X, dTByDX);
dThetaSum(transform.transfromedZ, dTByDX, thetaSumZ, thetaSumX);
tmp.assignTimes(dFactorByDX, thetaSumZ);
thetaSumX.assignPlus(tmp);
}
getContinousFactor(factor, continousCorrection);
thetaSumZ.assignTimes(continousCorrection);
thetaSumX.assignTimes(continousCorrection);
}
/**
* Evaluates the Riemann theta function, its first derivatives
* in the X and Y direction and,
* its second derivative into the same direction at Z.
*
*
*
*
*
*
* θ(z|B) = exp( factor ) · thetaSumZ
* DXθ(z|B) = exp( factor ) · thetaSumX
* DYθ(z|B) = exp( factor ) · thetaSumY
* D²X,Yθ(z|B) = exp( factor ) · thetaSumXY
* Evaluating the Riemann theta funciton and its first derivative simultanesly
* can be performed with almost no extra cost.
* @param Z argument vector
* @param X direction of derivative
* @param Y direction of derivative
* @param factor exponential part
* @param thetaSumZ oscillatory part
* @param thetaSumX oscillatory part of first derivative in X direction
* @param thetaSumY oscillatory part of first derivative in Y direction
* @param thetaSumXY oscillatory part of second derivative
*/
public final void ddTheta(
final ComplexVector Z,
final ComplexVector X,
final ComplexVector Y,
final Complex factor,
final Complex thetaSumZ,
final Complex thetaSumX,
final Complex thetaSumY,
final Complex thetaSumXY) {
if (modularIsId) {
transform.setZ(Z);
factor.assign(transform.factor);
ddThetaSum(transform.transfromedZ, X, Y, thetaSumZ, thetaSumX, thetaSumY, thetaSumXY);
ComplexVector.dotBilinear(X, transform.M, MX);
ComplexVector.dotBilinear(Y, transform.M, MY);
// compute thetaSumXY
tmp.assignTimes(MX, thetaSumY);
thetaSumXY.assignPlus(tmp);
tmp.assignTimes(MY, thetaSumX);
thetaSumXY.assignPlus(tmp);
tmp.assignTimes(MY, thetaSumZ);
tmp.assignTimes(MX);
thetaSumXY.assignPlus(tmp);
// compute thetaSumX
tmp.assignTimes(MX, thetaSumZ);
thetaSumX.assignPlus(tmp);
// compute thetaSumY
tmp.assignTimes(MY, thetaSumZ);
thetaSumY.assignPlus(tmp);
} else {
modular.setZ(Z);
transform.setZ(modular.tZ);
getCombinedFactor(factor);
getDerivativeOfCombinedFactor(X, dFactorByDX);
getDerivativeOfCombinedFactor(Y, dFactorByDY);
getDerivativeOfCombinedFactor(X, Y, ddFactorByDXDY);
getDervativeOfCombinedZTransformation(X, dTByDX);
getDervativeOfCombinedZTransformation(Y, dTByDY);
ddThetaSum(transform.transfromedZ, dTByDX, dTByDY, thetaSumZ, thetaSumX, thetaSumY, thetaSumXY);
// compute thetaSumXY
tmp.assignTimes(dFactorByDX, dFactorByDY);
tmp.assignPlus(ddFactorByDXDY);
tmp.assignTimes(thetaSumZ);
thetaSumXY.assignPlus(tmp);
tmp.assignTimes(dFactorByDY, thetaSumX);
thetaSumXY.assignPlus(tmp);
tmp.assignTimes(dFactorByDX, thetaSumY);
thetaSumXY.assignPlus(tmp);
// compute thetaSumX
tmp.assignTimes(dFactorByDX, thetaSumZ);
thetaSumX.assignPlus(tmp);
// compute thetaSumY
tmp.assignTimes(dFactorByDY, thetaSumZ);
thetaSumY.assignPlus(tmp);
}
getContinousFactor(factor, continousCorrection);
thetaSumZ.assignTimes(continousCorrection);
thetaSumX.assignTimes(continousCorrection);
thetaSumY.assignTimes(continousCorrection);
thetaSumXY.assignTimes(continousCorrection);
}
}