smithpoly.cpp

Go to the documentation of this file.
00001 /* /////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  File:       smithpoly.cpp
00004 //  Created:    2005/10/19    19:31
00005 //
00006 //  Author: Andrey Somsikov
00007 */
00008 
00009 #include "config.hpp"
00010 
00011 #if defined(ARAGELI_USE_COUNTERS)
00012 #include "counter.hpp"
00013 #endif
00014 
00015 #include <cstddef>
00016 #include <vector>
00017 
00018 #include "vector.hpp"
00019 
00020 //#define GCD_LOG
00021 
00022 namespace std
00023 {
00024 template <class T>
00025 std::ostream& operator << (std::ostream & s, std::vector<T> &v)
00026 {
00027     std::size_t i = v.size();
00028     s << "(";
00029     do {
00030         i--;
00031         s << v[i] << ", ";              
00032     } while(i != 0);
00033     s << ")";
00034     return s;
00035 }
00036 }
00037 
00038 // The following block was commented by Sergey Lyalin.
00039 //template <class T>
00040 //std::ostream& operator << (std::ostream & s, Arageli::vector<T> &v)
00041 //{
00042 //    std::size_t i = v.size();
00043 //    s << "(";
00044 //    do {
00045 //        i--;
00046 //        s << v[i] << ", ";                
00047 //    } while(i != 0);
00048 //    s << ")";
00049 //    return s;
00050 //}
00051 
00052 namespace Arageli
00053 {
00054 #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_smithpoly
00055 #if defined(ARAGELI_USE_COUNTERS)
00056     CCounter counter_big_int_sum_GC;
00057     CCounter counter_big_int_mul_GC;
00058     CCounter counter_big_int_div_GC;
00059 #endif
00060 #endif
00061 };
00062 
00063 #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
00064     defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_smithpoly)
00065 
00066 #include "smithpoly.hpp"
00067 #include "hermite.hpp"
00068 #include "rand.hpp"
00069 #include "sparse_polynom.hpp"
00070 #include <stdlib.h>
00071 #include <iostream>
00072 
00073 namespace Arageli
00074 {
00075 
00076 template <typename T>
00077 T normalize(const T &value)
00078 {
00079     if (0 == value.size())
00080         return value;
00081     T res;
00082     T::coef_type d = value.leading_monom().coef();
00083     for (T::monom_const_iterator it = value.monoms_begin();
00084         it != value.monoms_end(); it++)
00085         res += *it / T::monom(d);
00086     return res;
00087 }
00088 
00089 template <typename P>
00090 P gcd_gauss(P &f, P &g, 
00091             P *u = NULL, P *v = NULL)
00092 {
00093     typedef typename P::coef_type PT;
00094 
00095 #ifdef USE_POLYGCD_COUNTER
00096     counter_polygcd_count++;
00097     counter_polygcd_maxdeg += f.deg();
00098     counter_polygcd_maxdeg += g.deg();
00099 #endif
00100 
00101 #ifdef GCD_LOG
00102     std::cout << "START GCD COMPUTATION\n"
00103         << "f=" << f << "\n"
00104         << "g=" << g << "\n";
00105 #endif
00106 //    if (P(factory<PT>::null()) == f || P(factory<PT>::null()) == g)
00107     if (is_null(f.degree()) || 
00108         is_null(g.degree()))
00109     {
00110         if (is_null(f) && !is_null(g.degree()))
00111             return normalize(g);
00112         else if(is_null(g) && !is_null(f.degree()))
00113             return normalize(f);
00114 
00115         return P(factory<PT>::unit());
00116     }
00117         // should fill u, v;
00118     std::size_t i, j, k;
00119     std::size_t dSize = f.degree() + g.degree();
00120     matrix<PT> m(dSize, dSize);
00121     std::vector<PT> d(dSize);
00122     vector<std::size_t> uvColIndex(dSize);
00123 
00124     P::degree_type commonMult = std::min(
00125         f.monoms_begin()->degree(),
00126         g.monoms_begin()->degree());
00127 /*
00128     commonMult = 0;
00129     while (is_null(f[commonMult]) && is_null(g[commonMult]))
00130         commonMult++;
00131 */
00132 
00133     for (i = 0; i < m.nrows(); i++)
00134         for (j = 0; j < m.ncols(); j++)
00135             m(i, j) = factory<PT>::null();
00136 
00137     P::monom_iterator it;
00138 
00139     it = f.monoms_begin();
00140     for (i = 0; i < f.degree()+1; i++)
00141     {
00142         P::coef_type coef;
00143         if (it->degree() == P::degree_type(i))
00144         {
00145             coef = it->coef();
00146             it++;
00147         }
00148         else
00149             coef = factory<P::coef_type>::null();
00150         for (j = 0; j < g.degree(); j++)
00151             m(i + j, j) = coef;
00152     }
00153     it = g.monoms_begin();
00154     for (i = 0; i < g.degree()+1; i++)
00155     {
00156         P::coef_type coef;
00157         if (it->degree() == P::degree_type(i))
00158         {
00159             coef = it->coef();
00160             it++;
00161         }
00162         else
00163             coef = factory<P::coef_type>::null();
00164         for (j = 0; j < f.degree(); j++)
00165             m(i + j, g.degree() + j) = coef;
00166     }
00167 
00168 #ifdef GCD_LOG
00169 std::cout << "m = \n";
00170 output_aligned(std::cout, m, "|| ", " ||", "  ");
00171 std::cout << "\n";
00172 #endif
00173     matrix<PT> mTmp = m;
00174 
00175     PT mMult = factory<PT>::unit();
00176     i = m.nrows() - 1; // row
00177     for (;;)
00178     {
00179         // find left most non zero element in row i
00180         for (j = 0; j < m.ncols(); j++)
00181             if (!is_null(m.el(i, j)))
00182                 break;
00183         if (m.ncols() == j)
00184         {
00185             i++;
00186             // deg(d)=m.nrows()- i; coefficients of 
00187             //the u(x) and v(x) can be found from linear system
00188             break;
00189         }
00190         uvColIndex[i] = j;
00191         // d = (0, ..., 0, d_i:=1, 0, ..., 0)
00192 
00193 #ifdef GCD_LOG
00194 std::cout << "leading column = " << j << "\n";
00195 #endif
00196         for (k = 0; k < i; k++)
00197         {
00198             PT tmp = factory<PT>::opposite_unit()*m.el(k, j);
00199             m.mult_row(k, m.el(i, j));
00200             m.addmult_rows(k, i, tmp);
00201             m.div_row(k, mMult);
00202             d[k] = tmp/mMult;
00203         }
00204         d[i] = factory<PT>::unit();
00205         for (k = i+1; k < m.nrows(); k++)
00206         {
00207             PT tmp = factory<PT>::opposite_unit()*m.el(k, j);
00208             m.mult_row(k, m.el(i, j));
00209             m.addmult_rows(k, i, tmp);
00210             m.div_row(k, mMult);
00211             d[k] = tmp/mMult;
00212         }
00213 #ifdef GCD_LOG
00214 std::cout << "m = \n";
00215 output_aligned(std::cout, m, "|| ", " ||", "  ");
00216 std::cout << "\n";
00217 std::cout << "d = " << d << "\n";
00218 #endif
00219         mMult = m.el(i, j);
00220 
00221         if (commonMult == i)
00222             return P(factory<PT>::unit()) * P(P::monom(factory<PT>::unit(), commonMult));
00223 //break;
00224         i--;
00225     }
00226 
00227     vector<PT> uv(dSize);
00228     for (j = 0; j < dSize; j++)
00229         uv[j] = factory<PT>::null();
00230     for (j = i; j < dSize; j++)
00231         uv[uvColIndex[j]] = d[j]/m.el(j, uvColIndex[j]);
00232 
00233     for (j = 0; j < dSize; j++)
00234     {
00235         PT val = factory<PT>::null();
00236         for (i = 0; i < dSize; i++)
00237             val += mTmp.el(j, i) * uv[i];
00238         d[j] = val; // !!! should be d = mTmp*uv
00239     }
00240 
00241 #ifdef GCD_LOG
00242 std::cout << "uvColIndex = " << uvColIndex << "\n";
00243 std::cout << "uv = " << uv << "\n";
00244 std::cout << "d = " << d << "\n";
00245 #endif    
00246     P dNorm;
00247     for (i = 0; i < d.size(); i++)
00248         dNorm += P::monom(d[i], i);//= P(d, true);
00249     normalize(dNorm);
00250 #ifdef GCD_LOG
00251 std::cout << "dNorm = " << dNorm << "\n";
00252 #endif
00253 
00254 #if defined(ARAGELI_USE_COUNTERS)
00255 counter_big_int.suspend();
00256 counter_big_int_sum.suspend();
00257 counter_big_int_mul.suspend();
00258 counter_big_int_div.suspend();
00259 #endif
00260     if (f.degree() < dNorm.degree() ||
00261         g.degree() < dNorm.degree() ||
00262         P(factory<PT>::null()) != f % dNorm ||
00263         P(factory<PT>::null()) != g % dNorm)
00264     {
00265         std::cout << "gcd check fail:\n" << 
00266             "f=" << f << "\n" << 
00267             "g=" << g << "\n" <<
00268             "f % dNorm=" << f % dNorm << "\n" << 
00269             "g % dNorm=" << g % dNorm << "\n" << 
00270             "dNorm=" << dNorm << "\n"; 
00271 
00272         throw "gcd check fail";
00273     }
00274 #if defined(ARAGELI_USE_COUNTERS)
00275     counter_big_int.resume();
00276     counter_big_int_sum.resume();
00277     counter_big_int_mul.resume();
00278     counter_big_int_div.resume();
00279 #endif
00280 /*
00281     if (NULL != u)
00282     {
00283         vector<T> uData();
00284         u->initLe(uv.);
00285     }
00286 */
00287     if (0 == d.size())
00288         return P(factory<PT>::unit());
00289     
00290     P res;
00291     for (i = 0; i < d.size(); i++)
00292         res += P::monom(d[i], i);//= P(d, true);
00293     normalize(res);
00294 
00295     return res;
00296 }
00297 
00298 template <typename T>
00299 vector<T> makeGoodConditioning(matrix<T>& b, std::size_t col)
00300 {
00301     typedef matrix<T> M;
00302     typedef typename M::value_type MV;
00303     typedef typename M::size_type size_type;
00304 
00305     size_type i, j;
00306     vector<MV> res(b.ncols() - col);
00307 /*
00308 std::cout << "GC\n";
00309 output_aligned(std::cout, b, "|| ", " ||", "  ");
00310 std::cout << "\n";
00311 std::cout << "col = " << col << "\n";
00312 */
00313     for (j = col + 1; j < b.ncols(); j++)
00314     {
00315         MV gTarget = b.el(col, col);
00316         for (i = col+1; i < b.nrows(); i++)
00317         {
00318 //std::cout << "gTarget" << gTarget << "b" << b.el(i, col) << "\n";
00319             gTarget = gcd_gauss(gTarget, b.el(i, col));
00320         }
00321         MV g(gTarget);
00322         for (i = col; i < b.nrows(); i++)
00323         {
00324 //std::cout << "gTarget" << gTarget << "b" << b.el(i, j) << "\n";
00325             gTarget = gcd_gauss(gTarget, b.el(i, j));
00326         }
00327         MV a = factory<MV>::null();
00328         while (g != gTarget)
00329         {
00330             a+=factory<MV>::unit();
00331             b.addmult_cols(col, j, factory<MV>::unit());
00332             g = b.el(col, col);
00333             for (i = col+1; i < b.nrows(); i++)
00334                 g = gcd_gauss(g, b.el(i, col));
00335 //std::cout << "gTarget" << gTarget << " g" << g << " b=\n";
00336 //output_aligned(std::cout, b, "|| ", " ||", "  ");
00337 //std::cout << "\n";
00338 //std::cout << "a = " << a << "\n";
00339         }
00340         res[j-1 - col] = a;
00341     }
00342     return res;
00343 }
00344 
00345 template <typename T>
00346 matrix<T> smithPoly_rand(const matrix<T>& m)
00347 {
00348     typedef matrix<T> M;
00349     typedef typename M::value_type MV;
00350     typedef typename M::size_type size_type;
00351 
00352     size_type step = 0;
00353     size_type i, j;
00354     bool nonZeroFound;
00355     M s(m);
00356     M u(m.nrows(), eye);
00357     M v(m.ncols(), eye);
00358 
00359     for (;;)
00360     {
00361 //        std::cout << step++ << " ";
00362         size_type d = 0;
00363         size_type mSize = s.nrows();
00364         for (i = 0; i < s.nrows(); i++)
00365             for(j = 0; j < s.ncols(); j++)
00366                 // TO DO: use max function here
00367                 if (d < s.el(i, j).degree())
00368                     d = s.el(i, j).degree();
00369         MV maxRand(100000);
00370         maxRand *= MV(2)*MV(d)*MV(mSize)*MV(mSize)*MV(mSize);
00371 
00372         M c(s.nrows(), eye);
00373         for (i = 0; i < c.nrows(); i++)
00374             for(j = i+1; j < c.ncols(); j++)
00375                 c.el(i, j) = rand(maxRand);
00376 std::cout << step << " preconditioning: \n";
00377 output_aligned(std::cout, c, "|| ", " ||", "  ");
00378 
00379         s = s*c;
00380         s = hermite(s);
00381 std::cout << "after ref: \n";
00382 output_aligned(std::cout, s, "|| ", " ||", "  ");
00383         for (i = 0; i < s.nrows(); i++)
00384         {
00385             MV::monom::coef_type tmp(s.el(i, i).leading_monom().coef());
00386             for(j = i; j < s.ncols(); j++)
00387                 if (!is_null(s.el(i, j)))
00388                     s.el(i, j) /= tmp;
00389         }
00390 
00391 std::cout << "after ref normalized: \n";
00392 output_aligned(std::cout, s, "|| ", " ||", "  ");
00393         s.transpose();
00394         s = hermite(s);
00395         s.transpose();
00396 std::cout << "after cef: \n";
00397 output_aligned(std::cout, s, "|| ", " ||", "  ");
00398 step++;
00399 
00400         bool isInSmith = true;
00401         for (i = 0; i < s.nrows(); i++)
00402             for(j = 0; j < s.ncols(); j++)
00403                 if (i != j && !is_null(s.el(i, j)))
00404                 {
00405                     isInSmith = false;
00406                     break;
00407                 }
00408         if (!isInSmith)
00409             continue;
00410         for (i = 1; i < s.nrows(); i++)
00411         {
00412             if (!is_null(s.el(i, i) % s.el(i-1, i-1)))
00413             {
00414                 isInSmith = false;
00415                 break;
00416             }
00417         }
00418         if (isInSmith)
00419             break;
00420     }
00421 
00422     return s;
00423 }
00424 
00425 template <typename T, bool REFCNT>
00426 matrix<T, REFCNT> smithPoly_GC(const matrix<T, REFCNT>& m)
00427 {
00428     typedef matrix<T, REFCNT> M;
00429     typedef typename M::value_type MV;
00430     typedef typename M::size_type size_type;
00431 
00432     size_type i, j;
00433     bool nonZeroFound;
00434     matrix<T, REFCNT> s(m);
00435     matrix<T, REFCNT> u(m.nrows(), eye);
00436     matrix<T, REFCNT> v(m.ncols(), eye);
00437 
00438     j = 0;
00439     for (;;)
00440     {               
00441 //std::cout << "value=\n";
00442 //output_aligned(std::cout, s, "|| ", " ||", "  ");
00443 //std::cout << "\nu=\n";
00444 //output_aligned(std::cout, u, "|| ", " ||", "  ");
00445 //std::cout << "\nv=\n";
00446 //output_aligned(std::cout, v, "|| ", " ||", "  ");
00447 //std::cout << "\n";
00448 #if defined(ARAGELI_USE_COUNTERS)
00449         std::size_t sum = counter_big_int_sum.getSum();
00450         std::size_t mul = counter_big_int_mul.getSum();
00451         std::size_t div = counter_big_int_div.getSum();
00452 #endif
00453         vector<T> gc = makeGoodConditioning(s, j);
00454 #if defined(ARAGELI_USE_COUNTERS)
00455         counter_big_int_sum_GC += counter_big_int_sum.getSum() - sum;
00456         counter_big_int_mul_GC += counter_big_int_mul.getSum() - mul;
00457         counter_big_int_div_GC += counter_big_int_div.getSum() - div;
00458 #endif
00459         /*
00460         for (i = 0; i < j; i++)
00461         v.el(i, j) = factory<MV>::null();
00462         v.el(j, j) = factory<MV>::unit();
00463         */
00464         for (i = j+1; i < s.nrows(); i++)
00465             v.el(i, j) += gc[i-(j+1)];
00466 //std::cout << "gc " << j << ":\n";
00467 //output_aligned(std::cout, s, "|| ", " ||", "  ");
00468 //std::cout << "\n";
00469 
00470         size_type iPivot = j;
00471         T valPivot = s.el(j, j);
00472         for (i = j+1; i < s.nrows(); i++)
00473         {
00474 //std::cout << s.el(i, j) << "\n";
00475             if ((s.el(i, j).degree() < valPivot.degree() || is_null(valPivot)) && 
00476                 !is_null(s.el(i, j)))
00477             {
00478                 iPivot = i;
00479                 valPivot = s.el(i, j);
00480             }
00481         }
00482 //std::cout << "pivot:" << iPivot << ", " << j << "=" << valPivot << "\n";
00483         s.swap_rows(j, iPivot);
00484         u.swap_rows(j, iPivot);
00485 
00486         nonZeroFound = false;
00487         for (i = j+1; i < s.nrows(); i++)
00488         {
00489             if (!is_null(s.el(i, j)))
00490             {
00491                 T mult(factory<MV>::opposite_unit() * (s.el(i, j) / s.el(j, j)));
00492 //std::cout << "step" << i << "," << j << "*" << mult << ": \n";
00493 //output_aligned(std::cout, s, "|| ", " ||", "  ");
00494 //std::cout << "\n";
00495                 s.addmult_rows(i, j, mult);
00496                 u.addmult_rows(i, j, mult);
00497                 nonZeroFound = true;
00498             }
00499         }
00500         if (nonZeroFound)
00501             continue;
00502 
00503         j++;
00504         if (j >= s.ncols())
00505             break;
00506     }
00507 
00508     for (i = 0; i < s.nrows(); i++)
00509     {
00510         for (j = i+1; j < s.ncols(); j++)
00511         {
00512             T mult(factory<MV>::opposite_unit() * (s.el(i, j) / s.el(i, i)));
00513             s.addmult_cols(j, i, mult);
00514             v.addmult_cols(j, i, mult);
00515         }
00516     }
00517 
00518 /*
00519 if (s != u*m*v)
00520 {
00521     std::cout << "check failed\n";
00522     output_aligned(std::cout, u, "|| ", " ||", "  ");
00523     output_aligned(std::cout, m, "|| ", " ||", "  ");
00524     output_aligned(std::cout, v, "|| ", " ||", "  ");
00525     output_aligned(std::cout, u*m*v, "|| ", " ||", "  ");
00526     throw "check failed";
00527 }
00528 */
00529 
00530     return s;
00531 }
00532 
00533 typedef std::pair<std::size_t, std::size_t> TSmithPivot;
00534 
00535 template <class T>
00536 TSmithPivot findSmithPivot(const matrix<T> &m, std::size_t row, std::size_t col)
00537 {
00538     typedef matrix<T> M;
00539     typedef typename M::value_type MV;
00540     std::size_t rowMin = row;
00541     std::size_t colMin = col;
00542     T valueMin = m.el(rowMin, colMin);
00543     for (std::size_t i = row; i < m.nrows(); i++)
00544         for (std::size_t j = col; j < m.ncols(); j++)
00545         {
00546             if ((m.el(i, j).degree() < valueMin.degree() || is_null(valueMin)) && 
00547                 !is_null(m.el(i, j)))
00548             {
00549                 rowMin = i;
00550                 colMin = j;
00551                 valueMin = m.el(i, j);
00552             }
00553         }
00554     return TSmithPivot(rowMin, colMin);
00555 }
00556 
00557 template <typename T>
00558 matrix<T> smithPoly_std(const matrix<T> &m)
00559 {
00560     typedef matrix<T> M;
00561     typedef typename M::value_type MV;
00562     typedef typename M::size_type size_type;
00563 
00564     bool nonZeroFound;
00565     matrix<T> s(m);
00566     matrix<T> u(m.nrows(), eye);
00567     matrix<T> v(m.ncols(), eye);
00568 
00569     size_type i, j;
00570     size_type k = 0;
00571     for (;;)
00572     {
00573 //        tout << s << det(s) << " = det\n";
00574         TSmithPivot pivot = findSmithPivot<T>(s, k, k);
00575         if (!(k == pivot.first && k == pivot.second))
00576         {
00577             s.swap_rows(k, pivot.first);
00578             u.swap_rows(k, pivot.first);
00579             s.swap_cols(k, pivot.second);
00580             v.swap_cols(k, pivot.second);
00581         }
00582 /*
00583         if (s.el(k, k) < factory<MV>::null())
00584         {
00585             s.mult_row(k, factory<MV>::opposite_unit());
00586             u.mult_row(k, factory<MV>::opposite_unit());
00587         }
00588 */
00589 
00590         nonZeroFound = false;
00591         for (i = k+1; i < s.nrows(); i++)
00592         {
00593             if (!is_null(s.el(i, k)))
00594             {
00595 //                T tmp(s.el(i, k) / s.el(k, k));
00596 //                T mult(factory<MV>::opposite_unit() * tmp);
00597                 T mult(factory<MV>::opposite_unit() * (s.el(i, k) / s.el(k, k)));
00598                 s.addmult_rows(i, k, mult);
00599                 u.addmult_rows(i, k, mult);
00600                 nonZeroFound = true;
00601             }
00602         }
00603         if (nonZeroFound)
00604             continue;
00605 
00606         nonZeroFound = false;
00607         for (i = k+1; i < s.ncols(); i++)
00608         {
00609             if (!is_null(s.el(k, i)))
00610             {
00611                 T mult(factory<MV>::opposite_unit() * (s.el(k, i) / s.el(k, k)));
00612                 s.addmult_cols(i, k, mult);
00613                 v.addmult_cols(i, k, mult);
00614                 nonZeroFound = true;
00615             }
00616         }
00617         if (nonZeroFound)
00618             continue;
00619 
00620         bool nonDivFound = false;
00621         for (i = k; i < s.nrows(); i++)
00622         {
00623             for (j = k; j < s.ncols(); j++)
00624                 if (!is_null(s.el(i,j) % s.el(k,k)))
00625                 {
00626                     nonDivFound = true;
00627                     s.addmult_cols(k, j, factory<MV>::unit());
00628                     v.addmult_cols(k, j, factory<MV>::unit());
00629                     break;
00630                 }
00631             if (nonDivFound)
00632                 break;
00633         }
00634 
00635         if (!nonDivFound)
00636         {
00637             k++;
00638             if (k == s.nrows() || k == s.ncols())
00639                 break;
00640         }
00641     }
00642 
00643 
00644 /*
00645 #if defined(ARAGELI_USE_COUNTERS)
00646     counter_big_int.suspend();
00647     counter_big_int_sum.suspend();
00648     counter_big_int_mul.suspend();
00649     counter_big_int_div.suspend();
00650 #endif
00651     if (s != u*m*v)
00652     {
00653         std::cout << "check failed\n";
00654         std::cout << u << m << v << u*m*v;
00655         throw "check failed";
00656     }
00657 #if defined(ARAGELI_USE_COUNTERS)
00658     counter_big_int.resume();
00659     counter_big_int_sum.resume();
00660     counter_big_int_mul.resume();
00661     counter_big_int_div.resume();
00662 #endif
00663 */
00664 
00665     return s;
00666 }
00667 
00668 } // namespace Arageli
00669 
00670 #endif // #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
00671 
00672 /* End of file smithpoly.cpp */

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