sturm.cpp

Go to the documentation of this file.
00001 /*****************************************************************************
00002     
00003     sturm.cpp -- See declarations in sturm.hpp.
00004 
00005     This file is a part of the Arageli library.
00006 
00007     Copyright (C) Nikolai Yu. Zolotykh, 1999--2006
00008     Copyright (C) Sergey S. Lyalin, 2005 -- 2006
00009     University of Nizhni Novgorod, Russia -- 2006
00010 
00011 *****************************************************************************/
00012 
00013 #include "config.hpp"
00014 
00015 #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
00016     defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_STURM)
00017 
00018 #include <iostream> // only for debug output
00019 #include <list>
00020 #include <cstdlib>
00021 
00022 #include "interval.hpp"
00023 #include "cmp.hpp"
00024 #include "gcd.hpp"
00025 #include "polyalg.hpp"
00026 
00027 #include "sturm.hpp"
00028 
00029 
00030 namespace Arageli
00031 {
00032 
00033 namespace _Internal
00034 {
00035 
00036 template <typename T, typename N>
00037 struct sturm_segment
00038 {
00039     sturm_segment () {}
00040 
00041     sturm_segment
00042     (const T& limfirst, const T& limsecond, const N& changefirst, const N& changesecond)
00043     : lims(limfirst, limsecond), changes(changefirst, changesecond) {}
00044     
00045     interval<T> lims;
00046     interval<N> changes;
00047 };
00048 
00049 template <typename Out, typename T, typename N>
00050 Out& operator<< (Out& out, const sturm_segment<T, N>& x)
00051 {
00052     out << "(" << x.lims << ", " << x.changes << ")";
00053     return out;
00054 }
00055 
00056 }
00057 
00058 
00059 template <typename P, typename SS>
00060 void sturm_diff_sys (const P& p, SS& ss)
00061 {
00062     ss.reserve(p.degree() + 1);
00063 
00064     ss.push_back(p);
00065     ss.push_back(polynom_primpart(diff(p)));
00066     
00067     for(std::size_t i = 2; !ss[i-1].is_const(); ++i)
00068     {
00069         P pr, pq;
00070         typename P::coef_type mult;
00071         polynom_pseudodivide_simple(-ss[i-2], ss[i-1], pq, pr, mult);
00072         ss.push_back(polynom_primpart(pr));
00073         //ss.push_back(-ss[i-2] % ss[i-1]);
00074         
00075         if(ss[i].is_null())
00076         {
00077             ss.pop_back();
00078             ARAGELI_ASSERT_2(is_null(ss % safe_reference(ss.back())));
00079             ss /= safe_reference(ss.back());
00080             break;
00081         }
00082     }
00083 }
00084 
00085 
00086 template <typename SS, typename T>
00087 typename SS::size_type sturm_sign_changes (const SS& ss, const T& x, int signpx)
00088 {
00089     ARAGELI_ASSERT_0(!ss.is_empty());
00090     
00091     if(ss.size() == 1)return 0;
00092 
00093     typename SS::size_type changes = 0;
00094     
00095     for(typename SS::const_iterator i = ss.begin() + 1; i != ss.end(); ++i)
00096         if(int cursign = sign(i->subs(x)))
00097         {
00098             if(signpx != 0 && cursign != signpx)
00099                 changes++;
00100             signpx = cursign;
00101         }
00102 
00103     return changes;
00104 }
00105 
00106 
00107 template <typename SS, typename Seg>
00108 typename SS::size_type sturm_number_roots (const SS& ss, const Seg& seg)
00109 {
00110     ARAGELI_ASSERT_0(!ss.is_empty());
00111 
00112     typedef typename SS::size_type size_type;
00113     
00114     switch(seg.kind())
00115     {
00116         case interval_empty: return 0;
00117         case interval_point: return is_null(ss.front().subs(seg.left()));
00118         case interval_length:
00119         {
00120             int lsign = sign(ss.front().subs(seg.left()));
00121             size_type lchanges = sturm_sign_changes(ss, seg.left(), lsign);
00122             if(lchanges == 0)return is_null(lsign);
00123             size_type rchanges = sturm_sign_changes(ss, seg.right());
00124             ARAGELI_ASSERT_0(lchanges >= rchanges);
00125             return lchanges - rchanges + is_null(lsign);
00126         }
00127         default: ARAGELI_ASSERT_1(!"It is impossible.");
00128     }
00129 }
00130 
00131 
00132 template <typename SS, typename SegT, bool SegREFCNT>
00133 vector<typename SS::size_type> sturm_number_roots
00134 (
00135     const SS& ss,
00136     const vector<SegT, SegREFCNT>& lims
00137 )
00138 {
00139     vector<typename SS::size_type> res(lims.size());
00140     for(std::size_t i = 0; i < res.size(); ++i)
00141         res[i] = sturm_number_roots(ss, lims[i]);
00142     return res;
00143 }
00144 
00145 
00146 template <typename T, typename P, typename LIMS, typename SegBound>
00147 void sturm (const P& p, LIMS& lims, SegBound bs)
00148 {
00149     if(p.degree() <= 0)
00150     {
00151         lims.resize(0);
00152         return;
00153     }
00154 
00155     // Building of the Sturm system.
00156 
00157     typedef vector<P, false> SS;
00158     SS ss;
00159     sturm_diff_sys(p, ss);
00160 
00161     //output_aligned(cout << "\nSturm's system:\n", ss);
00162 
00163     typedef _Internal::sturm_segment<T, std::size_t> SSEG;
00164     SSEG t1;
00165 
00166     if(bs.first() > bs.second())
00167     {
00168         t1.lims.second() = root_upper_bound_cauchy<T>(p);
00169         t1.lims.first() = -t1.lims.second();
00170     }
00171     else
00172     {
00173         t1.lims.second() = bs.second();
00174         t1.lims.first() = bs.first();
00175     }
00176 
00177     // Copmutes sign changes at the ends of investigated segment.
00178 
00179     t1.changes.first() = sturm_sign_changes(ss, t1.lims.first());
00180     t1.changes.second() = sturm_sign_changes(ss, t1.lims.second());
00181 
00182     ARAGELI_ASSERT_1(t1.changes.first() >= t1.changes.second());
00183 
00184     if(t1.changes.first() == t1.changes.second())
00185     {// polynomial hasn't real roots
00186         lims.resize(0);
00187         return;
00188     }
00189 
00190     if(t1.changes.first() - t1.changes.second() == 1)
00191     {// polynomial has only one root
00192         lims.assign(1, t1.lims);
00193         return;
00194     }
00195 
00196     typedef std::list<SSEG> LSSEG;
00197     LSSEG lsseg;
00198 
00199     lsseg.push_back(t1);
00200     
00201 
00202     for(typename LSSEG::iterator lssegi = lsseg.begin(); lssegi != lsseg.end();)
00203     {
00204         if
00205         (
00206             lssegi->lims.first() == lssegi->lims.second() ||    // exact root
00207             lssegi->changes.first() - lssegi->changes.second() == 1 // one root
00208         )
00209         {
00210             ++lssegi;
00211             continue;
00212         }
00213 
00214         //std::cout << "\nlssegi->lims = " << lssegi->lims;
00215         //std::cout << "\nlssegi->changes = " << lssegi->changes;
00216 
00217         ARAGELI_ASSERT_1(lssegi->lims.first() < lssegi->lims.second());
00218         ARAGELI_ASSERT_1(lssegi->changes.first() - lssegi->changes.second() >= 2);
00219 
00220         T mid = (lssegi->lims.first() + lssegi->lims.second())/2;
00221 
00222         //cout << "\nmid = " << mid;
00223 
00224         ARAGELI_ASSERT_1(lssegi->lims.first() < mid);
00225         ARAGELI_ASSERT_1(mid < lssegi->lims.second());
00226 
00227         bool exact_root = is_null(ss[0].subs(mid));
00228         int changes = sturm_sign_changes(ss, mid);
00229 
00230         typename LSSEG::iterator new_lssegi = lssegi;
00231 
00232         if(lssegi->changes.first() == changes)
00233         {// there are no roots
00234             lssegi->lims.first() = mid;
00235         }
00236         else
00237         {
00238             new_lssegi = lsseg.insert
00239             (
00240                 lssegi,
00241                 SSEG
00242                 (
00243                     lssegi->lims.first(), mid,
00244                     lssegi->changes.first(), changes
00245                 )
00246             );
00247 
00248             lssegi->lims.first() = mid;
00249             lssegi->changes.first() = changes;
00250         }
00251 
00252         if(lssegi->changes.first() - changes == 1)
00253             new_lssegi = lssegi;    // pass segment with one root
00254 
00255         if(exact_root)
00256             lsseg.insert(lssegi, SSEG(mid, mid, changes, changes));
00257 
00258         if(changes == lssegi->changes.second())
00259         {
00260             if(new_lssegi == lssegi)
00261                 new_lssegi = lsseg.erase(lssegi);
00262             else
00263                 lsseg.erase(lssegi);
00264         }
00265 
00266         lssegi = new_lssegi;
00267     }
00268 
00269     // Here segments can overlapped and share roots.
00270     // Exact root location.
00271 
00272     //output_aligned
00273     //(
00274     //  std::cout << "\nlsseg =\n",
00275     //  vector<SSEG>(lsseg.size(), lsseg.begin(), fromseq)
00276     //);
00277 
00278     {
00279         typename LSSEG::iterator
00280             lssegi = lsseg.begin();
00281 
00282         for(;;)
00283         {
00284             typename LSSEG::iterator lssegend = lsseg.end();
00285             if(lssegi == lssegend)break;
00286             --lssegend;
00287             if(lssegi == lssegend)break;
00288 
00289             ARAGELI_ASSERT_1
00290             (
00291                 !is_null(p.subs(lssegi->lims.first())) ||
00292                 lssegi->changes.is_point()
00293             );
00294 
00295             if(lssegi->changes.is_point())
00296             {
00297                 typename LSSEG::iterator next = lssegi;
00298                 ++next;
00299                 ARAGELI_ASSERT_1(next != lsseg.end());
00300                 if(next->lims.first() == lssegi->lims.second())
00301                     if(next->changes.is_point())
00302                         lsseg.erase(next);
00303                     else
00304                     {
00305                         int lsign = -sign(p.subs(next->lims.second()));
00306                         ARAGELI_ASSERT_1(lsign);
00307                         
00308                         do
00309                         {
00310                             if(interval_root_dichotomy(p, lsign, next->lims))
00311                                 next->changes.first() = next->changes.second();
00312                         }while(next->lims.first() == lssegi->lims.second());
00313                     }
00314             }
00315             else
00316             {
00317                 typename LSSEG::iterator next = lssegi;
00318                 ++next;
00319                 ARAGELI_ASSERT_1(next != lsseg.end());
00320                 if(next->lims.first() == lssegi->lims.second())
00321                     if(is_null(p.subs(lssegi->lims.second())))
00322                     {
00323                         ARAGELI_ASSERT_1(next->changes.is_point());
00324                         lssegi = lsseg.erase(lssegi);
00325                         continue;
00326                     }
00327                     else
00328                     {
00329                         int lsign = sign(p.subs(lssegi->lims.first()));
00330                         ARAGELI_ASSERT_1(lsign);
00331 
00332                         do
00333                         {
00334                             //std::cout << " " << lssegi->lims << " ";
00335                             if(interval_root_dichotomy(p, lsign, lssegi->lims))
00336                                 lssegi->changes.first() = lssegi->changes.second();
00337                         }while(next->lims.first() == lssegi->lims.second());
00338                     }
00339             }
00340 
00341             ++lssegi;
00342         }
00343 
00344     }
00345 
00346     //output_aligned
00347     //(
00348     //  std::cout << "\nlsseg =\n",
00349     //  vector<SSEG>(lsseg.size(), lsseg.begin(), fromseq)
00350     //);
00351     
00352     ARAGELI_ASSERT_1(lsseg.size() <= p.degree());
00353     lims.resize(lsseg.size());
00354 
00355     typename LSSEG::const_iterator lssegi = lsseg.begin();
00356     for(std::size_t i = 0; i < lims.size(); ++i, ++lssegi)
00357         lims[i] = lssegi->lims;
00358 
00359     ARAGELI_ASSERT_2(is_unit(sturm_number_roots(ss, lims)));
00360 }
00361 
00362 
00363 template <typename P, typename Lims>
00364 bool interval_root_dichotomy (const P& p, int lsign, Lims& lims)
00365 {
00366     typedef typename Lims::value_type T;
00367 
00368     T mid = (lims.first() + lims.second())/2;
00369     ARAGELI_ASSERT_0(lims.first() < mid);
00370     ARAGELI_ASSERT_0(mid < lims.second());
00371     int msign = sign(p.subs(mid));
00372 
00373     if(is_null(msign))
00374     {// exact root
00375         lims.first() = lims.second() = mid;
00376         return true;
00377     }
00378     else if(lsign == msign)
00379         lims.first() = mid;
00380     else
00381         lims.second() = mid;
00382 
00383     return false;
00384 }
00385 
00386 
00387 template <typename P, typename Lims, typename T>
00388 bool interval_root_precise (const P& p, Lims& lims, const T& e)
00389 {
00390     T re = lims.length();   // real error
00391     if(is_null(re))return true;
00392 
00393     if(lims.length() > e)
00394     {
00395         int lsign = sign(p.subs(lims.first()));
00396 
00397         do
00398         {
00399             if(interval_root_dichotomy(p, lsign, lims))return true;
00400         }
00401         while(lims.length() > e);
00402     }
00403 
00404     return false;
00405 }
00406 
00407 
00408 template <typename T, typename P, typename Roots, typename Prec>
00409 void roots_poly_real
00410 (
00411     const P& p, Roots& roots,
00412     Prec& prec, const T& e
00413 )
00414 {
00415     ARAGELI_ASSERT_0(!p.is_null());
00416     typedef vector<interval<T> > Lims;
00417     Lims lims;
00418     sturm<T>(p, lims);
00419 
00420     // computes roots more exactly
00421     
00422     if(!is_null(e))
00423         for(typename Lims::iterator i = lims.begin(); i != lims.end(); ++i)
00424             interval_root_precise(p, *i, e);
00425 
00426     roots.resize(lims.size());
00427     prec.resize(lims.size());
00428 
00429     for(std::size_t i = 0; i < lims.size(); ++i)
00430     {
00431         roots[i] = lims[i].first();
00432         prec[i] = lims[i].length();
00433     }
00434 }
00435 
00436 
00437 }
00438 
00439 #endif // #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE

Generated on Thu Aug 31 17:38:11 2006 for Arageli by  doxygen 1.4.7