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

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