177 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
178 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
198 assert(nmul_firstpart == 0);
199 Treal delta, beta, trD2;
204 if (nmul >= maxmul) {
207 "Purification reached maxmul"
208 " without convergence", maxmul);
210 if (tracediff[nmul] > 0) {
215 D = -ONE * X * X + TWO * D;
218 D.frob_thresh(frob_trunc);
219 idemerror[nmul] = Tmatrix::frob_diff(D, X);
221 tracediff[nmul] = D.trace() - nocc;
225 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
228 trD2 = (tracediff[nmul] + nocc -
229 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
230 (1 - 2 * polys[nmul - 1]);
233 while (ind > 0 && polys[ind] == polys[ind - 1]) {
234 delta = delta + frob_trunc;
239 }
while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
240 (tracediff[nmul - 1] + nocc)) /
241 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
242 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
243 (tracediff[nmul - 1] + nocc)) /
244 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
252 if (tracediff[nmul - 1] > 0) {
255 D = -ONE * X * X + TWO * D;
262 D.frob_thresh(frob_trunc);
263 idemerror[nmul] = Tmatrix::frob_diff(D, X);
265 tracediff[nmul] = D.trace() - nocc;
267 nmul_firstpart = nmul;
271 if (nmul + 1 >= maxmul) {
274 "Purification reached maxmul"
275 " without convergence", maxmul);
277 if (tracediff[nmul] > 0) {
279 idemerror[nmul] = Tmatrix::frob_diff(D, X);
283 tracediff[nmul] = D.trace() - nocc;
285 D = -ONE * X * X + TWO * D;
286 idemerror[nmul] = Tmatrix::frob_diff(D, X);
289 tracediff[nmul] = D.trace() - nocc;
293 X = -ONE * D * D + TWO * X;
294 idemerror[nmul] = Tmatrix::frob_diff(D, X);
297 tracediff[nmul] = X.trace() - nocc;
300 idemerror[nmul] = Tmatrix::frob_diff(D, X);
303 tracediff[nmul] = D.trace() - nocc;
305 D.frob_thresh(frob_trunc);
307 }
while (idemerror[nmul - 1] > frob_trunc);
325 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
326 if (idemerror[mul] < 1.0 / 4) {
328 tmp =
bisection(homofun, (Treal)0, (Treal)1, tol);
333 homo = tmp > homo ? tmp : homo;
336 return (lmin - lmax) * homo + lmax;
343 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
344 if (idemerror[mul] < 1.0 / 4) {
346 tmp =
bisection(lumofun, (Treal)0, (Treal)1, tol);
351 lumo = tmp < lumo ? tmp : lumo;
354 return (lmin - lmax) * lumo + lmax;