if (norm > 1)
s = (int) Math.ceil(Num.log2(norm));
else
s = 0;
Functions F = Functions.functions; // alias F
// B <-- B/2^s
double v = 1.0 / Math.pow(2.0, s);
if (v <= 0)
throw new IllegalArgumentException (" v <= 0");
B.assign (F.mult(v));
DoubleFactory2D fac = DoubleFactory2D.dense;
final DoubleMatrix2D B0 = fac.identity(n); // B^0 = I
final DoubleMatrix2D B2 = alge.mult(B, B); // B^2
final DoubleMatrix2D B4 = alge.mult(B2, B2); // B^4
DoubleMatrix2D T = B2.copy(); // T = work matrix
DoubleMatrix2D W = B4.copy(); // W = work matrix
W.assign (F.mult(cPade[9])); // W <-- W*cPade[9]
W.assign (T, F.plusMult(cPade[7])); // W <-- W + T*cPade[7]
DoubleMatrix2D U = alge.mult(B4, W); // U <-- B4*W
// T = B2.copy();
W = B4.copy();
W.assign (F.mult(cPade[5])); // W <-- W*cPade[5]
W.assign (T, F.plusMult(cPade[3])); // W <-- W + T*cPade[3]
W.assign (B0, F.plusMult(cPade[1])); // W <-- W + B0*cPade[1]
U.assign (W, F.plus); // U <-- U + W
U = alge.mult(B, U); // U <-- B*U
// T = B2.copy();
W = B4.copy();
W.assign (F.mult(cPade[8])); // W <-- W*cPade[8]
W.assign (T, F. plusMult(cPade[6])); // W <-- W + T*cPade[6]
DoubleMatrix2D V = alge.mult(B4, W); // V <-- B4*W
// T = B2.copy();
W = B4.copy();
W.assign (F.mult(cPade[4])); // W <-- W*cPade[4]
W.assign (T, F.plusMult(cPade[2])); // W <-- W + T*cPade[2]
W.assign (B0, F.plusMult(cPade[0])); // W <-- W + B0*cPade[0]
V.assign (W, F.plus); // V <-- V + W
W = V.copy();
W.assign(U, F.plus); // W = V + U, Padé numerator
T = V.copy();
T.assign(U, F.minus); // T = V - U, Padé denominator
// Compute Padé approximant for exponential = W / T
LUDecomposition lu = new LUDecomposition(T);
B = lu.solve(W);
if (false) {
// This overflows for large |mu|
// B <-- B^(2^s)
for(int i = 0; i < s; i++)
B = alge.mult(B, B);
/*
if (bal > 0) {
throw new UnsupportedOperationException (" balancing");
} */
v = Math.exp(mu);
B.assign (F.mult(v)); // B <-- B*e^mu
} else {
// equivalent to B^(2^s) * e^mu, but only if no balancing
double r = mu * v;
r = Math.exp(r);
B.assign (F.mult(r));
for (int i = 0; i < s; i++)
B = alge.mult(B, B);
}
return B;