CHROMA
remez_gmp.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Remez algorithm for finding nth roots
3  */
4 
6 
7 #define JMAX 1000 //Maximum number of iterations of Newton's approximation
8 
9 namespace Chroma
10 {
11 
12  // Constructor
14  {
15  START_CODE();
16 
17  alloc = false;
18 
19  END_CODE();
20  }
21 
22 
23  // Constructor
24  RemezGMP::RemezGMP(const Real& lower, const Real& upper, long precision)
25  {
26  START_CODE();
27 
28  prec = precision;
30 
31  alloc = false;
32  apstrt = lower;
33  apend = upper;
34  apwidt = apend - apstrt;
35  QDPIO::cout << __func__ << ": approximation bounds are ["
36  << (double)apstrt << " , " << (double)apend << "]"<< std::endl;
37 
38  n = 0;
39  d = 0;
40 
41  tolerance = 1e-15;
42 
43  // Only require the approximation spread to be less than tolerance
44  if (sizeof(Real)==sizeof(double)) {
45  tolerance = 1e-15;
46  } else if (sizeof(Real)==sizeof(float)) {
47  tolerance = 1e-8;
48  }
49 
50  END_CODE();
51  }
52 
53  // Free memory and reallocate as necessary
54  void RemezGMP::allocate(int num_degree, int den_degree)
55  {
56  START_CODE();
57 
58  // Note use of new and delete in memory allocation - cannot run on qcdsp
59  param.resize(num_degree+den_degree+1);
60  roots.resize(num_degree);
61  poles.resize(den_degree);
62  xx.resize(num_degree+den_degree+3);
63  mm.resize(num_degree+den_degree+2);
64 
65  alloc = true;
66 
67  END_CODE();
68  }
69 
70  // Reset the bounds of the approximation
71  void RemezGMP::setBounds(const Real& lower, const Real& upper)
72  {
73  START_CODE();
74 
75  apstrt = lower;
76  apend = upper;
77  apwidt = apend - apstrt;
78 
79  END_CODE();
80  }
81 
82  // Generate the rational approximation x^(pnum/pden)
83  Real RemezGMP::generateApprox(int degree, unsigned long pnum, unsigned long pden)
84  {
85  return generateApprox(degree, degree, pnum, pden);
86  }
87 
88  // Generate the rational approximation x^(pnum/pden)
89  Real RemezGMP::generateApprox(int num_degree, int den_degree,
90  unsigned long pnum, unsigned long pden)
91  {
92  START_CODE();
93 
94  // Reallocate arrays, since degree has changed
95  if (num_degree != n || den_degree != d) allocate(num_degree,den_degree);
96 
97  step.resize(num_degree+den_degree+2);
98 
99  power_num = pnum;
100  power_den = pden;
101  spread = 1.0e37;
102  iter = 0;
103 
104  n = num_degree;
105  d = den_degree;
106  neq = n + d + 1;
107 
108  initialGuess();
109  stpini(step);
110 
111  while (spread > tolerance)
112  {
113  //iterate until convergance
114 
115  if (iter++%100==0)
116  QDPIO::cout << __func__ << ": Iteration " << iter-1
117  << " spread " << spread << " delta " << delta << std::endl;
118 
119  equations();
120 
121  if (delta < tolerance)
122  {
123  QDPIO::cerr << __func__ << ":Delta too small, try increasing precision" << std::endl;;
124  QDP_abort(1);
125  }
126 
127  search(step);
128 
129  }
130 
131  int sign;
132  Real error = (Real)getErr(mm[0],sign);
133  QDPIO::cout << __func__ << " Converged at " << iter << " iterations, error = " << error << std::endl;
134 
135  // Once the approximation has been generated, calculate the roots
136  if (!root())
137  {
138  QDPIO::cerr << __func__ << ": Root finding failed" << std::endl;
139  QDP_abort(1);
140  }
141 
142  END_CODE();
143 
144  // Return the maximum error in the approximation
145  return error;
146  }
147 
148  // Return the partial fraction expansion of the approximation x^(pnum/pden)
150  {
151  START_CODE();
152 
153  if (n!=d)
154  {
155  QDPIO::cerr << __func__ << ": Cannot handle case: Numerator degree neq Denominator degree" << std::endl;
156  QDP_abort(1);
157  }
158 
159  if (!alloc)
160  {
161  QDPIO::cerr << __func__ << ": Approximation not yet generated" << std::endl;
162  QDP_abort(1);
163  }
164 
165  RemezCoeff_t coeff;
166  coeff.res.resize(n);
167  coeff.pole.resize(d);
168 
169  multi1d<bigfloat> r(n);
170  multi1d<bigfloat> p(d);
171 
172  for (int i=0; i<n; i++) r[i] = roots[i];
173  for (int i=0; i<d; i++) p[i] = poles[i];
174 
175  // Perform a partial fraction expansion
176  pfe(r, p, norm);
177 
178  // Convert to Real and return
179  coeff.norm = (double)norm;
180  for (int i=0; i<n; i++) coeff.res[i] = (double)r[i];
181  for (int i=0; i<d; i++) coeff.pole[i] = (double)p[i];
182 
183  END_CODE();
184 
185  return coeff;
186  }
187 
188  // Return the partial fraction expansion of the approximation x^(-pnum/pden)
190  {
191  START_CODE();
192 
193  if (!alloc)
194  {
195  QDPIO::cerr << __func__ << ": Approximation not yet generated" << std::endl;
196  QDP_abort(1);
197  }
198 
199  RemezCoeff_t coeff;
200  coeff.res.resize(n);
201  coeff.pole.resize(d);
202 
203  multi1d<bigfloat> r(d);
204  multi1d<bigfloat> p(n);
205 
206  // Want the inverse function
207  for (int i=0; i<n; i++) {
208  r[i] = poles[i];
209  p[i] = roots[i];
210  }
211 
212  // Perform a partial fraction expansion
213  pfe(r, p, (bigfloat)1l/norm);
214 
215  // Convert to Real and return
216  coeff.norm = (double)((bigfloat)1l/(norm));
217  for (int i=0; i<n; i++) {
218  coeff.res[i] = (double)r[i];
219  coeff.pole[i] = (double)p[i];
220  }
221 
222  END_CODE();
223 
224  return coeff;
225  }
226 
227  // Initial values of maximal and minimal errors
229  {
230  START_CODE();
231 
232  // Supply initial guesses for solution points
233  long ncheb = neq; // Degree of Chebyshev error estimate
234  bigfloat a, r;
235 
236  // Find ncheb+1 extrema of Chebyshev polynomial
237 
238  a = ncheb;
239  mm[0] = apstrt;
240  for (long i = 1; i < ncheb; i++) {
241  r = 0.5 * (1 - cos((M_PI * i)/(double) a));
242  //r *= sqrt_bf(r);
243  r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
244  mm[i] = apstrt + r * apwidt;
245  }
246  mm[ncheb] = apend;
247 
248  a = 2.0 * ncheb;
249  for (long i = 0; i <= ncheb; i++) {
250  r = 0.5 * (1 - cos(M_PI * (2*i+1)/(double) a));
251  //r *= sqrt_bf(r); // Squeeze to low end of interval
252  r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
253  xx[i] = apstrt + r * apwidt;
254  }
255 
256  END_CODE();
257  }
258 
259  // Initialise step sizes
260  void RemezGMP::stpini(multi1d<bigfloat>& step)
261  {
262  START_CODE();
263 
264  if (step.size() == 0)
265  QDP_error_exit("%s: step not allocated", __func__);
266 
267  xx[neq+1] = apend;
268  delta = 0.25;
269  step[0] = xx[0] - apstrt;
270  for (int i = 1; i < neq; i++) step[i] = xx[i] - xx[i-1];
271  step[neq] = step[neq-1];
272 
273  END_CODE();
274  }
275 
276  // Search for error maxima and minima
277  void RemezGMP::search(multi1d<bigfloat>& step)
278  {
279  START_CODE();
280 
281  if (step.size() == 0)
282  QDP_error_exit("%s: step not allocated", __func__);
283 
284  bigfloat a, q, xm, ym, xn, yn, xx0, xx1;
285  int i, j, meq, emsign, ensign, steps;
286 
287  meq = neq + 1;
288  multi1d<bigfloat> yy(meq);
289 
290  bigfloat eclose = 1.0e30;
291  bigfloat farther = 0l;
292 
293  j = 1;
294  xx0 = apstrt;
295 
296  for (i = 0; i < meq; i++) {
297  steps = 0;
298  xx1 = xx[i]; // Next zero
299  if (i == meq-1) xx1 = apend;
300  xm = mm[i];
301  ym = getErr(xm,emsign);
302  q = step[i];
303  xn = xm + q;
304  if (xn < xx0 || xn >= xx1) { // Cannot skip over adjacent boundaries
305  q = -q;
306  xn = xm;
307  yn = ym;
308  ensign = emsign;
309  } else {
310  yn = getErr(xn,ensign);
311  if (yn < ym) {
312  q = -q;
313  xn = xm;
314  yn = ym;
315  ensign = emsign;
316  }
317  }
318 
319  while(yn >= ym) { // March until error becomes smaller.
320  if (++steps > 10) break;
321  ym = yn;
322  xm = xn;
323  emsign = ensign;
324  a = xm + q;
325  if (a == xm || a <= xx0 || a >= xx1) break;// Must not skip over the zeros either side.
326  xn = a;
327  yn = getErr(xn,ensign);
328  }
329 
330  mm[i] = xm; // Position of maximum
331  yy[i] = ym; // Value of maximum
332 
333  if (eclose > ym) eclose = ym;
334  if (farther < ym) farther = ym;
335 
336  xx0 = xx1; // Walk to next zero.
337  } // end of search loop
338 
339  q = (farther - eclose); // Decrease step size if error spread increased
340  if (eclose != 0.0) q /= eclose; // Relative error spread
341  if (q >= spread) delta *= 0.5; // Spread is increasing; decrease step size
342  spread = q;
343 
344  for (i = 0; i < neq; i++) {
345  q = yy[i+1];
346  if (q != 0.0) q = yy[i] / q - (bigfloat)1l;
347  else q = 0.0625;
348  if (q > (bigfloat)0.25) q = 0.25;
349  q *= mm[i+1] - mm[i];
350  step[i] = q * delta;
351  }
352  step[neq] = step[neq-1];
353 
354  for (i = 0; i < neq; i++) { // Insert new locations for the zeros.
355  xm = xx[i] - step[i];
356  if (xm <= apstrt) continue;
357  if (xm >= apend) continue;
358  if (xm <= mm[i]) xm = (bigfloat)0.5 * (mm[i] + xx[i]);
359  if (xm >= mm[i+1]) xm = (bigfloat)0.5 * (mm[i+1] + xx[i]);
360  xx[i] = xm;
361  }
362 
363  END_CODE();
364  }
365 
366  // Solve the equations
368  {
369  START_CODE();
370 
371  bigfloat x, y, z;
372  int i, j, ip;
373  bigfloat *aa;
374 
375  multi1d<bigfloat> AA(neq*neq);
376  multi1d<bigfloat> BB(neq);
377 
378  for (i = 0; i < neq; i++) { // set up the equations for solution by simq()
379  ip = neq * i; // offset to 1st element of this row of matrix
380  x = xx[i]; // the guess for this row
381  y = func(x); // right-hand-side std::vector
382 
383  z = (bigfloat)1l;
384  aa = &AA[ip];
385  for (j = 0; j <= n; j++) {
386  *aa++ = z;
387  z *= x;
388  }
389 
390  z = (bigfloat)1l;
391  for (j = 0; j < d; j++) {
392  *aa++ = -y * z;
393  z *= x;
394  }
395  BB[i] = y * z; // Right hand side std::vector
396  }
397 
398  // Solve the simultaneous linear equations.
399  if (simq(AA, BB, param, neq))
400  {
401  QDPIO::cerr << __func__ << ": simq failed" << std::endl;
402  QDP_abort(1);
403  }
404 
405  END_CODE();
406  }
407 
408  // Evaluate the rational form P(x)/Q(x) using coefficients
409  // from the solution std::vector param
411  {
412  START_CODE();
413 
414  bigfloat yn, yd;
415  int i;
416 
417  // Work backwards toward the constant term.
418  yn = param[n]; // Highest order numerator coefficient
419  for (i = n-1; i >= 0; i--) yn = x * yn + param[i];
420  yd = x + param[n+d]; // Highest degree coefficient = 1.0
421  for (i = n+d-1; i > n; i--) yd = x * yd + param[i];
422 
423  END_CODE();
424 
425  return(yn/yd);
426  }
427 
428  // Compute size and sign of the approximation error at x
429  bigfloat RemezGMP::getErr(const bigfloat& x, int& sign)
430  {
431  bigfloat e, f;
432 
433  f = func(x);
434  e = approx(x) - f;
435  if (f != 0) e /= f;
436  if (e < (bigfloat)0.0) {
437  sign = -1;
438  e = -e;
439  }
440  else
441  sign = 1;
442 
443  return(e);
444  }
445 
446  // Calculate function required for the approximation
448  {
449  START_CODE();
450 
451  bigfloat y,dy,f=1l,df;
452 
453  // initial guess to accelerate convergance
454  y = (bigfloat)pow((double)x,(double)((bigfloat)1l/(bigfloat)power_den));
455  while (abs_bf(f)>(bigfloat)1l/pow_bf((bigfloat)10,prec)) { // approx good to 10^(-prec)
456  f = pow_bf(y,power_den) - x;
457  df = (bigfloat)power_den*pow_bf(y,power_den-1);// need power_den-1 because of diff
458  dy = f/df;
459  y -= dy;
460  }
461 
462  END_CODE();
463 
464  return pow_bf(y,power_num);
465  }
466 
467  // Solve the system AX=B
468  int RemezGMP::simq(multi1d<bigfloat>& A, multi1d<bigfloat>& B, multi1d<bigfloat>& X, int n)
469  {
470  START_CODE();
471 
472  if (A.size() == 0 || B.size() == 0)
473  QDP_error_exit("%s: A or B not allocated", __func__);
474 
475  int i, j, ij, ip, ipj, ipk, ipn;
476  int idxpiv, iback;
477  int k, kp, kp1, kpk, kpn;
478  int nip, nkp, nm1;
479  bigfloat em, q, rownrm, big, size, pivot, sum;
480  bigfloat *aa;
481 
482  multi1d<int> IPS(neq); // simq() work std::vector
483 
484  nm1 = n - 1;
485  // Initialize IPS and X
486 
487  ij = 0;
488  for (i = 0; i < n; i++) {
489  IPS[i] = i;
490  rownrm = 0.0;
491  for(j = 0; j < n; j++) {
492  q = abs_bf(A[ij]);
493  if(rownrm < q) rownrm = q;
494  ++ij;
495  }
496  if (rownrm == (bigfloat)0l)
497  {
498  QDPIO::cout << "simq rownrm=0\n" << std::endl;
499  return(1);
500  }
501  X[i] = (bigfloat)1.0 / rownrm;
502  }
503 
504  for (k = 0; k < nm1; k++) {
505  big = 0.0;
506  idxpiv = 0;
507 
508  for (i = k; i < n; i++) {
509  ip = IPS[i];
510  ipk = n*ip + k;
511  size = abs_bf(A[ipk]) * X[ip];
512  if (size > big) {
513  big = size;
514  idxpiv = i;
515  }
516  }
517 
518  if (big == (bigfloat)0l)
519  {
520  QDPIO::cout << "simq big=0" << std::endl;
521  return(2);
522  }
523  if (idxpiv != k) {
524  j = IPS[k];
525  IPS[k] = IPS[idxpiv];
526  IPS[idxpiv] = j;
527  }
528  kp = IPS[k];
529  kpk = n*kp + k;
530  pivot = A[kpk];
531  kp1 = k+1;
532  for (i = kp1; i < n; i++) {
533  ip = IPS[i];
534  ipk = n*ip + k;
535  em = -A[ipk] / pivot;
536  A[ipk] = -em;
537  nip = n*ip;
538  nkp = n*kp;
539  aa = &A[nkp+kp1];
540  for (j = kp1; j < n; j++) {
541  ipj = nip + j;
542  A[ipj] = A[ipj] + em * *aa++;
543  }
544  }
545  }
546  kpn = n * IPS[n-1] + n - 1; // last element of IPS[n] th row
547  if (A[kpn] == (bigfloat)0l)
548  {
549  QDPIO::cout << "simq A[kpn]=0" << std::endl;
550  return(3);
551  }
552 
553 
554  ip = IPS[0];
555  X[0] = B[ip];
556  for (i = 1; i < n; i++) {
557  ip = IPS[i];
558  ipj = n * ip;
559  sum = 0.0;
560  for (j = 0; j < i; j++) {
561  sum += A[ipj] * X[j];
562  ++ipj;
563  }
564  X[i] = B[ip] - sum;
565  }
566 
567  ipn = n * IPS[n-1] + n - 1;
568  X[n-1] = X[n-1] / A[ipn];
569 
570  for (iback = 1; iback < n; iback++) {
571  //i goes (n-1),...,1
572  i = nm1 - iback;
573  ip = IPS[i];
574  nip = n*ip;
575  sum = 0.0;
576  aa = &A[nip+i+1];
577  for (j= i + 1; j < n; j++)
578  sum += *aa++ * X[j];
579  X[i] = (X[i] - sum) / A[nip+i];
580  }
581 
582  END_CODE();
583 
584  return(0);
585  }
586 
587  // Calculate the roots of the approximation
589  {
590  START_CODE();
591 
592  long i,j;
593  bigfloat x,dx=0.05;
594  bigfloat upper=1, lower=-100000;
595  bigfloat tol = 1e-20;
596 
597  multi1d<bigfloat> poly(neq+1);
598 
599  // First find the numerator roots
600  for (i=0; i<=n; i++) poly[i] = param[i];
601  for (i=n-1; i>=0; i--) {
602  roots[i] = rtnewt(poly,i+1,lower,upper,tol);
603  if (roots[i] == 0.0)
604  {
605  QDPIO::cerr << __func__ << ": Failure to converge on root" << std::endl;
606  return 0;
607  }
608  poly[0] = -poly[0]/roots[i];
609  for (j=1; j<=i; j++) poly[j] = (poly[j-1] - poly[j])/roots[i];
610  }
611 
612  // Now find the denominator roots
613  poly[d] = 1l;
614  for (i=0; i<d; i++) poly[i] = param[n+1+i];
615  for (i=d-1; i>=0; i--) {
616  poles[i]=rtnewt(poly,i+1,lower,upper,tol);
617  if (poles[i] == 0.0) {
618  QDPIO::cerr << __func__ << ": Failure to converge on root" << std::endl;
619  return 0;
620  }
621  poly[0] = -poly[0]/poles[i];
622  for (j=1; j<=i; j++) poly[j] = (poly[j-1] - poly[j])/poles[i];
623  }
624 
625  norm = param[n];
626  QDPIO::cout << __func__ << ": Normalisation constant is " << (double)norm << std::endl;
627  for (i=0; i<n; i++)
628  QDPIO::cout << "root[" << i << "] = " << (double)roots[i] << std::endl;
629  for (i=0; i<d; i++)
630  QDPIO::cout << "pole[" << i << "] = " << (double)poles[i] << std::endl;
631 
632  END_CODE();
633 
634  return 1;
635  }
636 
637  // Evaluate the polynomial
638  bigfloat RemezGMP::polyEval(const bigfloat& x, const multi1d<bigfloat>& poly, long size)
639  {
640  bigfloat f = poly[size];
641  for (int i=size-1; i>=0; i--) f = f*x + poly[i];
642  return f;
643  }
644 
645  // Evaluate the differential of the polynomial
646  bigfloat RemezGMP::polyDiff(const bigfloat& x, const multi1d<bigfloat>& poly, long size)
647  {
648  bigfloat df = (bigfloat)size*poly[size];
649  for (int i=size-1; i>0; i--)
650  df = df*x + (bigfloat)i*poly[i];
651  return df;
652  }
653 
654  // Newton's method to calculate roots
655  bigfloat RemezGMP::rtnewt(const multi1d<bigfloat>& poly, long i,
656  const bigfloat& x1, const bigfloat& x2, const bigfloat& xacc)
657  {
658  START_CODE();
659 
660  int j;
661  bigfloat df, dx, f, rtn;
662 
663  rtn=(bigfloat)0.5*(x1+x2);
664  for (j=1; j<=JMAX;j++)
665  {
666  f = polyEval(rtn, poly, i);
667  df = polyDiff(rtn, poly, i);
668  dx = f/df;
669  rtn -= dx;
670  if ((x1-rtn)*(rtn-x2) < (bigfloat)0.0)
671  QDPIO::cerr << __func__ << ": Jumped out of brackets in rtnewt" << std::endl;
672  if (abs_bf(dx) < xacc) return rtn;
673  }
674  QDPIO::cerr << __func__ << ": Maximum number of iterations exceeded in rtnewt" << std::endl;
675 
676  END_CODE();
677  return 0.0;
678  }
679 
680  // Evaluate the partial fraction expansion of the rational function
681  // with res roots and poles poles. Result is overwritten on input
682  // arrays.
683  void RemezGMP::pfe(multi1d<bigfloat>& res, multi1d<bigfloat>& poles, const bigfloat& norm)
684  {
685  START_CODE();
686 
687 // QDPIO::cout << __func__ << " : enter" << std::endl;
688 
689  if (res.size() == 0 || poles.size() == 0)
690  QDP_error_exit("%s: res or poles not allocated", __func__);
691 
692  int i,j,small;
693  bigfloat temp;
694  multi1d<bigfloat> numerator(n);
695  multi1d<bigfloat> denominator(d);
696 
697  // Construct the polynomials explicitly
698  for (i=1; i<n; i++) {
699  numerator[i] = 0l;
700  denominator[i] = 0l;
701  }
702  numerator[0]=1l;
703  denominator[0]=1l;
704 
705  for (j=0; j<n; j++) {
706  for (i=n-1; i>=0; i--) {
707  numerator[i] *= -res[j];
708  denominator[i] *= -poles[j];
709  if (i>0) {
710  numerator[i] += numerator[i-1];
711  denominator[i] += denominator[i-1];
712  }
713  }
714  }
715 
716  // Convert to proper fraction form.
717  // Fraction is now in the form 1 + n/d, where O(n)+1=O(d)
718  for (i=0; i<n; i++) numerator[i] -= denominator[i];
719 
720  // Find the residues of the partial fraction expansion and absorb the
721  // coefficients.
722  for (i=0; i<n; i++) {
723  res[i] = 0l;
724  for (j=n-1; j>=0; j--) {
725  res[i] = poles[i]*res[i]+numerator[j];
726  }
727  for (j=n-1; j>=0; j--) {
728  if (i!=j) res[i] /= poles[i]-poles[j];
729  }
730  res[i] *= norm;
731  }
732 
733  // res now holds the residues
734  j = 0;
735  for (i=0; i<n; i++) poles[i] = -poles[i];
736 
737  // Move the ordering of the poles from smallest to largest
738  for (j=0; j<n; j++)
739  {
740  small = j;
741  for (i=j+1; i<n; i++) {
742  if (poles[i] < poles[small]) small = i;
743  }
744  if (small != j) {
745  temp = poles[small];
746  poles[small] = poles[j];
747  poles[j] = temp;
748  temp = res[small];
749  res[small] = res[j];
750  res[j] = temp;
751  }
752  QDPIO::cout << __func__ << ": Residue = " << (double)res[j] << " Pole = " << (double)poles[j] << std::endl;
753  }
754 
755 // QDPIO::cout << __func__ << " : exit" << std::endl;
756 
757  END_CODE();
758  }
759 
760 
761  // Given a partial fraction expansion, evaluate it at x
762  Real RemezGMP::evalPFE(const Real& x, const RemezCoeff_t& coeff)
763  {
764  if ((coeff.res.size() != coeff.pole.size()) || coeff.res.size() == 0)
765  {
766  QDPIO::cerr << __func__ << ": invalid res and pole" << std::endl;
767  QDP_abort(1);
768  }
769 
770  Real f = coeff.norm;
771  for (int i=0; i < coeff.res.size(); ++i)
772  f += coeff.res[i] / (x + coeff.pole[i]);
773 
774  return f;
775  }
776 
777 } // end namespace Chroma
778 
bigfloat polyDiff(const bigfloat &x, const multi1d< bigfloat > &poly, long size)
Evaluate the differential of the polynomial.
Definition: remez_gmp.cc:646
bigfloat delta
Definition: remez_gmp.h:77
bigfloat getErr(const bigfloat &x, int &sign)
Compute size and sign of the approximation error at x.
Definition: remez_gmp.cc:429
bool alloc
Flag to determine whether the arrays have been allocated.
Definition: remez_gmp.h:72
int simq(multi1d< bigfloat > &A, multi1d< bigfloat > &B, multi1d< bigfloat > &X, int n)
Solve the system AX=B.
Definition: remez_gmp.cc:468
bigfloat func(const bigfloat &x)
Calculate function required for the approximation.
Definition: remez_gmp.cc:447
void allocate(int num_degree, int den_degree)
Free memory and reallocate as necessary.
Definition: remez_gmp.cc:54
int neq
The number of equations we must solve at each iteration (n+d+1)
Definition: remez_gmp.h:80
bigfloat approx(const bigfloat &x)
Definition: remez_gmp.cc:410
void initialGuess()
Initial values of maximal and minmal errors.
Definition: remez_gmp.cc:228
multi1d< bigfloat > xx
Definition: remez_gmp.h:76
bigfloat rtnewt(const multi1d< bigfloat > &poly, long i, const bigfloat &x1, const bigfloat &x2, const bigfloat &xacc)
Newton's method to calculate roots.
Definition: remez_gmp.cc:655
RemezCoeff_t getIPFE()
Return the partial fraction expansion of the approximation x^(-pnum/pden)
Definition: remez_gmp.cc:189
multi1d< bigfloat > mm
Definition: remez_gmp.h:76
unsigned long power_den
Definition: remez_gmp.h:69
Real evalPFE(const Real &x, const RemezCoeff_t &coeff)
Given a partial fraction expansion, evaluate it at x.
Definition: remez_gmp.cc:762
bigfloat tolerance
Definition: remez_gmp.h:77
bigfloat apwidt
Definition: remez_gmp.h:65
multi1d< bigfloat > roots
Definition: remez_gmp.h:58
bigfloat norm
Definition: remez_gmp.h:59
int n
The numerator and denominator degree (n=d)
Definition: remez_gmp.h:62
void stpini(multi1d< bigfloat > &step)
Initialise step sizes.
Definition: remez_gmp.cc:260
multi1d< bigfloat > poles
Definition: remez_gmp.h:58
void pfe(multi1d< bigfloat > &res, multi1d< bigfloat > &poles, const bigfloat &norm)
Definition: remez_gmp.cc:683
multi1d< bigfloat > step
Definition: remez_gmp.h:76
bigfloat polyEval(const bigfloat &x, const multi1d< bigfloat > &poly, long size)
Evaluate the polynomial.
Definition: remez_gmp.cc:638
multi1d< bigfloat > param
Definition: remez_gmp.h:58
RemezGMP()
Default constructor.
Definition: remez_gmp.cc:13
RemezCoeff_t getPFE()
Return the partial fraction expansion of the approximation x^(pnum/pden)
Definition: remez_gmp.cc:149
unsigned long power_num
the numerator and denominator of the power we are approximating
Definition: remez_gmp.h:68
bigfloat apstrt
The bounds of the approximation.
Definition: remez_gmp.h:65
bigfloat spread
Definition: remez_gmp.h:77
void search(multi1d< bigfloat > &step)
Search for error maxima and minima.
Definition: remez_gmp.cc:277
void setBounds(const Real &lower, const Real &upper)
Reset the bounds of the approximation.
Definition: remez_gmp.cc:71
long prec
The precision of the GNU MP library.
Definition: remez_gmp.h:83
int root()
Calculate the roots of the approximation.
Definition: remez_gmp.cc:588
Real generateApprox(int num_degree, int den_degree, unsigned long power_num, unsigned long power_den)
Generate the rational approximation x^(pnum/pden)
Definition: remez_gmp.cc:89
bigfloat apend
Definition: remez_gmp.h:65
void equations()
Solve the equations.
Definition: remez_gmp.cc:367
Bigfloat.
Definition: bigfloat.h:23
static void setDefaultPrecision(unsigned long dprec)
Definition: bigfloat.h:36
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
int y
Definition: meslate.cc:35
int z
Definition: meslate.cc:36
int x
Definition: meslate.cc:34
Double q
Definition: mesq.cc:17
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
int i
Definition: pbg5p_w.cc:55
Complex a
Definition: invbicg.cc:95
START_CODE()
A(A, psi, r, Ncb, PLUS)
int k
Definition: invbicg.cc:119
int l
Definition: pade_trln_w.cc:111
Double sum
Definition: qtopcor.cc:37
int norm
Definition: qtopcor.cc:35
#define JMAX
Definition: remez_gmp.cc:7
Remez algorithm for finding nth roots.
Convenient structure to package Remez coeffs.
Definition: remez_coeff.h:19
multi1d< Real > pole
Definition: remez_coeff.h:21
multi1d< Real > res
Definition: remez_coeff.h:20
#define M_PI