src/2geom/poly-laguerre-solve.cpp

   1 #include <2geom/poly-laguerre-solve.h>
   2 #include <iterator>
   3
   4
   5 namespace Geom {
   6
   7 typedef std::complex<double> cdouble;
   8
   9 cdouble laguerre_internal_complex(Poly const & p,
  10                                   double x0,
  11                                   double tol,
  12                                   bool & quad_root) {
  13     cdouble a = 2*tol;
  14     cdouble xk = x0;
  15     double n = p.degree();
  16     quad_root = false;
  17     const unsigned shuffle_rate = 10;
  18     //static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
  19     unsigned shuffle_counter = 0;
  20     while(std::norm(a) > (tol*tol)) {
  21         //std::cout << "xk = " << xk << std::endl;
  22         cdouble b = p.back();
  23         cdouble d = 0, f = 0;
  24         double err = abs(b);
  25         double abx = abs(xk);
  26         for(int j = p.size()-2; j >= 0; j--) {
  27             f = xk*f + d;
  28             d = xk*d + b;
  29             b = xk*b + p[j];
  30             err = abs(b) + abx*err;
  31         }
  32
  33         err *= 1e-7; // magic epsilon for convergence, should be computed from tol
  34
  35         cdouble px = b;
  36         if(abs(b) < err)
  37             return xk;
  38         //if(std::norm(px) < tol*tol)
  39         //    return xk;
  40         cdouble G = d / px;
  41         cdouble H = G*G - f / px;
  42
  43         //std::cout << "G = " << G << "H = " << H;
  44         cdouble radicand = (n - 1)*(n*H-G*G);
  45         //assert(radicand.real() > 0);
  46         if(radicand.real() < 0)
  47             quad_root = true;
  48         //std::cout << "radicand = " << radicand << std::endl;
  49         if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
  50             a = - sqrt(radicand);
  51         else
  52             a = sqrt(radicand);
  53         //std::cout << "a = " << a << std::endl;
  54         a = n / (a + G);
  55         //std::cout << "a = " << a << std::endl;
  56         //if(shuffle_counter % shuffle_rate == 0)
  57             //a *= shuffle[shuffle_counter / shuffle_rate];
  58         xk -= a;
  59         shuffle_counter++;
  60         if(shuffle_counter >= 90)
  61             break;
  62     }
  63     //std::cout << "xk = " << xk << std::endl;
  64     return xk;
  65 }
  66
  67 double laguerre_internal(Poly const & p,
  68                          Poly const & pp,
  69                          Poly const & ppp,
  70                          double x0,
  71                          double tol,
  72                          bool & quad_root) {
  73     double a = 2*tol;
  74     double xk = x0;
  75     double n = p.degree();
  76     quad_root = false;
  77     while(a*a > (tol*tol)) {
  78         //std::cout << "xk = " << xk << std::endl;
  79         double px = p(xk);
  80         if(px*px < tol*tol)
  81             return xk;
  82         double G = pp(xk) / px;
  83         double H = G*G - ppp(xk) / px;
  84
  85         //std::cout << "G = " << G << "H = " << H;
  86         double radicand = (n - 1)*(n*H-G*G);
  87         assert(radicand > 0);
  88         //std::cout << "radicand = " << radicand << std::endl;
  89         if(G < 0) // here we try to maximise the denominator avoiding cancellation
  90             a = - sqrt(radicand);
  91         else
  92             a = sqrt(radicand);
  93         //std::cout << "a = " << a << std::endl;
  94         a = n / (a + G);
  95         //std::cout << "a = " << a << std::endl;
  96         xk -= a;
  97     }
  98     //std::cout << "xk = " << xk << std::endl;
  99     return xk;
 100 }
 101
 102
 103 std::vector<cdouble >
 104 laguerre(Poly p, const double tol) {
 105     std::vector<cdouble > solutions;
 106     //std::cout << "p = " << p << " = ";
 107     while(p.size() > 1)
 108     {
 109         double x0 = 0;
 110         bool quad_root = false;
 111         cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
 112         //if(abs(sol) > 1) break;
 113         Poly dvs;
 114         if(quad_root) {
 115             dvs.push_back((sol*conj(sol)).real());
 116             dvs.push_back(-(sol + conj(sol)).real());
 117             dvs.push_back(1.0);
 118             //std::cout << "(" <<  dvs << ")";
 119             //solutions.push_back(sol);
 120             //solutions.push_back(conj(sol));
 121         } else {
 122             //std::cout << sol << std::endl;
 123             dvs.push_back(-sol.real());
 124             dvs.push_back(1.0);
 125             solutions.push_back(sol);
 126             //std::cout << "(" <<  dvs << ")";
 127         }
 128         Poly r;
 129         p = divide(p, dvs, r);
 130         //std::cout << r << std::endl;
 131     }
 132     return solutions;
 133 }
 134
 135 std::vector<double>
 136 laguerre_real_interval(Poly const & /*ply*/,
 137                        const double /*lo*/, const double /*hi*/,
 138                        const double /*tol*/)
 139 {
 140     /* not implemented*/
 141     assert(false);
 142     std::vector<double> result;
 143     return result;
 144 }
 145
 146 } // namespace Geom
 147
 148 /*
 149   Local Variables:
 150   mode:c++
 151   c-file-style:"stroustrup"
 152   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 153   indent-tabs-mode:nil
 154   fill-column:99
 155   End:
 156 */
 157 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :