Changeset 91788c0 in git


Ignore:
Timestamp:
Nov 29, 2011, 2:36:04 PM (12 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
eefc3a498c04d352c93280e8e8d0cafe97f4e754
Parents:
6deedd6650206f3f5a123f8dfbd41e449fd35f93
git-author:
Martin Lee <martinlee84@web.de>2011-11-29 14:36:04+01:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2011-12-12 18:00:42+01:00
Message:
add: code for sparse heuristic lifting a la Lucks/Wang
moved: code for other sparse heuristic to facSparseHensel
rm: unused function leadingCoeffReconstruction from facFqFactorize.cc
Location:
factory
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • factory/Makefile.am

    r6deedd r91788c0  
    6969                facHensel.cc \
    7070                facIrredTest.cc \
     71                facSparseHensel.cc \
    7172                fieldGCD.cc \
    7273                ffops.cc \
    … …  
    152153                facHensel.h \
    153154                facIrredTest.h \
     155                facSparseHensel.h \
    154156                fieldGCD.h \
    155157                ffops.h \

  • factory/facFactorize.cc

    r6deedd r91788c0  
    2929#include "algext.h"
    3030#include "cf_reval.h"
     31#include "facSparseHensel.h"
    3132
    3233#ifdef HAVE_NTL

  • factory/facFqFactorize.cc

    r6deedd r91788c0  
    3131#include "cf_map_ext.h"
    3232#include "algext.h"
     33#include "facSparseHensel.h"
    3334
    3435#ifdef HAVE_NTL
    … …  
    3940TIMING_DEFINE_PRINT(fac_hensel_lift)
    4041TIMING_DEFINE_PRINT(fac_factor_recombination)
    41 
    42 CFList
    43 recombination (const CFList& factors1, const CFList& factors2, int s, int thres,
    44                const CanonicalForm& evalPoint, const Variable& x);
    45 
    46 static
    47 int comp (const CanonicalForm& A, const CanonicalForm& B)
    48 {
    49   if (A.inCoeffDomain() && !B.inCoeffDomain())
    50     return -1;
    51   else if (!A.inCoeffDomain() && B.inCoeffDomain())
    52     return 1;
    53   else if (A.inCoeffDomain() && B.inCoeffDomain())
    54     return 0;
    55   else if (degree (A, 1) > degree (B, 1))
    56     return 1;
    57   else if (degree (A, 1) < degree (B, 1))
    58     return -1;
    59   // here A and B are not in CoeffDomain
    60   if (degree (A) > degree (B))
    61     return 1;
    62   else if (degree (A) < degree (B))
    63     return -1;
    64   return 0;
    65 }
    66 
    67 static
    68 void swap (CFArray& A, int i, int j)
    69 {
    70   CanonicalForm tmp= A[i];
    71   A[i]= A[j];
    72   A[j]= tmp;
    73 }
    74 
    75 static
    76 void quickSort (int lo, int hi, CFArray& A)
    77 {
    78   int i= lo, j= hi;
    79   CanonicalForm tmp= A[(lo+hi)/2];
    80   while (i <= j)
    81   {
    82     while (comp (A [i], tmp) < 0 && i < hi) i++;
    83     while (comp (tmp, A[j]) < 0 && j > lo) j--;
    84     if (i <= j)
    85     {
    86       swap (A, i, j);
    87       i++;
    88       j--;
    89     }
    90   }
    91   if (lo < j) quickSort (lo, j, A);
    92   if (i < hi) quickSort (i, hi, A);
    93 }
    94 
    95 static
    96 void sort (CFArray& A)
    97 {
    98   quickSort (0, A.size() - 1, A);
    99 }
    100 
    101 CFList
    102 findNormalizingFactor1 (const CFList& biFactors, const CanonicalForm& evalPoint,
    103                         CFList& uniFactors)
    104 {
    105   CFList result;
    106   CanonicalForm tmp;
    107   for (CFListIterator i= biFactors; i.hasItem(); i++)
    108   {
    109     tmp= i.getItem() (evalPoint);
    110     uniFactors.append (tmp /Lc (tmp));
    111     result.append (Lc (tmp));
    112   }
    113   return result;
    114 }
    115 
    116 CFList
    117 findNormalizingFactor2 (CFList& biFactors, const CanonicalForm& evalPoint,
    118                         const CFList& uniFactors)
    119 {
    120   CFList result;
    121   CFList uniBiFactors= biFactors;
    122   CFList newBiFactors;
    123   CFList l;
    124   int pos;
    125   CFListIterator iter;
    126   for (iter= uniBiFactors; iter.hasItem(); iter++)
    127   {
    128     iter.getItem()= iter.getItem() (evalPoint);
    129     l.append (Lc (iter.getItem()));
    130     iter.getItem() /= Lc (iter.getItem());
    131   }
    132   for (CFListIterator i= uniFactors; i.hasItem(); i++)
    133   {
    134     pos= findItem (uniBiFactors, i.getItem());
    135     newBiFactors.append (getItem (biFactors, pos));
    136     result.append (getItem (l, pos));
    137   }
    138   biFactors= newBiFactors;
    139   return result;
    140 }
    141 
    142 CFArray
    143 getTerms (const CanonicalForm& F)
    144 {
    145   CFArray result= CFArray (size (F));
    146   int j= 0;
    147   Variable x= F.mvar();
    148   Variable y= Variable (1);
    149   CFIterator k;
    150   for (CFIterator i= F; i.hasTerms(); i++)
    151   {
    152     for (k= i.coeff(); k.hasTerms(); k++, j++)
    153       result[j]= k.coeff()*power (x,i.exp())*power (y,k.exp());
    154   }
    155   sort (result);
    156   return result;
    157 }
    158 
    159 void getTerms (const CFList& F, CFArray* result)
    160 {
    161   int j= 0;
    162   for (CFListIterator i= F; i.hasItem(); i++, j++)
    163     result[j]= getTerms (i.getItem());
    164 }
    165 
    166 CFArray
    167 evaluate (const CFArray& A, const CanonicalForm& eval, const Variable& y)
    168 {
    169   int j= A.size();
    170   CFArray result= CFArray (j);
    171   for (int i= 0; i < j; i++)
    172     result [i]= A[i] (eval, y);
    173   return result;
    174 }
    175 
    176 CFArray*
    177 evaluate (CFArray* const& A, int sizeA, const CanonicalForm& eval,
    178           const Variable& y)
    179 {
    180   CFArray* result= new CFArray [sizeA];
    181   for (int i= 0; i < sizeA; i++)
    182     result[i]= evaluate (A[i], eval, y);
    183   return result;
    184 }
    185 
    186 CFList normalize (const CFList& L, const CFList& normalizingFactor)
    187 {
    188   CFList result;
    189   CFListIterator j= normalizingFactor;
    190   for (CFListIterator i= L; i.hasItem(); i++, j++)
    191     result.append (i.getItem() / j.getItem());
    192   return result;
    193 }
    194 
    195 int search (const CFArray& A, const CanonicalForm& F, int i, int j)
    196 {
    197   for (; i < j; i++)
    198     if (A[i] == F)
    199       return i;
    200   return -1;
    201 }
    202 
    203 CanonicalForm patch (const CanonicalForm& F1, const CanonicalForm& F2,
    204                      const CanonicalForm& eval)
    205 {
    206   CanonicalForm result= F1;
    207   if (F2.level() != 1)
    208   {
    209     int d= degree (F2);
    210     result *= power (F2.mvar(), d);
    211     result /= power (eval, d);
    212   }
    213   return result;
    214 }
    215 
    216 CFList
    217 sparseHeuristic (const CanonicalForm& A, const CFList& biFactors,
    218                  CFList*& moreBiFactors, const CFList& evaluation,
    219                  int minFactorsLength)
    220 {
    221   int j= A.level() - 1;
    222   int i;
    223 
    224   //initialize storage
    225   CFArray *** storeFactors= new CFArray** [j];
    226   for (i= 0; i < j; i++)
    227     storeFactors [i]= new CFArray* [2];
    228 
    229   CFArray eval= CFArray (j);
    230   i= j - 1;
    231   for (CFListIterator iter= evaluation; iter.hasItem(); iter++,i--)
    232     eval[i]= iter.getItem();
    233   storeFactors [0] [0]= new CFArray [minFactorsLength];
    234   storeFactors [0] [1]= new CFArray [minFactorsLength];
    235   for (i= 1; i < j; i++)
    236   {
    237     storeFactors[i] [0]= new CFArray [minFactorsLength];
    238     storeFactors[i] [1]= new CFArray [minFactorsLength];
    239   }
    240   //
    241 
    242   CFList * normalizingFactors= new CFList [j];
    243   CFList uniFactors;
    244   normalizingFactors [0]= findNormalizingFactor1 (biFactors,
    245                                               evaluation.getLast(), uniFactors);
    246   for (i= j - 1; i > 0; i--)
    247   {
    248     if (moreBiFactors[i-1].length() != minFactorsLength)
    249     {
    250       moreBiFactors[i-1]=
    251         recombination (moreBiFactors [i-1], uniFactors, 1,
    252                        moreBiFactors[i-1].length()-uniFactors.length()+1,
    253                        eval[i], Variable (i + 2)
    254                       );
    255     }
    256     normalizingFactors [i]= findNormalizingFactor2 (moreBiFactors [i - 1],
    257                                                     eval[i], uniFactors);
    258   }
    259 
    260   CFList tmp;
    261   tmp= normalize (biFactors, normalizingFactors[0]);
    262   getTerms (tmp, storeFactors [0] [0]);
    263   storeFactors [0] [1]= evaluate (storeFactors [0] [0], minFactorsLength,
    264                                   evaluation.getLast(), Variable (2));
    265   for (i= j - 1; i > 0; i--)
    266   {
    267     tmp= normalize (moreBiFactors [i-1], normalizingFactors [i]);
    268     getTerms (tmp, storeFactors [i] [0]);
    269     storeFactors [i] [1]= evaluate (storeFactors [i] [0], minFactorsLength,
    270                                     eval[i], Variable (i + 2));
    271   }
    272 
    273 
    274   int k, l, m, mm, count, sizeOfUniFactors= 0;
    275   int*** seperator= new int** [j];
    276   Variable x= Variable (1);
    277 
    278   for (i= 0; i < j; i++)
    279     seperator [i]= new int* [minFactorsLength];
    280   for (k= 0; k < minFactorsLength; k++)
    281   {
    282     for (i= 0; i < j; i++)
    283     {
    284       count= 0;
    285       for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
    286       {
    287         if (degree (storeFactors[i][0][k][l], x) <
    288             degree (storeFactors[i][0][k][l+1], x))
    289           count++;
    290       }
    291       if (i == 0)
    292         sizeOfUniFactors= count;
    293       else if (sizeOfUniFactors != count)
    294       {
    295         for (m= 0; m < j; m++)
    296         {
    297           delete [] storeFactors [m] [0];
    298           delete [] storeFactors [m] [1];
    299           delete [] storeFactors [m];
    300           for (mm= 0; mm < k; mm++)
    301             delete [] seperator [m][mm];
    302           delete [] seperator [m];
    303         }
    304         delete [] storeFactors;
    305         for (m= 0; m < i; m++)
    306           delete [] seperator [m][k];
    307         delete [] seperator;
    308         return CFList();
    309       }
    310       seperator [i][k]= new int [count + 3];
    311       seperator [i][k][0]= count + 1;
    312       seperator [i][k][1]= 0;
    313       count= 2;
    314       for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
    315       {
    316         if (degree (storeFactors[i][0][k][l], x) <
    317             degree (storeFactors[i][0][k][l+1], x))
    318         {
    319           seperator[i][k][count]=l + 1;
    320           count++;
    321         }
    322       }
    323       seperator [i][k][count]= storeFactors[i][0][k].size();
    324     }
    325   }
    326 
    327   CanonicalForm tmp1, factor, quot;
    328   CanonicalForm B= A;
    329   CFList result;
    330   int maxTerms, n, index1, index2, mmm, found, columns, oneCount;
    331   int ** mat;
    332 
    333   for (k= 0; k < minFactorsLength; k++)
    334   {
    335     factor= 0;
    336     sizeOfUniFactors= seperator [0][k][0];
    337     for (n= 1; n <= sizeOfUniFactors; n++)
    338     {
    339       columns= 0;
    340       maxTerms= 1;
    341       index1= j - 1;
    342       for (i= j - 1; i >= 0; i--)
    343       {
    344         if (maxTerms < seperator[i][k][n+1]-seperator[i][k][n])
    345         {
    346           maxTerms= seperator[i][k][n + 1]-seperator[i][k][n];
    347           index1= i;
    348         }
    349       }
    350       for (i= j - 1; i >= 0; i--)
    351       {
    352         if (i == index1)
    353           continue;
    354         columns += seperator [i][k][n+1]-seperator[i][k][n];
    355       }
    356       mat= new int *[maxTerms];
    357       mm= 0;
    358       for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
    359       {
    360         tmp1= storeFactors [index1][1][k][m];
    361         mat[mm]= new int [columns];
    362         for (i= 0; i < columns; i++)
    363           mat[mm][i]= 0;
    364         index2= 0;
    365         for (i= j - 1; i >= 0; i--)
    366         {
    367           if (i == index1)
    368             continue;
    369           found= -1;
    370           if ((found= search (storeFactors[i][1][k], tmp1,
    371                               seperator[i][k][n], seperator[i][k][n+1])) >= 0)
    372             mat[mm][index2 + found - seperator [i][k][n]]= 1;
    373           index2 += seperator [i][k][n+1]-seperator[i][k][n];
    374         }
    375       }
    376 
    377       index2= 0;
    378       for (i= j - 1; i >= 0; i--)
    379       {
    380         if (i == index1)
    381           continue;
    382         oneCount= 0;
    383         for (mm= 0; mm < seperator [i][k][n + 1] - seperator [i][k][n]; mm++)
    384         {
    385           for (m= 0; m < maxTerms; m++)
    386           {
    387             if (mat[m][mm+index2] == 1)
    388               oneCount++;
    389           }
    390         }
    391         if (oneCount == seperator [i][k][n+1]-seperator[i][k][n] - 1)
    392         {
    393           for (mm= 0; mm < seperator [i][k][n+1]-seperator[i][k][n]; mm++)
    394           {
    395             oneCount= 0;
    396             for (m= 0; m < maxTerms; m++)
    397               if (mat[m][mm+index2] == 1)
    398                 oneCount++;
    399             if (oneCount > 0)
    400               continue;
    401             for (m= 0; m < maxTerms; m++)
    402             {
    403               oneCount= 0;
    404               for (mmm= 0; mmm < seperator[i][k][n+1]-seperator[i][k][n]; mmm++)
    405               {
    406                 if (mat[m][mmm+index2] == 1)
    407                   oneCount++;
    408               }
    409               if (oneCount > 0)
    410                 continue;
    411               mat[m][mm+index2]= 1;
    412             }
    413           }
    414         }
    415         index2 += seperator [i][k][n+1] - seperator [i][k][n];
    416       }
    417 
    418       //read off solution
    419       mm= 0;
    420       for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
    421       {
    422         tmp1= storeFactors [index1][0][k][m];
    423         index2= 0;
    424         for (i= j - 1; i > -1; i--)
    425         {
    426           if (i == index1)
    427             continue;
    428           for (mmm= 0; mmm < seperator [i][k][n+1]-seperator[i][k][n]; mmm++)
    429             if (mat[mm][mmm+index2] == 1)
    430               tmp1= patch (tmp1, storeFactors[i][0][k][seperator[i][k][n]+mmm],
    431                            eval[i]);
    432           index2 += seperator [i][k][n+1]-seperator[i][k][n];
    433         }
    434         factor += tmp1;
    435       }
    436 
    437       for (m= 0; m < maxTerms; m++)
    438         delete [] mat [m];
    439       delete [] mat;
    440     }
    441 
    442     if (fdivides (factor, B, quot))
    443     {
    444       result.append (factor);
    445       B= quot;
    446       if (result.length() == biFactors.length() - 1)
    447       {
    448         result.append (quot);
    449         break;
    450       }
    451     }
    452   }
    453 
    454   //delete
    455   for (i= 0; i < j; i++)
    456   {
    457     delete [] storeFactors [i] [0];
    458     delete [] storeFactors [i] [1];
    459     delete [] storeFactors [i];
    460     for (k= 0; k < minFactorsLength; k++)
    461       delete [] seperator [i][k];
    462     delete [] seperator [i];
    463   }
    464   delete [] seperator;
    465   delete [] storeFactors;
    466   //
    467 
    468   return result;
    469 }
    47042
    47143static inline
    … …  
    15911163}
    15921164
    1593 CFList
    1594 leadingCoeffReconstruction (const CanonicalForm& F, const CFList& factors,
    1595                             const CFList& M)
    1596 {
    1597   CanonicalForm quot, buf= F;
    1598   CanonicalForm LCBuf= LC (buf, 1);
    1599 
    1600   CFList result;
    1601 
    1602   CanonicalForm tmp;
    1603   for (CFListIterator i= factors; i.hasItem(); i++)
    1604   {
    1605     tmp= i.getItem();
    1606     tmp= mulMod (tmp, LCBuf, M);
    1607     tmp= tmp/content (tmp, 1);
    1608     if (fdivides (tmp, buf, quot))
    1609     {
    1610       buf= quot;
    1611       result.append (tmp);
    1612       LCBuf= LC (buf, 1);
    1613     }
    1614     else //no one-to-one correspondence
    1615       return CFList();
    1616   }
    1617 
    1618   return result;
    1619 }
    1620 
    16211165void
    16221166gcdFreeBasis (CFFList& factors1, CFFList& factors2)

  • factory/facFqFactorize.h

    r6deedd r91788c0  
    264264                    );
    265265
     266/// recombination of bivariate factors @a factors1 s. t. the result evaluated
     267/// at @a evalPoint coincides with @a factors2
     268CFList
     269recombination (const CFList& factors1,        ///<[in] list of bivariate factors
     270               const CFList& factors2,        ///<[in] list univariate factors
     271               int s,                         ///<[in] algorithm starts checking
     272                                              ///<  subsets of size s
     273               int thres,                     ///<[in] threshold for the size of
     274                                              ///<  subsets which are checked
     275               const CanonicalForm& evalPoint,///<[in] evaluation point
     276               const Variable& x              ///<[in] second variable of
     277                                              ///< bivariate factors
     278              );
     279
    266280/// Lift bound adaption. Essentially an early factor detection but only the lift
    267281/// bound is adapted.
Note: See TracChangeset for help on using the changeset viewer.