![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index 921ec3e3b36ae81112929dfd5f986cd214c4bd3a..5b59690b742aec0ef565959b2813472086831ea9 100644 (file)
-#include "poly-laguerre-solve.h"\r
-#include <iterator>\r
-\r
-typedef std::complex<double> cdouble;\r
-\r
-cdouble laguerre_internal_complex(Poly const & p,\r
- double x0,\r
- double tol,\r
- bool & quad_root) {\r
- cdouble a = 2*tol;\r
- cdouble xk = x0;\r
- double n = p.degree();\r
- quad_root = false;\r
- const unsigned shuffle_rate = 10;\r
- static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};\r
- unsigned shuffle_counter = 0;\r
- while(std::norm(a) > (tol*tol)) {\r
- //std::cout << "xk = " << xk << std::endl;\r
- cdouble b = p.back();\r
- cdouble d = 0, f = 0;\r
- double err = abs(b);\r
- double abx = abs(xk);\r
- for(int j = p.size()-2; j >= 0; j--) {\r
- f = xk*f + d;\r
- d = xk*d + b;\r
- b = xk*b + p[j];\r
- err = abs(b) + abx*err;\r
- }\r
-\r
- err *= 1e-7; // magic epsilon for convergence, should be computed from tol\r
-\r
- cdouble px = b;\r
- if(abs(b) < err)\r
- return xk;\r
- //if(std::norm(px) < tol*tol)\r
- // return xk;\r
- cdouble G = d / px;\r
- cdouble H = G*G - f / px;\r
-\r
- //std::cout << "G = " << G << "H = " << H;\r
- cdouble radicand = (n - 1)*(n*H-G*G);\r
- //assert(radicand.real() > 0);\r
- if(radicand.real() < 0)\r
- quad_root = true;\r
- //std::cout << "radicand = " << radicand << std::endl;\r
- if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation\r
- a = - sqrt(radicand);\r
- else\r
- a = sqrt(radicand);\r
- //std::cout << "a = " << a << std::endl;\r
- a = n / (a + G);\r
- //std::cout << "a = " << a << std::endl;\r
- if(shuffle_counter % shuffle_rate == 0)\r
- ;//a *= shuffle[shuffle_counter / shuffle_rate];\r
- xk -= a;\r
- shuffle_counter++;\r
- if(shuffle_counter >= 90)\r
- break;\r
- }\r
- //std::cout << "xk = " << xk << std::endl;\r
- return xk;\r
-}\r
-\r
-double laguerre_internal(Poly const & p,\r
- Poly const & pp,\r
- Poly const & ppp,\r
- double x0,\r
- double tol,\r
- bool & quad_root) {\r
- double a = 2*tol;\r
- double xk = x0;\r
- double n = p.degree();\r
- quad_root = false;\r
- while(a*a > (tol*tol)) {\r
- //std::cout << "xk = " << xk << std::endl;\r
- double px = p(xk);\r
- if(px*px < tol*tol)\r
- return xk;\r
- double G = pp(xk) / px;\r
- double H = G*G - ppp(xk) / px;\r
-\r
- //std::cout << "G = " << G << "H = " << H;\r
- double radicand = (n - 1)*(n*H-G*G);\r
- assert(radicand > 0);\r
- //std::cout << "radicand = " << radicand << std::endl;\r
- if(G < 0) // here we try to maximise the denominator avoiding cancellation\r
- a = - sqrt(radicand);\r
- else\r
- a = sqrt(radicand);\r
- //std::cout << "a = " << a << std::endl;\r
- a = n / (a + G);\r
- //std::cout << "a = " << a << std::endl;\r
- xk -= a;\r
- }\r
- //std::cout << "xk = " << xk << std::endl;\r
- return xk;\r
-}\r
-\r
-\r
-std::vector<cdouble >\r
-laguerre(Poly p, const double tol) {\r
- std::vector<cdouble > solutions;\r
- //std::cout << "p = " << p << " = ";\r
- while(p.size() > 1)\r
- {\r
- double x0 = 0;\r
- bool quad_root = false;\r
- cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);\r
- //if(abs(sol) > 1) break;\r
- Poly dvs;\r
- if(quad_root) {\r
- dvs.push_back((sol*conj(sol)).real());\r
- dvs.push_back(-(sol + conj(sol)).real());\r
- dvs.push_back(1.0);\r
- //std::cout << "(" << dvs << ")";\r
- //solutions.push_back(sol);\r
- //solutions.push_back(conj(sol));\r
- } else {\r
- //std::cout << sol << std::endl;\r
- dvs.push_back(-sol.real());\r
- dvs.push_back(1.0);\r
- solutions.push_back(sol);\r
- //std::cout << "(" << dvs << ")";\r
- }\r
- Poly r;\r
- p = divide(p, dvs, r);\r
- //std::cout << r << std::endl;\r
- }\r
- return solutions;\r
-}\r
-\r
-std::vector<double>\r
-laguerre_real_interval(Poly const & ply,\r
- const double lo, const double hi,\r
- const double tol) {\r
-}\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+#include <2geom/poly-laguerre-solve.h>
+#include <iterator>
+
+
+namespace Geom {
+
+typedef std::complex<double> cdouble;
+
+cdouble laguerre_internal_complex(Poly const & p,
+ double x0,
+ double tol,
+ bool & quad_root) {
+ cdouble a = 2*tol;
+ cdouble xk = x0;
+ double n = p.degree();
+ quad_root = false;
+ const unsigned shuffle_rate = 10;
+ //static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
+ unsigned shuffle_counter = 0;
+ while(std::norm(a) > (tol*tol)) {
+ //std::cout << "xk = " << xk << std::endl;
+ cdouble b = p.back();
+ cdouble d = 0, f = 0;
+ double err = abs(b);
+ double abx = abs(xk);
+ for(int j = p.size()-2; j >= 0; j--) {
+ f = xk*f + d;
+ d = xk*d + b;
+ b = xk*b + p[j];
+ err = abs(b) + abx*err;
+ }
+
+ err *= 1e-7; // magic epsilon for convergence, should be computed from tol
+
+ cdouble px = b;
+ if(abs(b) < err)
+ return xk;
+ //if(std::norm(px) < tol*tol)
+ // return xk;
+ cdouble G = d / px;
+ cdouble H = G*G - f / px;
+
+ //std::cout << "G = " << G << "H = " << H;
+ cdouble radicand = (n - 1)*(n*H-G*G);
+ //assert(radicand.real() > 0);
+ if(radicand.real() < 0)
+ quad_root = true;
+ //std::cout << "radicand = " << radicand << std::endl;
+ if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
+ a = - sqrt(radicand);
+ else
+ a = sqrt(radicand);
+ //std::cout << "a = " << a << std::endl;
+ a = n / (a + G);
+ //std::cout << "a = " << a << std::endl;
+ //if(shuffle_counter % shuffle_rate == 0)
+ //a *= shuffle[shuffle_counter / shuffle_rate];
+ xk -= a;
+ shuffle_counter++;
+ if(shuffle_counter >= 90)
+ break;
+ }
+ //std::cout << "xk = " << xk << std::endl;
+ return xk;
+}
+
+double laguerre_internal(Poly const & p,
+ Poly const & pp,
+ Poly const & ppp,
+ double x0,
+ double tol,
+ bool & quad_root) {
+ double a = 2*tol;
+ double xk = x0;
+ double n = p.degree();
+ quad_root = false;
+ while(a*a > (tol*tol)) {
+ //std::cout << "xk = " << xk << std::endl;
+ double px = p(xk);
+ if(px*px < tol*tol)
+ return xk;
+ double G = pp(xk) / px;
+ double H = G*G - ppp(xk) / px;
+
+ //std::cout << "G = " << G << "H = " << H;
+ double radicand = (n - 1)*(n*H-G*G);
+ assert(radicand > 0);
+ //std::cout << "radicand = " << radicand << std::endl;
+ if(G < 0) // here we try to maximise the denominator avoiding cancellation
+ a = - sqrt(radicand);
+ else
+ a = sqrt(radicand);
+ //std::cout << "a = " << a << std::endl;
+ a = n / (a + G);
+ //std::cout << "a = " << a << std::endl;
+ xk -= a;
+ }
+ //std::cout << "xk = " << xk << std::endl;
+ return xk;
+}
+
+
+std::vector<cdouble >
+laguerre(Poly p, const double tol) {
+ std::vector<cdouble > solutions;
+ //std::cout << "p = " << p << " = ";
+ while(p.size() > 1)
+ {
+ double x0 = 0;
+ bool quad_root = false;
+ cdouble sol = laguerre_internal_complex(p, x0, tol, quad_root);
+ //if(abs(sol) > 1) break;
+ Poly dvs;
+ if(quad_root) {
+ dvs.push_back((sol*conj(sol)).real());
+ dvs.push_back(-(sol + conj(sol)).real());
+ dvs.push_back(1.0);
+ //std::cout << "(" << dvs << ")";
+ //solutions.push_back(sol);
+ //solutions.push_back(conj(sol));
+ } else {
+ //std::cout << sol << std::endl;
+ dvs.push_back(-sol.real());
+ dvs.push_back(1.0);
+ solutions.push_back(sol);
+ //std::cout << "(" << dvs << ")";
+ }
+ Poly r;
+ p = divide(p, dvs, r);
+ //std::cout << r << std::endl;
+ }
+ return solutions;
+}
+
+std::vector<double>
+laguerre_real_interval(Poly const & /*ply*/,
+ const double /*lo*/, const double /*hi*/,
+ const double /*tol*/)
+{
+ /* not implemented*/
+ assert(false);
+ std::vector<double> result;
+ return result;
+}
+
+} // namespace Geom
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :