Changeset 806c18 in git
- Timestamp:
- Nov 15, 2010, 4:34:57 PM (13 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
- Children:
- 7c3bca08c96331a56864c1d35b8c2e8ff2e0be89
- Parents:
- c840d97af622b4e4da8761738b540e21144f716b
- Location:
- factory
- Files:
-
- 69 edited
Legend:
- Unmodified
- Added
- Removed
-
factory/DegreePattern.cc
rc840d9 r806c18 1 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 2 * Computer Algebra System SINGULAR 3 3 \*****************************************************************************/ 4 4 /** @file DegreePattern.cc 5 * 6 * This file provides functions for manipulating DegreePatterns 5 * 6 * This file provides functions for manipulating DegreePatterns 7 7 * 8 8 * @author Martin Lee … … 15 15 #include <config.h> 16 16 17 #include "DegreePattern.h" 17 #include "DegreePattern.h" 18 18 #include "cf_iter.h" 19 19 #include "ftmpl_functions.h" … … 21 21 22 22 23 DegreePattern::DegreePattern (const CFList& l) 24 { 23 DegreePattern::DegreePattern (const CFList& l) 24 { 25 25 m_data = NULL; 26 26 … … 29 29 else 30 30 { 31 31 32 32 Variable x= Variable (1); 33 33 int p= getCharacteristic(); … … 42 42 int j= 0; 43 43 for (CFIterator i= buf; i.hasTerms(); i++, j++) 44 ; 44 ; 45 45 46 46 ASSERT ( j > 1, "j > 1 expected" ); 47 47 48 48 m_data = new Pattern( j - 1 ); 49 49 … … 62 62 void DegreePattern::intersect (const DegreePattern& degPat) 63 63 { 64 if (degPat.getLength() < getLength()) 64 if (degPat.getLength() < getLength()) 65 65 { 66 66 DegreePattern bufDeg= *this; … … 72 72 int length= tmin (getLength(), degPat.getLength()); 73 73 int* buf= new int [length]; 74 for (int i= 0; i < length; i++) 74 for (int i= 0; i < length; i++) 75 75 { 76 if (degPat.find ((*this)[i])) 77 { 76 if (degPat.find ((*this)[i])) 77 { 78 78 buf[i]= (*this)[i]; 79 79 count++; … … 86 86 init (count); 87 87 count= 0; 88 for (int i= 0; i < length; i++) 88 for (int i= 0; i < length; i++) 89 89 { 90 if (buf[i] != -1) 90 if (buf[i] != -1) 91 91 { 92 92 (*this) [count]= buf[i]; … … 96 96 delete[] buf; 97 97 } 98 98 99 99 void DegreePattern::refine () 100 100 { … … 107 107 for (int i= 0; i < getLength(); i++) 108 108 buf[i]= -1; 109 for (int i= 1; i < getLength(); i++) 109 for (int i= 1; i < getLength(); i++) 110 110 { 111 111 pos= (*this).find (d - (*this)[i]); 112 if (pos) 113 { 112 if (pos) 113 { 114 114 buf[i]= (*this)[i]; 115 115 count++; … … 127 127 init (count); 128 128 count= 0; 129 for (int i= 0; i < length; i++) 129 for (int i= 0; i < length; i++) 130 130 { 131 if (buf[i] != -1) 131 if (buf[i] != -1) 132 132 { 133 133 (*this)[count]= buf[i]; -
factory/DegreePattern.h
rc840d9 r806c18 1 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 2 * Computer Algebra System SINGULAR 3 3 \*****************************************************************************/ 4 4 /** @file DegreePattern.h 5 * 6 * This file provides a class to handle degree patterns. 5 * 6 * This file provides a class to handle degree patterns. 7 7 * 8 8 * @author Martin Lee … … 25 25 26 26 /** @class DegreePattern DegreePattern.h "factory/DegreePattern.h" 27 * 27 * 28 28 * DegreePattern provides a functionality to create, intersect and refine 29 29 * degree patterns. 30 30 * 31 * 31 * 32 32 */ 33 33 class DegreePattern … … 52 52 ASSERT ( m_data != NULL, "non-null pointer expected"); 53 53 ASSERT ( m_data->m_refCounter == 0, "ref count of 0 expected"); 54 if( m_data->m_pattern != NULL ) 54 if( m_data->m_pattern != NULL ) 55 55 delete[] m_data->m_pattern; 56 56 m_data->m_pattern = NULL; 57 57 58 58 delete m_data; 59 59 m_data = NULL; … … 67 67 if( (--m_data->m_refCounter) < 1 ) 68 68 release(); 69 69 70 70 m_data = new Pattern(n); 71 71 } … … 94 94 /// operator [] 95 95 /// 96 /// @return @a operator[] returns the element at @a index 96 /// @return @a operator[] returns the element at @a index 97 97 inline int operator[] (const int index ///< [in] some int >= 0, < getLength() 98 98 ) const … … 106 106 /// operator [] 107 107 /// 108 /// @return @a operator[] sets the element at @a index 108 /// @return @a operator[] sets the element at @a index 109 109 inline int& operator[] (const int index ///< [in] some int >= 0, < getLength() 110 110 ) … … 121 121 /// copy constructor 122 122 DegreePattern (const DegreePattern& degPat ///< [in] some degree pattern 123 ): m_data( degPat.m_data ) 123 ): m_data( degPat.m_data ) 124 124 { 125 125 ASSERT( degPat.m_data != NULL, "non-null pointer expected" ); … … 130 130 DegreePattern (const CFList& l ///< [in] some list of (univariate) polys 131 131 ); 132 132 133 133 /// assignment 134 134 DegreePattern& operator= (const DegreePattern& degPat ///< [in] some degree … … 148 148 149 149 /// destructor 150 ~DegreePattern () 150 ~DegreePattern () 151 151 { 152 152 ASSERT( m_data != NULL, "non-null pointer expected" ); 153 if( (--m_data->m_refCounter) < 1 ) 153 if( (--m_data->m_refCounter) < 1 ) 154 154 release(); 155 155 } 156 156 157 /// find an element @a x 157 /// find an element @a x 158 158 /// 159 /// @return @a find returns the index + 1 of @a x, if @a x is an element of 159 /// @return @a find returns the index + 1 of @a x, if @a x is an element of 160 160 /// the degree pattern, 0 otherwise 161 161 int find (const int x ///< [in] some int … … 167 167 return 0; 168 168 }; 169 170 /// intersect two degree patterns 169 170 /// intersect two degree patterns 171 171 void intersect (const DegreePattern& degPat ///< [in] some degree pattern 172 172 ); 173 173 /// Refine a degree pattern. Assumes that (*this)[0]:= @a d is the degree 174 174 /// of the poly to be factored. Now for every other entry @a a there should be 175 /// some entry @a b such that @a a+b= d. Elements which do not satisfy this 175 /// some entry @a b such that @a a+b= d. Elements which do not satisfy this 176 176 /// relation are removed. 177 177 void refine (); -
factory/ExtensionInfo.cc
rc840d9 r806c18 1 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 2 * Computer Algebra System SINGULAR 3 3 \*****************************************************************************/ 4 4 /** @file ExtensionInfo.cc 5 * 6 * This file provides member functions for ExtensionInfo 5 * 6 * This file provides member functions for ExtensionInfo 7 7 * 8 8 * @author Martin Lee … … 28 28 } 29 29 30 ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 31 const CanonicalForm& gamma, const CanonicalForm& 30 ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 31 const CanonicalForm& gamma, const CanonicalForm& 32 32 delta, const int nGFDegree, const char cGFName, 33 33 const bool extension) … … 42 42 } 43 43 44 ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 45 const CanonicalForm& gamma, const CanonicalForm& 44 ExtensionInfo::ExtensionInfo (const Variable& alpha, const Variable& beta, 45 const CanonicalForm& gamma, const CanonicalForm& 46 46 delta) 47 47 { … … 66 66 } 67 67 68 #if 0 68 69 ExtensionInfo::ExtensionInfo (const Variable& alpha) 69 70 { … … 76 77 m_extension= true; 77 78 } 79 #endif 78 80 79 81 ExtensionInfo::ExtensionInfo (const int nGFDegree, const char cGFName, const -
factory/ExtensionInfo.h
rc840d9 r806c18 1 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 2 * Computer Algebra System SINGULAR 3 3 \*****************************************************************************/ 4 4 /** @file ExtensionInfo.h 5 * 5 * 6 6 * This file provides a class to store information about finite fields and 7 * extensions thereof. 8 * 7 * extensions thereof. 8 * 9 9 * 10 10 * @author Martin Lee … … 23 23 24 24 /** @class ExtensionInfo ExtensionInfo.h "factory/ExtensionInfo.h" 25 * ExtensionInfo contains information about extension. 26 * If @a m_extension is true we are in an extension of some initial field. 27 * If the initial field is \f$ F_p \f$ and we pass to \f$ F_p (\alpha) \f$ 28 * then @a m_alpha is an algebraic variable, @a m_beta= Variable(1), 25 * ExtensionInfo contains information about extension. 26 * If @a m_extension is true we are in an extension of some initial field. 27 * If the initial field is \f$ F_p \f$ and we pass to \f$ F_p (\alpha) \f$ 28 * then @a m_alpha is an algebraic variable, @a m_beta= Variable(1), 29 29 * @a m_gamma= @a m_delta= 1, @a m_GFDegree= 0, @a m_GFName= 'Z'. If we pass 30 * to some GF (p^k) then @a m_alpha= Variable (1), @a m_beta= Variable(1), 31 * @a m_gamma= @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'. 32 * @n If the initial field is \f$ F_p (\epsilon) \f$, then @a m_beta= 33 * \f$ \epsilon \f$, @a m_alpha an algebraic variable defining an extension of 34 * \f$ F_p (\epsilon) \f$, @a m_gamma is a primitive element of 35 * \f$ F_p (\alpha) \f$, @a m_delta is a primitive element of 36 * \f$ F_p (\beta) \f$, @a m_GFDegree= 0, @a m_GFName= 'Z'. 37 * @n If the initial field is GF(p^k), then @a m_alpha= Variable (1), 38 * @a m_beta= Variable (1), @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree()= k, 39 * @a m_GFName= gf_name of the initial field. 40 * @n If @a m_extension is false and the current field is \f$ F_p \f$ then 41 * @a m_alpha= Variable (1), @a m_beta= Variable (1), @a m_gamma= 1, 42 * @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'. 43 * @n If the current field is \f$ F_p (\alpha) \f$ then 44 * @a m_alpha is some algebraic variable, @a m_beta= Variable (1), 45 * @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree= 0, @a m_GFName= 'Z'. 46 * @n If the current field is GF then @a m_alpha= Variable (1), 47 * @a m_beta= Variable (1), @a m_gamma= 1, @a m_delta= 1, 30 * to some GF (p^k) then @a m_alpha= Variable (1), @a m_beta= Variable(1), 31 * @a m_gamma= @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'. 32 * @n If the initial field is \f$ F_p (\epsilon) \f$, then @a m_beta= 33 * \f$ \epsilon \f$, @a m_alpha an algebraic variable defining an extension of 34 * \f$ F_p (\epsilon) \f$, @a m_gamma is a primitive element of 35 * \f$ F_p (\alpha) \f$, @a m_delta is a primitive element of 36 * \f$ F_p (\beta) \f$, @a m_GFDegree= 0, @a m_GFName= 'Z'. 37 * @n If the initial field is GF(p^k), then @a m_alpha= Variable (1), 38 * @a m_beta= Variable (1), @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree()= k, 39 * @a m_GFName= gf_name of the initial field. 40 * @n If @a m_extension is false and the current field is \f$ F_p \f$ then 41 * @a m_alpha= Variable (1), @a m_beta= Variable (1), @a m_gamma= 1, 42 * @a m_delta= 1, @a m_GFDegree= 1, @a m_GFName= 'Z'. 43 * @n If the current field is \f$ F_p (\alpha) \f$ then 44 * @a m_alpha is some algebraic variable, @a m_beta= Variable (1), 45 * @a m_gamma= 1, @a m_delta= 1, @a m_GFDegree= 0, @a m_GFName= 'Z'. 46 * @n If the current field is GF then @a m_alpha= Variable (1), 47 * @a m_beta= Variable (1), @a m_gamma= 1, @a m_delta= 1, 48 48 * @a m_GFDegree= getGFDegree(), @a m_GFName= gf_name. 49 49 * 50 * @sa facFqBivar.h, facFqFactorize.h 50 * @sa facFqBivar.h, facFqFactorize.h 51 51 */ 52 52 class ExtensionInfo 53 53 { 54 54 private: 55 /// an algebraic variable or Variable (1) 55 /// an algebraic variable or Variable (1) 56 56 Variable m_alpha; 57 57 /// an algebraic variable or Variable (1) 58 Variable m_beta; 58 Variable m_beta; 59 59 /// a primitive element of \f$ F_p (\alpha) \f$ or 1 60 60 CanonicalForm m_gamma; … … 69 69 public: 70 70 /// \f$ F_p \f$ as initial field, if @a extension is true we are in some GF 71 ExtensionInfo (const bool extension ///< [in] some bool 71 ExtensionInfo (const bool extension ///< [in] some bool 72 72 ); 73 /// Construct an @a ExtensionInfo 74 ExtensionInfo (const Variable& alpha, ///< [in] some algebraic variable 73 /// Construct an @a ExtensionInfo 74 ExtensionInfo (const Variable& alpha, ///< [in] some algebraic variable 75 75 const Variable& beta, ///< [in] some algebraic variable 76 76 const CanonicalForm& gamma, ///< [in] some primitive element … … 85 85 /// \f$ F_p (\beta) \f$ as initial field and switch to an extension given by 86 86 /// @a alpha, needs primitive elements @a gamma and @a delta for maps 87 /// between \f$ F_p (\alpha) \subset F_p (\beta) \f$ 87 /// between \f$ F_p (\alpha) \subset F_p (\beta) \f$ 88 88 ExtensionInfo (const Variable& alpha, ///< [in] some algebraic variable 89 89 const Variable& beta, ///< [in] some algebraic variable … … 94 94 ); 95 95 /// \f$ F_p (\alpha) \f$ as initial field, if @a extension is false. 96 /// Else initial field is \f$ F_p \f$ 96 /// Else initial field is \f$ F_p \f$ 97 97 ExtensionInfo (const Variable& alpha, ///< [in] some algebraic variable 98 98 const bool extension ///< [in] some bool 99 ); 99 ); 100 100 101 ExtensionInfo (const Variable& alpha ///< [in] some algebraic variable 101 //ExtensionInfo (const Variable& alpha ///< [in] some algebraic variable 102 // ); 103 104 /// GF as initial field 105 ExtensionInfo (const int nGFDegree, ///< [in] GF degree of initial field 106 const char cGFName, ///< [in] name of GF variable 107 const bool extension ///< [in] some bool 102 108 ); 103 104 /// GF as initial field105 ExtensionInfo (const int nGFDegree, ///< [in] GF degree of initial field106 const char cGFName, ///< [in] name of GF variable107 const bool extension ///< [in] some bool108 );109 109 110 110 /// getter -
factory/NTLconvert.cc
rc840d9 r806c18 47 47 #define Free(A,L) free(A) 48 48 #endif 49 49 50 50 void out_cf(const char *s1,const CanonicalForm &f,const char *s2); 51 51 52 52 53 53 int fac_NTL_char=-1; // the current characterstic for NTL calls … … 625 625 626 626 //reverse the list to obtain the correct string 627 //NTL_SNS 627 //NTL_SNS 628 628 for (int i=l-1;i>=0;i--) // l ist the position of \0 629 629 { … … 1069 1069 } 1070 1070 1071 CanonicalForm convertNTLzz_pEX2CF (zz_pEX f, Variable x, Variable alpha) 1071 CanonicalForm convertNTLzz_pEX2CF (zz_pEX f, Variable x, Variable alpha) 1072 1072 { 1073 1073 CanonicalForm bigone; … … 1076 1076 bigone= 0; 1077 1077 bigone.mapinto(); 1078 for (int j=0;j<deg(f)+1;j++) 1078 for (int j=0;j<deg(f)+1;j++) 1079 1079 { 1080 1080 if (coeff(f,j)!=0) -
factory/abs_fac.cc
rc840d9 r806c18 36 36 37 37 if( F.isZero() ) return 0; 38 if( F.inBaseDomain() ) return F; 39 40 if( level(F) < 0 ) return 1; 38 if( F.inBaseDomain() ) return F; 39 40 if( level(F) < 0 ) return 1; 41 41 42 42 r = LC( F, x); … … 47 47 { 48 48 t = LC(g, x); 49 if( t == 1 || t == -1 ) return 1; 49 if( t == 1 || t == -1 ) return 1; 50 50 r = MYGCD( r, t); 51 if( r == 1 ) return 1; 52 g = g - power(x,degree(g,x))*t; 51 if( r == 1 ) return 1; 52 g = g - power(x,degree(g,x))*t; 53 53 } 54 54 return r; … … 613 613 CFFListIterator J=Result; 614 614 for ( ; J.hasItem(); J++) 615 { 615 { 616 616 Factor_Norm = J.getItem().factor(); 617 617 Factor_Norm = Factor_Norm(x+k*l,x); // die Störungen werden rückgänig gemacht … … 861 861 F = A[i]; 862 862 863 if( degree(F,x) == 0 ) 864 if( c.level() < 0 ) return 1; else return c; 863 if( degree(F,x) == 0 ) 864 if( c.level() < 0 ) return 1; else return c; 865 865 866 866 F = gamma*F/LC(F, x); … … 908 908 if( g != g(0,i) ) 909 909 L.append(i); 910 910 911 911 int nvg = L.length(); 912 912 913 913 914 914 915 915 if( f.level() < 0 && g.level() < 0 ) { ; 916 916 return 1; } 917 917 918 918 CFArray A; 919 919 920 920 i=0; 921 921 int r; … … 940 940 if( nvf <= 1 && nvg <=1 ) 941 941 { 942 gamma = MyGCDmod( LC(F,x), LC(G,x) ); 942 gamma = MyGCDmod( LC(F,x), LC(G,x) ); 943 943 c = MyGCDmod( Cf, Cg ); 944 944 } … … 956 956 F = A[i]; 957 957 958 if( degree(F,x) == 0 ) 959 if( c.level() < 0 ) return 1; else return c; 958 if( degree(F,x) == 0 ) 959 if( c.level() < 0 ) return 1; else return c; 960 960 961 961 F = gamma*F/LC(F, x); … … 1007 1007 { 1008 1008 CanonicalForm qi = MYGCD( temp2, temp4); 1009 if( qi != 1 ) L.append( CFFactor( qi, i ) ); 1009 if( qi != 1 ) L.append( CFFactor( qi, i ) ); 1010 1010 i++; 1011 1011 temp2 = temp2/qi; … … 1037 1037 for ( CFIterator I = F; I.hasTerms(); ++I ) fillVarsRec( I.coeff(), vrs ); 1038 1038 1039 N.append( CFFactor(F,1) ); 1039 N.append( CFFactor(F,1) ); 1040 1040 1041 1041 int i = n+1; … … 1043 1043 while( i >= 0 ) 1044 1044 { 1045 b = 0; 1045 b = 0; 1046 1046 1047 1047 if( i == 0 ){ v = mvar(F); b=1 ;} 1048 1048 else 1049 if( vrs[i] != 0 ){ b=1; v= Variable(i);} 1050 if( vrs[i] == 0 ) i--; 1049 if( vrs[i] != 0 ){ b=1; v= Variable(i);} 1050 if( vrs[i] == 0 ) i--; 1051 1051 1052 1052 if( b ) 1053 { 1053 { 1054 1054 for( CFFListIterator J = L; J.hasItem(); J++ ) 1055 1055 { … … 1065 1065 } 1066 1066 if( N.length() == L.length() ) i -= 1; 1067 L=N; 1067 L=N; 1068 1068 } 1069 1069 } -
factory/algext.cc
rc840d9 r806c18 569 569 other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg) 570 570 other = leadDeg(Dp,other); 571 571 572 572 if(isEqual(bound, other, 1, mv)) // equal 573 573 { -
factory/bifac.cc
rc840d9 r806c18 4 4 #include "bifacConfig.h" 5 5 6 #define BIFAC_BASIS_OF_G_CHECK 6 #define BIFAC_BASIS_OF_G_CHECK 1 7 7 void Reduce( bool ); 8 8 CanonicalForm Bigcd( const CanonicalForm& f, const CanonicalForm& g); … … 28 28 29 29 //================================================== 30 class PolyVector 30 class PolyVector 31 31 //================================================== 32 32 { … … 53 53 int correction = 1; // univariate polynomials 54 54 if( n==0) correction = n+1; 55 55 56 56 value = new CanonicalForm[m*(n+1)+n+1]; 57 57 for(int i=0; i<=m*(n+1)+n; i++) value[i]=0; 58 58 59 60 for ( CFIterator i = f; i.hasTerms(); i++ ) { 61 62 if( i.coeff().mvar().level()< 0 ){ 63 64 else{ 65 59 60 for ( CFIterator i = f; i.hasTerms(); i++ ) { 61 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ){ 62 if( i.coeff().mvar().level()< 0 ){ 63 value[ 0*(n+1) + i.exp()*correction ] = j.coeff();} 64 else{ 65 value[ j.exp()*(n+1) + i.exp()*correction ] = j.coeff();}}} 66 66 } 67 67 } … … 84 84 s << "["; 85 85 for (int j=0;j<=V.n;j++) 86 86 s << V.value[i*(V.n+1)+j] << ", "; 87 87 s << "]\n"; 88 88 } … … 120 120 //--<>--------------------------------- 121 121 { 122 } 122 } 123 123 124 124 … … 137 137 138 138 // // === Datei löschen === 139 140 // ofstream* aus = new ofstream(name, ios::out); 139 140 // ofstream* aus = new ofstream(name, ios::out); 141 141 // delete aus; 142 142 … … 144 144 // // === Jetzt immer nur anhängen === 145 145 146 // aus = new ofstream(name, ios::app); 146 // aus = new ofstream(name, ios::app); 147 147 // *aus << "// Zeilen Spalten\n" 148 148 // << "// x-Koord. y-Koord. Wert\n"; … … 157 157 // *aus << i << " " << j << " " << M(i+1,j+1) << endl;; 158 158 // delete aus; 159 // } 159 // } 160 160 161 161 //======================================================= … … 163 163 //======================================================= 164 164 { 165 165 ; 166 166 } 167 167 … … 174 174 175 175 for ( CFIterator i = f; i.hasTerms(); i++ ) 176 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 176 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 177 177 z++; 178 178 return( z ); 179 179 } 180 180 181 181 //======================================================= 182 182 void BIFAC::biGanzMachen( CanonicalForm & f ) … … 188 188 189 189 for ( CFIterator i = f; i.hasTerms(); i++ ) 190 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 190 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 191 191 { 192 192 if( !init ) 193 { 194 195 193 { 194 ggT = j.coeff(); 195 init = true; 196 196 } 197 else 198 197 else 198 ggT = gcd(j.coeff(), ggT); 199 199 } 200 200 f /= ggT; … … 203 203 204 204 //======================================================= 205 void BIFAC::biNormieren( CanonicalForm & f ) 205 void BIFAC::biNormieren( CanonicalForm & f ) 206 206 //======================================================= 207 207 { … … 209 209 { 210 210 for ( CFIterator i = f; i.hasTerms(); i++ ) 211 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 212 213 214 215 216 211 for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ) 212 if( j.coeff().den() != 1 ) 213 { 214 f *= j.coeff().den(); 215 biNormieren( f ); 216 } 217 217 biGanzMachen( f ); 218 } 218 } 219 219 else 220 220 { … … 242 242 for(i=0; i<=m-1; i++) 243 243 for(j=0; j<=n; j++) 244 244 g += A(k, i*(n+1)+j+1)* power(x,i) * power(y,j); 245 245 Lg.append(g); 246 246 } … … 255 255 h=0; 256 256 for(i=0; i<=m; i++) 257 258 257 for(j=0; j<n; j++) 258 h += A(k, i*n+j+1 +m*(n+1))* power(x,i) * power(y,j); 259 259 Lh.append(h); 260 260 } 261 261 262 262 // === Is the solution correct? === 263 263 CFListIterator itg=Lg; … … 268 268 h = ith.getItem(); 269 269 ff = f*(deriv(g,y)-deriv(h,x)) +h*deriv(f,x) -g*deriv(f,y); 270 if( !ff.isZero()) { 270 if( !ff.isZero()) { 271 271 #ifndef NOSTREAMIO 272 AUSGABE_ERR("* Falsche Polynome!"); 273 exit (1); 272 AUSGABE_ERR("* Falsche Polynome!"); 273 exit (1); 274 274 #else 275 275 printf("wrong polys\n"); … … 277 277 #endif 278 278 } 279 } 279 } 280 280 } 281 281 /////////// END VALIDATION //////////////////////////////////// … … 302 302 int n = degree(f,y); 303 303 int r,s, ii,jj; 304 304 305 305 306 306 // ======= Creation of the system of linear equations for G ============= … … 309 309 310 310 CFMatrix M(rows, columns); // Remember: The first index is (1,1) -- not (0,0)! 311 311 312 312 for ( CFIterator i = f; i.hasTerms(); i++ ) // All coeffizients of y 313 313 { … … 319 319 // Now we regard g_{ii,jj) 320 320 for( ii=0; ii<m; ii++) 321 322 323 324 325 321 for( jj=0; jj<=n; jj++) 322 { 323 if( s>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) += -j.coeff() * s; 324 if( jj>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) += j.coeff() * jj; 325 } 326 326 327 327 // Now we regard h_{ii,jj} 328 328 for( ii=0; ii<=m; ii++) 329 330 331 332 333 329 for( jj=0; jj<n; jj++) 330 { 331 if( r>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n)+jj +m*(n+1) +1) += j.coeff() * r; 332 if( ii>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n) +jj +m*(n+1) +1) += -j.coeff() * ii; 333 } 334 334 } 335 335 } … … 364 364 tmp =0; 365 365 for(int j=1; j<=columns; j++) 366 366 tmp += M(i,j) * basis(k,j); 367 367 if( tmp!= 0) { 368 368 exit(17); 369 369 } 370 370 } … … 374 374 375 375 //======================================================= 376 // Compute a r x r - matrix A=(a_ij) for 376 // Compute a r x r - matrix A=(a_ij) for 377 377 // gg_i = SUM a_ij * g_j * f_x (mod f) 378 // Return a list consisting of 378 // Return a list consisting of 379 379 // r x (r+1) Matrix A 380 // the last columns contains only the indices of the 380 // the last columns contains only the indices of the 381 381 // first r linear independent lines 382 382 // REMARK: this is used by BIFAC::createEg but NOT by createEgUni!! … … 390 390 int r = G.length(); // number of factors 391 391 392 LGS L(r,r,true); 393 // LGS L(r,r); 392 LGS L(r,r,true); 393 // LGS L(r,r); 394 394 CFMatrix Z(1,r); 395 395 CFMatrix A(r,r+2); // the last two column contain the bi-degree … … 402 402 CanonicalForm q; 403 403 404 for( CFListIterator it=G; it.hasItem(); it++, i++){ 404 for( CFListIterator it=G; it.hasItem(); it++, i++){ 405 405 406 406 gifx[i].init( (it.getItem()*fx)%f ); … … 410 410 411 411 e=1; // row number of A 412 n=0; // 412 n=0; // 413 413 m=0; // 414 414 while (L.rank() != r ) … … 461 461 g += rand_coeff1[i] * it.getItem(); 462 462 } 463 463 464 464 delete[] rand_coeff1; 465 465 … … 483 483 484 484 bool suitable1 = false; // Is Eg by chance unsuitable? 485 bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? 485 bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? 486 486 487 487 // === (0) Preparation === … … 511 511 { 512 512 513 suitable2 = false; 513 suitable2 = false; 514 514 // === (1) Creating g === 515 515 while ( !suitable2 ) … … 519 519 // for( CFListIterator it=G; it.hasItem(); it++, i++) 520 520 // { 521 // 522 // 523 // 521 // gi[i] = it.getItem(); 522 // rand_coeff[i] = RANDOM.generate().intval(); 523 // g += rand_coeff[i] * it.getItem(); 524 524 // } 525 525 g = create_g( G ); 526 526 527 527 // === (2) Computing g_i * g === 528 529 for(i=0; i<r; i++){ 530 531 528 // 529 for(i=0; i<r; i++){ 530 531 ggi[i] = (g*gi[i])%f; // seite 10 532 532 } 533 533 534 534 // === Check if all polynomials are <> 0 === 535 535 suitable2 = true; // It should be fine, but ... 536 if( g.isZero() ) 537 536 if( g.isZero() ) 537 suitable2 = false; 538 538 // else 539 // 540 // 541 // 539 // for(i=0; i<r; i++) 540 // if( ggi[i].isZero() ) 541 // suitable2 = false; 542 542 543 543 } // end of Žwhile ( !suitable2 )Ž 544 544 545 545 // === (3) Computing Eg(x) === 546 546 547 547 for(i=0;i<r;i++) // Get Polynomials as vectors 548 548 v_ggi[i].init(ggi[i]); 549 549 550 550 // Matrix A 551 for(i=1; i<=r; i++) 551 for(i=1; i<=r; i++) 552 552 for( j=1; j<=r; j++) 553 553 { 554 555 556 554 A(i,j) = 0; 555 for( e=1; e<=r; e++) 556 { 557 557 558 558 … … 562 562 // 563 563 564 564 } 565 565 } 566 566 … … 605 605 CanonicalForm ff, ffx,g, gg, Eg; 606 606 607 607 608 608 bool suitable1 = false; // Is Eg unsuitable? 609 bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? 609 bool suitable2 = false; // Is on of g*g_i or g_i*f_x zero? 610 610 bool suitable3 = false; // Is 'konst' unsuitable? 611 611 … … 627 627 CFMatrix Z(1,r); // `Vector` for data transportation 628 628 CFMatrix AA(m,r); // but first we generate AA. 629 CFMatrix AI(r,r+1); // 630 LGS L(r,r,true); 629 CFMatrix AI(r,r+1); // 630 LGS L(r,r,true); 631 631 IntRandom RANDOM(S); 632 632 … … 642 642 ff = f(konst,'y'); 643 643 ffx = fx(konst,'y'); 644 644 645 645 if( gcd(ff, ffx) == 1) 646 646 suitable3 = true; … … 649 649 konst *= -1; 650 650 if( konst >= 0 ) 651 651 konst++; 652 652 } 653 653 } … … 660 660 for( CFListIterator it=G; it.hasItem(); it++, i++) 661 661 { 662 gi[i] = it.getItem()(konst,'y'); 662 gi[i] = it.getItem()(konst,'y'); 663 663 } 664 664 … … 666 666 // = (3) Compute the matrices 'AA' and 'AI' = 667 667 // =============================================== 668 668 669 669 670 670 for( i=0; i<r; i++) // First store all coeffizients in AA. … … 673 673 //biNormieren(ggi[i]); 674 674 for ( CFIterator j = ggi[i]; j.hasTerms(); j++ ) 675 AA( j.exp()+1, i+1) = j.coeff(); 675 AA( j.exp()+1, i+1) = j.coeff(); 676 676 } 677 677 … … 684 684 ASSERT( i<=m, "Too few linear independent rows!"); 685 685 686 for (k=1; k<=r; k++) 686 for (k=1; k<=r; k++) 687 687 Z(1,k) = AA(i,k); 688 688 if( L.new_row(Z,0) ) // linear independent row? … … 701 701 // = (4) Big loop to find a suitable 'Eg(x) = 702 702 // ============================================== 703 703 704 704 while ( !suitable1 ) // Is Eg(x) suitable? -> Check at the end of this procedure! 705 705 { … … 711 711 // rand_coeff[0] = 0; 712 712 // rand_coeff[1] = 4; 713 713 714 714 715 715 while ( !suitable2 ) … … 720 720 for( CFListIterator it=G; it.hasItem(); it++, i++) 721 721 { 722 723 722 rand_coeff[i] = RANDOM.generate().intval(); 723 g += rand_coeff[i] * it.getItem(); 724 724 } 725 725 gg = g(konst,'y'); // univariate! … … 728 728 // === (ii) Check if all polynomials are <> 0 === 729 729 suitable2 = true; // It should be fine, but ... 730 if( gg.isZero() ) 731 730 if( gg.isZero() ) 731 suitable2 = false; 732 732 // else 733 // 734 // 735 // 733 // for(i=0; i<r; i++) 734 // if( ggi[i].isZero() ) 735 // suitable2 = false; 736 736 } // end of Žwhile ( !suitable2 )Ž 737 737 738 738 // createRg(g,f); 739 739 … … 741 741 // = (b) Compute matrix 'A' = 742 742 // =============================================== 743 for(i=1; i<=r; i++) 744 { 745 for( ii=1; ii<=m; ii++) 746 743 for(i=1; i<=r; i++) 744 { 745 for( ii=1; ii<=m; ii++) 746 AA (ii,1) = 0; // !! Redefinition of AA !! 747 747 for ( CFIterator j = ggi[i-1]; j.hasTerms(); j++ ) 748 AA( j.exp()+1, 1) = j.coeff(); 748 AA( j.exp()+1, 1) = j.coeff(); 749 749 750 750 for( ii=1; ii<=r; ii++) 751 751 { 752 753 754 A(i,ii) += ( AI(ii,k ) * AA( AI(k, r+1 ).intval(),1) ); 752 A(i,ii) = 0; 753 for( k=1; k<=r; k++) 754 A(i,ii) += ( AI(ii,k ) * AA( AI(k, r+1 ).intval(),1) ); 755 755 } 756 756 } … … 761 761 // = (c) Compute Eg(x) and check it = 762 762 // =============================================== 763 763 764 764 Eg = determinant(A,r); 765 765 if( gcd(Eg, deriv(Eg,x)) == 1 ) … … 769 769 } // end of Žwhile ( !suitable1 )Ž 770 770 771 771 772 772 // ============================================== 773 773 // = (5) Prepare for leaving = … … 777 777 delete[] ggi; 778 778 delete[] rand_coeff; 779 779 780 780 CFList LL; 781 781 LL.append(Eg); … … 805 805 // =============================================== 806 806 807 CanonicalForm alpha=1; 808 809 while( resultant( f, fx, x)(alpha) == 0 ) 810 { 811 //while( resultant( f, fx, x)(alpha).inCoeffDomain() != true ) 807 CanonicalForm alpha=1; 808 809 while( resultant( f, fx, x)(alpha) == 0 ) 810 { 811 //while( resultant( f, fx, x)(alpha).inCoeffDomain() != true ) 812 812 //alpha +=1; 813 813 } … … 817 817 // = (2) Find a suitable constant = 818 818 // =============================================== 819 819 820 820 Rg = resultant( f(alpha,y), g(alpha,y)-z*fx(alpha,y), x); 821 821 822 822 823 823 CFList LL; … … 838 838 CanonicalForm tmp; 839 839 840 factorsUni = AbsFactorize(ff); 840 factorsUni = AbsFactorize(ff); 841 841 842 842 for( CFFListIterator l=factorsUni; l.hasItem(); l++) … … 854 854 //======================================================= 855 855 CanonicalForm BIFAC::RationalFactor (CanonicalForm phi, CanonicalForm ff, \ 856 856 CanonicalForm fx, CanonicalForm g) 857 857 //======================================================= 858 858 { … … 866 866 867 867 hh = Bigcd(ff, h); 868 868 869 869 return(hh); 870 870 } … … 880 880 ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1 881 881 CanonicalForm phi = i.getItem().factor(); 882 882 883 883 if( ! phi.inBaseDomain()) 884 884 { … … 908 908 else 909 909 root = -tailcoeff(fac)/lc(fac); 910 911 910 911 912 912 AbsFac.append( f1(root,e) ); 913 913 AbsFac.append( i.getItem().exp() * exponent); … … 933 933 CanonicalForm h, h_abs, h_res, h_rat; 934 934 CanonicalForm fx = deriv(ff,x); 935 935 936 936 937 937 for( CFFListIterator i=Phis; i.hasItem(); i++) … … 939 939 ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1 940 940 phi = i.getItem().factor(); 941 941 942 942 if( ! phi.inBaseDomain()) 943 943 { … … 946 946 if( phi.degree() == 1 ) 947 947 { 948 949 950 951 952 953 954 955 956 gl_AL.append(h); // Factor of degree 1 957 958 948 if( taildegree(phi) > 0 ) // case: phi = a * x 949 h = gcd( ff,g ); 950 else // case: phi = a * x + c 951 { 952 h = gcd( ff, g+tailcoeff(phi)/lc(phi)*fx); 953 } 954 955 //biNormieren( h ); 956 gl_AL.append(h); // Factor of degree 1 957 gl_AL.append(exponent); // Multiplicity (exponent) 958 gl_AL.append(0); // No field extension 959 959 } else 960 960 { 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 // === all knowledge of absolute factors. === 993 h_rat = RationalFactor(phi, ff,fx, g); 994 995 996 997 961 // === Case 2: phi has degree > 1 === 962 e=rootOf(phi, 'e'); 963 h = gcd( ff, g-e*fx); 964 //biNormieren( h ); 965 966 AbsFac = getAbsoluteFactors(h, phi); 967 for( CFListIterator l=AbsFac; l.hasItem(); l++) 968 gl_AL.append( l.getItem() ); 969 970 971 // === (1) Get the rational factor by multi- === 972 // === plication of the absolute factor. === 973 h_abs=1; 974 ii = 0; 975 976 for( CFListIterator l=AbsFac; l.hasItem(); l++) 977 { 978 ii++; 979 if (ii%3 == 1 ) 980 h_abs *= l.getItem(); 981 } 982 //biNormieren( h_abs ); 983 984 985 // === (2) Compute the rational factor === 986 // === by using the resultant. === 987 h_res = resultant(phi(z,x), h(z,e), z); 988 //biNormieren( h_res ); 989 990 991 // === (3) Compute the rational factor by ignoring === 992 // === all knowledge of absolute factors. === 993 h_rat = RationalFactor(phi, ff,fx, g); 994 //biNormieren( h_rat ); 995 996 ASSERT( (h_abs == h_res) && (h_res == h_rat), "Wrong rational factor ?!?"); 997 h = h_abs; 998 998 } 999 999 // End of absolute factorization. … … 1005 1005 1006 1006 1007 //====================================================== 1007 //====================================================== 1008 1008 // Factorization of a squarefree bivariate polynomial 1009 1009 // in which every factor appears only once. 1010 1010 // Do we need a complete factorization ('absolute' is true) 1011 1011 // or only a rational factorization ('absolute' false)? 1012 //====================================================== 1012 //====================================================== 1013 1013 void BIFAC::bifacSqrFree(CanonicalForm ff) 1014 1014 //======================================================= … … 1035 1035 // LL = createEg(G,ff); 1036 1036 // LL = createEgUni(G,ff); // Hier ist noch ein FEHLER !!!! 1037 1037 1038 1038 LL = createRg( G, ff); // viel langsamer als EgUni 1039 1040 1039 1040 1041 1041 Eg = LL.getFirst(); 1042 Eg = Eg/LC(Eg); 1043 1042 Eg = Eg/LC(Eg); 1043 1044 1044 g = LL.getLast(); 1045 1045 1046 1046 // g = G.getFirst(); 1047 1047 1048 1049 CFFList PHI = AbsFactorize( Eg ); 1050 1051 CFFListIterator J=PHI; 1052 1053 1054 1048 1049 CFFList PHI = AbsFactorize( Eg ); 1050 1051 CFFListIterator J=PHI; 1052 CanonicalForm Eg2=1; 1053 for ( ; J.hasItem(); J++) 1054 { Eg2 = Eg2 * J.getItem().factor(); } 1055 1055 1056 1056 // === Is Eg(x) irreducible ? === 1057 1057 anz=0; 1058 1058 1059 // PHI = AbsFactorize( Eg) ; 1060 1061 1062 for( CFFListIterator i=PHI; i.hasItem(); i++) { 1059 // PHI = AbsFactorize( Eg) ; 1060 // 1061 1062 for( CFFListIterator i=PHI; i.hasItem(); i++) { 1063 1063 if( !i.getItem().factor().inBaseDomain()) 1064 1065 1064 anz++; 1065 } 1066 1066 1067 1067 /* if( absolute ) // Only for a absolute factorization 1068 1068 AbsoluteFactorization( PHI,ff, g); 1069 else // only for a rational factorization 1069 else // only for a rational factorization 1070 1070 { */ 1071 1071 if( anz==1 ){ ; 1072 1073 else 1074 1072 gl_RL.append( CFFactor(ff,exponent));} 1073 else 1074 RationalFactorizationOnly( PHI,ff, g); 1075 1075 /* } */ 1076 1076 } … … 1107 1107 // =============================================== 1108 1108 1109 1110 1111 // 1112 // 1109 1110 CFFList Q =Mysqrfree(f); 1111 // 1112 // cout << Q << endl; 1113 1113 // 1114 1114 … … 1123 1123 { 1124 1124 1125 1126 1127 1125 if( i.getItem().factor().level() < 0 ) ; 1126 else 1127 { 1128 1128 if( ( degree(i.getItem().factor(),x) == 0 || degree( i.getItem().factor(),y) == 0) ) { 1129 1129 // case: univariate … … 1131 1131 else // case: bivariate 1132 1132 { 1133 exponent = i.getItem().exp(); // global variable 1134 1135 1136 1137 bifacSqrFree(i.getItem().factor()/dumm ); 1133 exponent = i.getItem().exp(); // global variable 1134 CanonicalForm dumm = i.getItem().factor(); 1135 dumm = dumm.LC(); 1136 if( dumm.level() > 0 ){ dumm = 1; } 1137 bifacSqrFree(i.getItem().factor()/dumm ); 1138 1138 } 1139 1139 }} … … 1156 1156 1157 1157 if( min >= 32003 ) return ( 32003 ); // this is the maximum 1158 1158 1159 1159 // Find the smallest poosible prime 1160 1160 while ( cf_getPrime(nr) < min) { nr++; } … … 1184 1184 //cout << gl_RL<<endl; 1185 1185 1186 1187 1188 1189 1190 1191 1192 c = c(W,y); 1193 c = c(y,x); 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1186 if( sw ) 1187 { 1188 Variable W('W'); 1189 for( CFFListIterator i=gl_RL; i.hasItem(); i++) 1190 { 1191 c = i.getItem().factor(); 1192 c = c(W,y); 1193 c = c(y,x); 1194 c = c(x,W); 1195 aL.append( CFFactor( c, i.getItem().exp() )); 1196 } 1197 1198 f = f(W,y); f=f(y,x); f=f(x,W); 1199 } 1200 else aL = gl_RL; 1201 1202 gl_RL = aL; 1203 1204 //cout << aL; 1205 1205 1206 1206 … … 1217 1217 1218 1218 // cout << "\n* Test auf Korrektheit ..."; 1219 1220 1219 1220 1221 1221 for( CFFListIterator i=aL; i.hasItem(); i++) 1222 1222 { 1223 1223 ff *= power(i.getItem().factor(), i.getItem().exp() ); 1224 1224 // cout << " ff = " << ff 1225 // 1225 // << "\n a^b = " << i.getItem().factor() << " ^ " << i.getItem().exp() << endl; 1226 1226 } 1227 1227 c = f.LC()/ff.LC(); … … 1231 1231 1232 1232 // cout << "\n\nOriginal f = " << f << "\n\nff = " << ff 1233 // 1233 // << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl; 1234 1234 // cout << "degree 0: " << c << endl; 1235 1236 1235 1236 1237 1237 #ifndef NOSTREAMIO 1238 if( f != ff ) cout << "\n\nOriginal f = " << f << "\n\nff = " << ff 1239 1238 if( f != ff ) cout << "\n\nOriginal f = " << f << "\n\nff = " << ff 1239 << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl; 1240 1240 #endif 1241 1241 ASSERT( f==ff, "Wrong rational factorization. Abborting!"); 1242 1242 // cout << " [OK]\n"; 1243 1243 1244 1244 } 1245 1245 //--<>--------------------------------- … … 1254 1254 1255 1255 ASSERT( ch==0 && !isOn(SW_RATIONAL), "Integer numbers not allowed" ); 1256 1256 1257 1257 1258 1258 // === Check the characteristic === 1259 if( ch != 0 ) 1259 if( ch != 0 ) 1260 1260 { 1261 1261 ch2 = findCharacteristic(f); … … 1267 1267 // PROVISORISCH 1268 1268 //cerr << "\n Characteristic is too small!" 1269 1269 // << "\n The result might be wrong!\n\n"; 1270 1270 exit(1); 1271 1271 … … 1273 1273 } 1274 1274 1275 1276 1277 1278 1279 1280 f = f(W,x); f = f(x,y); f=f(y,W); 1281 1282 1283 1284 1285 1275 Variable W('W'); 1276 CanonicalForm l; 1277 int sw = 0; 1278 1279 if( degree(f,x) < degree(f,y) ) { 1280 f = f(W,x); f = f(x,y); f=f(y,W); 1281 sw = 1; 1282 } 1283 l = f.LC(); 1284 1285 if( l.level()<0 ) { f = f/f.LC(); gl_RL.append( CFFactor(l,1) ); } 1286 1286 1287 1287 -
factory/bifac.h
rc840d9 r806c18 44 44 //////////////////////////////////////////////////////////////// 45 45 public: 46 //////////////////////////////////////////////////////////////// 46 //////////////////////////////////////////////////////////////// 47 47 48 48 // === KONST-/ DESTRUKTOREN ==== … … 78 78 void unifac(CanonicalForm f, int grad); 79 79 CanonicalForm RationalFactor (CanonicalForm phi, CanonicalForm f, \ 80 80 CanonicalForm fx, CanonicalForm g); 81 81 void RationalFactorizationOnly (CFFList Phis, CanonicalForm f, CanonicalForm g); 82 82 CFList getAbsoluteFactors (CanonicalForm f1, CanonicalForm phi); … … 84 84 void bifacSqrFree( CanonicalForm f ); 85 85 void bifacMain(CanonicalForm f); 86 86 87 87 88 88 // === Variable ======= 89 89 CFFList gl_RL; // where to store the rational factorization 90 90 CFList gl_AL; // where to store the absolute factorization 91 bool absolute; // Compute an absolute factorization as well? 92 int exponent; // 91 bool absolute; // Compute an absolute factorization as well? 92 int exponent; // 93 93 }; 94 94 -
factory/bifacConfig.h
rc840d9 r806c18 1 /* =================================================================== 1 /* =================================================================== 2 2 GLOBAL COMPILE OPTIONS FOR BIFAC 3 3 =================================================================== */ … … 40 40 41 41 42 /* =================================================================== 42 /* =================================================================== 43 43 GLOBAL COMPILE OPTIONS FOR MULTIFAC 44 44 =================================================================== */ -
factory/canonicalform.cc
rc840d9 r806c18 1442 1442 return g.value->bgcdcoeff( f.value ); 1443 1443 else if ( what == INTMARK && ! cf_glob_switches.isOn( SW_RATIONAL ) ) 1444 1444 { 1445 1445 // calculate gcd using standard integer 1446 1446 // arithmetic … … 1452 1452 // swap fInt and gInt 1453 1453 if ( gInt > fInt ) 1454 1454 { 1455 1455 int swap = gInt; 1456 1456 gInt = fInt; … … 1460 1460 // now, 0 <= gInt <= fInt. Start the loop. 1461 1461 while ( gInt ) 1462 1462 { 1463 1463 // calculate (fInt, gInt) = (gInt, fInt%gInt) 1464 1464 int r = fInt % gInt; … … 1469 1469 return CanonicalForm( fInt ); 1470 1470 } 1471 1471 else 1472 1472 // we do not go for maximal speed for these stupid 1473 1473 // special cases -
factory/canonicalform.h
rc840d9 r806c18 49 49 inline int is_imm ( const InternalCF * const ptr ) 50 50 { 51 // returns 0 if ptr is not immediate 51 // returns 0 if ptr is not immediate 52 52 return ( ((int)((long)ptr)) & 3 ); 53 53 } … … 325 325 { 326 326 if ( f.level() > 0 ) 327 327 return power( f.mvar(), f.degree() ) * f.LC(); 328 328 else 329 329 return f; 330 330 } 331 331 -
factory/cf_binom.cc
rc840d9 r806c18 34 34 35 35 if ( ! initialized ) { 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 36 initialized = true; 37 ptZ = new CFArray[MAXPT+1]; 38 ptF = new CFArray[MAXPT+1]; 39 int i, j; 40 ptZ[0] = CFArray(1); ptZ[0][0] = 1; 41 ptF[0] = CFArray(1); 42 for ( i = 1; i <= INITPT; i++ ) { 43 ptF[i] = CFArray(i+1); 44 ptZ[i] = CFArray(i+1); 45 (ptZ[i])[0] = 1; 46 for ( j = 1; j < i; j++ ) 47 (ptZ[i])[j] = (ptZ[i-1])[j-1] + (ptZ[i-1])[j]; 48 (ptZ[i])[i] = 1; 49 } 50 for ( i = INITPT+1; i <= MAXPT; i++ ) { 51 ptF[i] = CFArray(i+1); 52 ptZ[i] = CFArray(i+1); 53 } 54 ptZmax = INITPT; 55 ptFmax = 0; 56 56 } 57 57 } … … 61 61 { 62 62 if ( n == 0 ) 63 63 return 1; 64 64 else if ( n == 1 ) 65 65 return x + a; 66 66 else if ( getCharacteristic() == 0 ) { 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 67 if ( n <= MAXPT ) { 68 if ( n > ptZmax ) { 69 int i, j; 70 for ( i = ptZmax+1; i <= n; i++ ) { 71 (ptZ[i])[0] = 1; 72 for ( j = 1; j < i; j++ ) 73 (ptZ[i])[j] = (ptZ[i-1])[j-1] + (ptZ[i-1])[j]; 74 (ptZ[i])[i] = 1; 75 } 76 ptZmax = n; 77 } 78 CanonicalForm result = 0, apower = 1; 79 int k; 80 for ( k = n; k >= 0; k-- ) { 81 result += power( x, k ) * apower * (ptZ[n])[k]; 82 if ( k != 0 ) 83 apower *= a; 84 } 85 return result; 86 } 87 else { 88 CanonicalForm result = binomialpower( x, a, MAXPT ); 89 CanonicalForm xa = x + a; 90 int i; 91 for ( i = MAXPT; i < n; i++ ) 92 result *= xa; 93 return result; 94 } 95 95 } 96 96 else { 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 97 if ( getCharacteristic() != charac || gfdeg != getGFDegree() ) { 98 ptFmax = 0; 99 charac = getCharacteristic(); 100 gfdeg = getGFDegree(); 101 (ptF[0])[0] = 1; 102 } 103 if ( n <= MAXPT ) { 104 if ( n > ptFmax ) { 105 int i, j; 106 for ( i = ptFmax+1; i <= n; i++ ) { 107 (ptF[i])[0] = 1; 108 for ( j = 1; j < i; j++ ) 109 (ptF[i])[j] = (ptF[i-1])[j-1] + (ptF[i-1])[j]; 110 (ptF[i])[i] = 1; 111 } 112 ptFmax = n; 113 } 114 CanonicalForm result = 0, apower = 1; 115 int k; 116 for ( k = n; k >= 0; k-- ) { 117 result += power( x, k ) * apower * (ptF[n])[k]; 118 if ( k != 0 ) 119 apower *= a; 120 } 121 return result; 122 } 123 else { 124 CanonicalForm result = binomialpower( x, a, MAXPT ); 125 CanonicalForm xa = x + a; 126 int i; 127 for ( i = MAXPT; i < n; i++ ) 128 result *= xa; 129 return result; 130 } 131 131 } 132 132 } -
factory/cf_char.cc
rc840d9 r806c18 38 38 theCharacteristic = c; 39 39 ff_big = c > cf_getSmallPrime( cf_getNumSmallPrimes()-1 ); 40 40 if (c > 536870909) factoryError("characteristic is too large(max is 2^29)"); 41 41 ff_setprime( c ); 42 42 resetFPT(); -
factory/cf_chinese.cc
rc840d9 r806c18 162 162 { 163 163 //assume(P>0); 164 // assume !isOn(SW_RATIONAL): mod is a no-op otherwise 164 // assume !isOn(SW_RATIONAL): mod is a no-op otherwise 165 165 if (N<0) N +=P; 166 166 CanonicalForm A,B,C,D,E; … … 171 171 if (2*N*N<P) 172 172 { 173 174 175 173 On(SW_RATIONAL); 174 N /=B; 175 Off(SW_RATIONAL); 176 176 return(N); 177 177 } -
factory/cf_cyclo.cc
rc840d9 r806c18 11 11 * (c) by The SINGULAR Team, see LICENSE file 12 12 * 13 * @internal 13 * @internal 14 14 * @version \$Id$ 15 15 * … … 31 31 /// function may fail if integer contains primes which exceed the largest prime 32 32 /// in our table 33 int* integerFactorizer (const long integer, int& length, bool& fail) 33 int* integerFactorizer (const long integer, int& length, bool& fail) 34 34 { 35 ASSERT (integer != 0 && integer != 1 && integer != -1, 36 "non-zero non-unit expected"); 35 ASSERT (integer != 0 && integer != 1 && integer != -1, 36 "non-zero non-unit expected"); 37 37 int* result; 38 38 length= 0; … … 43 43 44 44 int exp= 0; 45 while ((i != 1) && (i%2 == 0)) 45 while ((i != 1) && (i%2 == 0)) 46 46 { 47 47 i /= 2; 48 48 exp++; 49 49 } 50 if (exp != 0) 50 if (exp != 0) 51 51 { 52 52 result= new int [exp]; … … 54 54 result[k]= 2; 55 55 length += exp; 56 } 56 } 57 57 if (i == 1) return result; 58 58 … … 61 61 int* buf; 62 62 int next_prime; 63 while ((i != 1) && (j < 31937)) 63 while ((i != 1) && (j < 31937)) 64 64 { 65 65 next_prime= cf_getPrime (j); 66 while ((i != 1) && (i%next_prime == 0)) 66 while ((i != 1) && (i%next_prime == 0)) 67 67 { 68 68 i /= next_prime; 69 69 exp++; 70 } 71 if (exp != 0) 70 } 71 if (exp != 0) 72 72 { 73 buf= result; 73 buf= result; 74 74 result= new int [length + exp]; 75 for (int k= 0; k < length; k++) 75 for (int k= 0; k < length; k++) 76 76 result [k]= buf[k]; 77 77 for (int k= 0; k < exp; k++) 78 78 result [k + length]= next_prime; 79 length += exp; 79 length += exp; 80 80 } 81 81 exp= 0; … … 96 96 int* buf; 97 97 result[0]= factors [0]; 98 for (int i= 1; i < factors_length; i++) 98 for (int i= 1; i < factors_length; i++) 99 99 { 100 if (factors[i - 1] != factors[i]) 100 if (factors[i - 1] != factors[i]) 101 101 { 102 102 buf= result; … … 104 104 for (int j= 0; j < length; j++) 105 105 result[j]= buf [j]; 106 result[length]= factors[i]; 107 length++; 106 result[length]= factors[i]; 107 length++; 108 108 } 109 109 } … … 113 113 /// compute the n-th cyclotomic polynomial, 114 114 /// function may fail if integer_factorizer fails to factorize n 115 CanonicalForm cyclotomicPoly (int n, bool& fail) 115 CanonicalForm cyclotomicPoly (int n, bool& fail) 116 116 { 117 117 fail= false; … … 133 133 { 134 134 result= result (power (x, distinct_factors[i]), x)/result; 135 prod *= distinct_factors[i]; 135 prod *= distinct_factors[i]; 136 136 } 137 137 return result (power (x, n/prod), x); … … 140 140 #ifdef HAVE_NTL 141 141 /// checks if alpha is a primitive element, alpha is assumed to be an algebraic 142 /// variable over some finite prime field 143 bool isPrimitive (const Variable& alpha, bool& fail) 142 /// variable over some finite prime field 143 bool isPrimitive (const Variable& alpha, bool& fail) 144 144 { 145 145 int p= getCharacteristic(); -
factory/cf_cyclo.h
rc840d9 r806c18 13 13 * (c) by The SINGULAR Team, see LICENSE file 14 14 * 15 * @internal 15 * @internal 16 16 * @version \$Id$ 17 17 * -
factory/cf_eval.cc
rc840d9 r806c18 15 15 { 16 16 if ( this != &e ) { 17 17 values = e.values; 18 18 } 19 19 return *this; … … 24 24 { 25 25 if ( f.inCoeffDomain() || f.level() < values.min() ) 26 26 return f; 27 27 else if ( f.level() < values.max() ) 28 28 return evalCF( f, values, values.min(), f.level() ); 29 29 else 30 30 return evalCF( f, values, values.min(), values.max() ); 31 31 } 32 32 … … 35 35 { 36 36 if ( i > j ) 37 37 return f; 38 38 return evalCF( f, values, i, j ); 39 39 } … … 44 44 int n = values.max(); 45 45 for ( int i = values.min(); i <= n; i++ ) 46 46 values[i] += 1; 47 47 } 48 48 … … 60 60 { 61 61 if ( m > n ) 62 62 return f; 63 63 else { 64 65 66 67 68 69 64 CanonicalForm result = f; 65 while ( n >= m ) { 66 result = result( a[n], Variable( n ) ); 67 n--; 68 } 69 return result; 70 70 } 71 71 // iterated method turned out to be faster than -
factory/cf_factor.cc
rc840d9 r806c18 129 129 { 130 130 mpz_t m; 131 131 gmp_numerator(f,m); 132 132 char * str = new char[mpz_sizeinbase( m, 10 ) + 2]; 133 133 str = mpz_get_str( str, 10, m ); 134 134 printf("%s",str); 135 135 delete[] str; 136 136 mpz_clear(m); 137 137 } 138 138 else if (f.inQ()) 139 139 { 140 140 mpz_t m; 141 141 gmp_numerator(f,m); 142 142 char * str = new char[mpz_sizeinbase( m, 10 ) + 2]; 143 143 str = mpz_get_str( str, 10, m ); 144 144 printf("%s/",str); 145 145 delete[] str; 146 146 mpz_clear(m); 147 147 gmp_denominator(f,m); 148 148 str = new char[mpz_sizeinbase( m, 10 ) + 2]; … … 150 150 printf("%s",str); 151 151 delete[] str; 152 152 mpz_clear(m); 153 153 } 154 154 #else -
factory/cf_factory.cc
rc840d9 r806c18 29 29 { 30 30 if ( currenttype == IntegerDomain ) 31 32 33 34 31 if ( value >= MINIMMEDIATE && value <= MAXIMMEDIATE ) 32 return int2imm( value ); 33 else 34 return new InternalInteger( value ); 35 35 // else if ( currenttype == RationalDomain ) 36 // 37 // 38 // 39 // 36 // if ( value >= MINIMMEDIATE && value <= MAXIMMEDIATE ) 37 // return int2imm( value ); 38 // else 39 // return new InternalRational( value ); 40 40 else if ( currenttype == FiniteFieldDomain ) 41 41 return int2imm_p( ff_norm( value ) ); 42 42 else if ( currenttype == GaloisFieldDomain ) 43 43 return int2imm_gf( gf_int2gf( value ) ); 44 44 else if ( currenttype == PrimePowerDomain ) 45 46 else { 47 48 45 return new InternalPrimePower( value ); 46 else { 47 ASSERT( 0, "illegal basic domain!" ); 48 return 0; 49 49 } 50 50 } … … 54 54 { 55 55 if ( type == IntegerDomain ) 56 57 58 59 56 if ( value >= MINIMMEDIATE && value <= MAXIMMEDIATE ) 57 return int2imm( value ); 58 else 59 return new InternalInteger( value ); 60 60 // else if ( type == RationalDomain ) 61 // 62 // 63 // 64 // 61 // if ( value >= MINIMMEDIATE && value <= MAXIMMEDIATE ) 62 // return int2imm( value ); 63 // else 64 // return new InternalRational( value ); 65 65 else if ( type == FiniteFieldDomain ) 66 66 return int2imm_p( ff_norm( value ) ); 67 67 else if ( type == GaloisFieldDomain ) 68 68 return int2imm_gf( gf_int2gf( value ) ); 69 69 else if ( type == PrimePowerDomain ) 70 71 else { 72 73 70 return new InternalPrimePower( value ); 71 else { 72 ASSERT1( 0, "illegal basic domain (type = %d)!", type ); 73 return 0; 74 74 } 75 75 } … … 79 79 { 80 80 if ( currenttype == IntegerDomain ) { 81 82 83 84 85 86 87 88 81 InternalInteger * dummy = new InternalInteger( str ); 82 if ( dummy->is_imm() ) { 83 InternalCF * res = int2imm( dummy->intval() ); 84 delete dummy; 85 return res; 86 } 87 else 88 return dummy; 89 89 } 90 90 // else if ( currenttype == RationalDomain ) { 91 // 92 // 93 // 94 // 95 // 96 // 97 // 98 // 91 // InternalRational * dummy = new InternalRational( str ); 92 // if ( dummy->is_imm() ) { 93 // InternalCF * res = int2imm( dummy->intval() ); 94 // delete dummy; 95 // return res; 96 // } 97 // else 98 // return dummy; 99 99 // } 100 100 else if ( currenttype == FiniteFieldDomain ) { 101 102 103 104 101 InternalInteger * dummy = new InternalInteger( str ); 102 InternalCF * res = int2imm_p( dummy->intmod( ff_prime ) ); 103 delete dummy; 104 return res; 105 105 } 106 106 else if ( currenttype == GaloisFieldDomain ) { 107 108 109 110 107 InternalInteger * dummy = new InternalInteger( str ); 108 InternalCF * res = int2imm_gf( gf_int2gf( dummy->intmod( ff_prime ) ) ); 109 delete dummy; 110 return res; 111 111 } 112 112 else if ( currenttype == PrimePowerDomain ) 113 114 else { 115 116 113 return new InternalPrimePower( str ); 114 else { 115 ASSERT( 0, "illegal basic domain!" ); 116 return 0; 117 117 } 118 118 } … … 122 122 { 123 123 if ( currenttype == IntegerDomain ) { 124 125 126 127 128 129 130 131 124 InternalInteger * dummy = new InternalInteger( str, base ); 125 if ( dummy->is_imm() ) { 126 InternalCF * res = int2imm( dummy->intval() ); 127 delete dummy; 128 return res; 129 } 130 else 131 return dummy; 132 132 } 133 133 // else if ( currenttype == RationalDomain ) { 134 // 135 // 136 // 137 // 138 // 139 // 140 // 141 // 134 // InternalRational * dummy = new InternalRational( str ); 135 // if ( dummy->is_imm() ) { 136 // InternalCF * res = int2imm( dummy->intval() ); 137 // delete dummy; 138 // return res; 139 // } 140 // else 141 // return dummy; 142 142 // } 143 143 else if ( currenttype == FiniteFieldDomain ) { 144 145 146 147 144 InternalInteger * dummy = new InternalInteger( str, base ); 145 InternalCF * res = int2imm_p( dummy->intmod( ff_prime ) ); 146 delete dummy; 147 return res; 148 148 } 149 149 else if ( currenttype == GaloisFieldDomain ) { 150 151 152 153 150 InternalInteger * dummy = new InternalInteger( str, base ); 151 InternalCF * res = int2imm_gf( gf_int2gf( dummy->intmod( ff_prime ) ) ); 152 delete dummy; 153 return res; 154 154 } 155 155 else if ( currenttype == PrimePowerDomain ) 156 157 else { 158 159 156 return new InternalPrimePower( str, base ); 157 else { 158 ASSERT( 0, "illegal basic domain!" ); 159 return 0; 160 160 } 161 161 } … … 165 165 { 166 166 if ( type == IntegerDomain ) { 167 168 169 170 171 172 173 174 167 InternalInteger * dummy = new InternalInteger( str ); 168 if ( dummy->is_imm() ) { 169 InternalCF * res = int2imm( dummy->intval() ); 170 delete dummy; 171 return res; 172 } 173 else 174 return dummy; 175 175 } 176 176 // else if ( type == RationalDomain ) { 177 // 178 // 179 // 180 // 181 // 182 // 183 // 184 // 177 // InternalRational * dummy = new InternalRational( str ); 178 // if ( dummy->is_imm() ) { 179 // InternalCF * res = int2imm( dummy->intval() ); 180 // delete dummy; 181 // return res; 182 // } 183 // else 184 // return dummy; 185 185 // } 186 186 else if ( type == FiniteFieldDomain ) { 187 188 189 190 187 InternalInteger * dummy = new InternalInteger( str ); 188 InternalCF * res = int2imm( dummy->intmod( ff_prime ) ); 189 delete dummy; 190 return res; 191 191 } 192 192 else if ( type == GaloisFieldDomain ) { 193 194 195 196 193 InternalInteger * dummy = new InternalInteger( str ); 194 InternalCF * res = int2imm_gf( gf_int2gf( dummy->intmod( ff_prime ) ) ); 195 delete dummy; 196 return res; 197 197 } 198 198 else if ( type == PrimePowerDomain ) 199 200 else { 201 202 199 return new InternalPrimePower( str ); 200 else { 201 ASSERT( 0, "illegal basic domain!" ); 202 return 0; 203 203 } 204 204 } … … 208 208 { 209 209 if ( nonimm ) 210 211 212 213 214 215 216 217 218 else 219 210 if ( type == IntegerDomain ) 211 return new InternalInteger( value ); 212 else if ( type == RationalDomain ) 213 return new InternalRational( value ); 214 else { 215 ASSERT( 0, "illegal basic domain!" ); 216 return 0; 217 } 218 else 219 return CFFactory::basic( type, value ); 220 220 } 221 221 … … 224 224 { 225 225 if ( currenttype != IntegerDomain ) { 226 227 228 } 229 else 230 226 InternalPrimePower * dummy = new InternalPrimePower( num ); 227 return (InternalCF*)(dummy->normalize_myself()); 228 } 229 else 230 return new InternalInteger( num ); 231 231 } 232 232 … … 242 242 { 243 243 if ( normalize ) { 244 245 246 } 247 else 248 244 InternalRational * result = new InternalRational( num, den ); 245 return result->normalize_myself(); 246 } 247 else 248 return new InternalRational( num, den ); 249 249 } 250 250 … … 253 253 { 254 254 if ( v.level() == LEVELBASE ) 255 256 else 257 255 return c.getval(); 256 else 257 return new InternalPoly( v, exp, c ); 258 258 } 259 259 … … 262 262 { 263 263 if ( v.level() == LEVELBASE ) 264 265 else 266 264 return CFFactory::basic( 1 ); 265 else 266 return new InternalPoly( v, exp, 1 ); 267 267 } 268 268 … … 272 272 MP_INT dummy; 273 273 if ( value->levelcoeff() == IntegerDomain ) 274 274 mpz_init_set( &dummy, &InternalInteger::MPI( value ) ); 275 275 else if ( symmetric ) { 276 277 278 279 280 281 } 282 else 283 276 mpz_init( &dummy ); 277 if ( mpz_cmp( &InternalPrimePower::primepowhalf, &InternalPrimePower::MPI( value ) ) < 0 ) 278 mpz_sub( &dummy, &InternalPrimePower::MPI( value ), &InternalPrimePower::primepow ); 279 else 280 mpz_set( &dummy, &InternalPrimePower::MPI( value ) ); 281 } 282 else 283 mpz_init_set( &dummy, &InternalPrimePower::MPI( value ) ); 284 284 return dummy; 285 285 } -
factory/cf_gcd_charp.cc
rc840d9 r806c18 41 41 42 42 43 static CanonicalForm 43 static CanonicalForm 44 44 contentWRT0(const CanonicalForm & F, const int lev, const int lev0) 45 45 // Computes the content of a polynomial, considering the variables of level … … 64 64 { 65 65 CanonicalForm cc=i.coeff(); 66 if (cc.level() > lev0) 66 if (cc.level() > lev0) 67 67 pol = contentWRT0(cc, lev, lev0-1 ); 68 68 else … … 79 79 } 80 80 81 static CanonicalForm 81 static CanonicalForm 82 82 contentWRT(const CanonicalForm & F, const int lev) 83 83 // Computes the content of a polynomial, considering the variables of level … … 108 108 { 109 109 CanonicalForm cc=i.coeff(); 110 if (cc.level() > lev) 110 if (cc.level() > lev) 111 111 pol = contentWRT0(cc, lev, cc.level()); 112 112 else … … 364 364 // cout << "pow: " << ipower(p, k) << endl; 365 365 //} 366 366 367 367 CFMap M,N; 368 368 compress(A,B,M,N); … … 538 538 On(SW_USE_GCD_P); 539 539 return temp; 540 } 540 } 541 541 } 542 542 else … … 560 560 On(SW_USE_GCD_P); 561 561 return temp; 562 } 562 } 563 563 } 564 564 } -
factory/cf_gcd_smallp.cc
rc840d9 r806c18 6 6 * @date 22.10.2009 7 7 * 8 * This file implements the GCD of two polynomials over \f$ F_{p} \f$ , 9 * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms 8 * This file implements the GCD of two polynomials over \f$ F_{p} \f$ , 9 * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms 10 10 * for Computer Algebra" by Geddes, Czapor, Labahnn 11 11 * … … 13 13 * (c) by The SINGULAR Team, see LICENSE file 14 14 * 15 * @internal 15 * @internal 16 16 * @version \$Id$ 17 17 * … … 46 46 TIMING_DEFINE_PRINT(newton_interpolation); 47 47 48 /// compressing two polynomials F and G, M is used for compressing, 48 /// compressing two polynomials F and G, M is used for compressing, 49 49 /// N to reverse the compression 50 50 static inline 51 51 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, 52 CFMap & N, bool& topLevel) 53 { 52 CFMap & N, bool& topLevel) 53 { 54 54 int n= tmax (F.level(), G.level()); 55 55 int * degsf= new int [n + 1]; … … 58 58 for (int i = 0; i <= n; i++) 59 59 degsf[i]= degsg[i]= 0; 60 60 61 61 degsf= degrees (F, degsf); 62 62 degsg= degrees (G, degsg); 63 63 64 64 int both_non_zero= 0; 65 65 int f_zero= 0; … … 67 67 int both_zero= 0; 68 68 69 if (topLevel) 70 { 71 for (int i= 1; i <= n; i++) 72 { 73 if (degsf[i] != 0 && degsg[i] != 0) 69 if (topLevel) 70 { 71 for (int i= 1; i <= n; i++) 72 { 73 if (degsf[i] != 0 && degsg[i] != 0) 74 74 { 75 75 both_non_zero++; 76 76 continue; 77 77 } 78 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 78 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 79 79 { 80 80 f_zero++; 81 81 continue; 82 82 } 83 if (degsg[i] == 0 && degsf[i] && i <= F.level()) 83 if (degsg[i] == 0 && degsf[i] && i <= F.level()) 84 84 { 85 85 g_zero++; … … 88 88 } 89 89 90 if (both_non_zero == 0) 90 if (both_non_zero == 0) 91 91 { 92 92 delete [] degsf; … … 98 98 int k= 1; 99 99 int l= 1; 100 for (int i= 1; i <= n; i++) 101 { 102 if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) 103 { 104 if (k + both_non_zero != i) 100 for (int i= 1; i <= n; i++) 101 { 102 if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) 103 { 104 if (k + both_non_zero != i) 105 105 { 106 106 M.newpair (Variable (i), Variable (k + both_non_zero)); … … 109 109 k++; 110 110 } 111 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 112 { 113 if (l + g_zero + both_non_zero != i) 111 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 112 { 113 if (l + g_zero + both_non_zero != i) 114 114 { 115 115 M.newpair (Variable (i), Variable (l + g_zero + both_non_zero)); … … 119 119 } 120 120 } 121 121 122 122 // sort Variables x_{i} in increasing order of 123 // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) 123 // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) 124 124 int m= tmin (F.level(), G.level()); 125 125 int max_min_deg; … … 127 127 l= 0; 128 128 int i= 1; 129 while (k > 0) 129 while (k > 0) 130 130 { 131 131 max_min_deg= tmin (degsf[i], degsg[i]); 132 while (max_min_deg == 0) 132 while (max_min_deg == 0) 133 133 { 134 134 i++; 135 135 max_min_deg= tmin (degsf[i], degsg[i]); 136 136 } 137 for (int j= i + 1; j <= m; j++) 138 { 139 if (tmin (degsf[j],degsg[j]) >= max_min_deg) 137 for (int j= i + 1; j <= m; j++) 138 { 139 if (tmin (degsf[j],degsg[j]) >= max_min_deg) 140 140 { 141 141 max_min_deg= tmin (degsf[j], degsg[j]); 142 l= j; 142 l= j; 143 143 } 144 144 } 145 if (l != 0) 146 { 147 if (l != k) 145 if (l != 0) 146 { 147 if (l != k) 148 148 { 149 149 M.newpair (Variable (l), Variable(k)); … … 153 153 l= 0; 154 154 } 155 else 155 else 156 156 { 157 157 degsf[l]= 0; … … 159 159 l= 0; 160 160 } 161 } 162 else if (l == 0) 163 { 164 if (i != k) 161 } 162 else if (l == 0) 163 { 164 if (i != k) 165 165 { 166 166 M.newpair (Variable (i), Variable (k)); … … 169 169 degsg[i]= 0; 170 170 } 171 else 171 else 172 172 { 173 173 degsf[i]= 0; … … 175 175 } 176 176 i++; 177 } 177 } 178 178 k--; 179 179 } 180 180 } 181 else 181 else 182 182 { 183 183 //arrange Variables such that no gaps occur 184 for (int i= 1; i <= n; i++) 185 { 186 if (degsf[i] == 0 && degsg[i] == 0) 184 for (int i= 1; i <= n; i++) 185 { 186 if (degsf[i] == 0 && degsg[i] == 0) 187 187 { 188 188 both_zero++; 189 189 continue; 190 190 } 191 else 192 { 193 if (both_zero != 0) 191 else 192 { 193 if (both_zero != 0) 194 194 { 195 195 M.newpair (Variable (i), Variable (i - both_zero)); … … 207 207 208 208 int 209 substituteCheck (const CanonicalForm& F, const CanonicalForm& G) 209 substituteCheck (const CanonicalForm& F, const CanonicalForm& G) 210 210 { 211 211 if (F.inCoeffDomain() || G.inCoeffDomain()) … … 214 214 if (degree (F, x) <= 1 || degree (G, x) <= 1) 215 215 return 0; 216 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive 217 CanonicalForm g= swapvar (G, G.mvar(), x); 216 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive 217 CanonicalForm g= swapvar (G, G.mvar(), x); 218 218 int sizef= 0; 219 int sizeg= 0; 219 int sizeg= 0; 220 220 for (CFIterator i= f; i.hasTerms(); i++, sizef++) 221 221 { … … 240 240 expg [j]= i.exp(); 241 241 } 242 242 243 243 int indf= sizef - 1; 244 244 int indg= sizeg - 1; … … 247 247 if (expg[indg] == 0) 248 248 indg--; 249 249 250 250 if (expg[indg] != expf [indf] || (expg[indg] == 1 && expf[indf] == 1)) 251 251 { … … 265 265 } 266 266 } 267 267 268 268 for (int i= indg - 1; i >= 0; i--) 269 269 { … … 282 282 283 283 // substiution 284 void 285 subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A, 284 void 285 subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A, 286 286 CanonicalForm& B, const int d 287 287 ) … … 293 293 return; 294 294 } 295 295 296 296 CanonicalForm C= 0; 297 CanonicalForm D= 0; 297 CanonicalForm D= 0; 298 298 Variable x= Variable (1); 299 299 CanonicalForm f= swapvar (F, x, F.mvar()); … … 307 307 } 308 308 309 CanonicalForm 309 CanonicalForm 310 310 reverseSubst (const CanonicalForm& F, const int d) 311 311 { 312 312 if (d == 1) 313 313 return F; 314 314 315 315 Variable x= Variable (1); 316 316 if (degree (F, x) <= 0) … … 320 320 for (CFIterator i= f; i.hasTerms(); i++) 321 321 result += i.coeff()*power (f.mvar(), d*i.exp()); 322 return swapvar (result, x, F.mvar()); 323 } 324 325 static inline CanonicalForm 322 return swapvar (result, x, F.mvar()); 323 } 324 325 static inline CanonicalForm 326 326 uni_content (const CanonicalForm & F); 327 327 … … 335 335 if (F.level() != x.level() && F.isUnivariate()) 336 336 return F.genOne(); 337 337 338 338 if (x.level() != 1) 339 339 { … … 344 344 else 345 345 return uni_content (F); 346 } 347 348 /// compute the content of F, where F is considered as an element of 349 /// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$ 350 static inline CanonicalForm 351 uni_content (const CanonicalForm & F) 352 { 346 } 347 348 /// compute the content of F, where F is considered as an element of 349 /// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$ 350 static inline CanonicalForm 351 uni_content (const CanonicalForm & F) 352 { 353 353 if (F.inBaseDomain()) 354 354 return F.genOne(); … … 360 360 361 361 int l= F.level(); 362 if (l == 2) 362 if (l == 2) 363 363 return content(F); 364 else 364 else 365 365 { 366 366 CanonicalForm pol, c = 0; 367 367 CFIterator i = F; 368 for (; i.hasTerms(); i++) 369 { 370 pol= i.coeff(); 368 for (; i.hasTerms(); i++) 369 { 370 pol= i.coeff(); 371 371 pol= uni_content (pol); 372 372 c= gcd (c, pol); … … 378 378 } 379 379 380 CanonicalForm 381 extractContents (const CanonicalForm& F, const CanonicalForm& G, 382 CanonicalForm& contentF, CanonicalForm& contentG, 380 CanonicalForm 381 extractContents (const CanonicalForm& F, const CanonicalForm& G, 382 CanonicalForm& contentF, CanonicalForm& contentG, 383 383 CanonicalForm& ppF, CanonicalForm& ppG, const int d) 384 384 { … … 407 407 /// is dp. 408 408 static inline 409 CanonicalForm uni_lcoeff (const CanonicalForm& F) 409 CanonicalForm uni_lcoeff (const CanonicalForm& F) 410 410 { 411 411 if (F.level() <= 1) 412 return F; 412 return F; 413 413 else 414 414 { … … 418 418 { 419 419 if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg) 420 return uni_lcoeff (i.coeff()); 420 return uni_lcoeff (i.coeff()); 421 421 } 422 422 } … … 441 441 } 442 442 443 /// compute a random element a of \f$ F_{p}(\alpha ) \f$ , 443 /// compute a random element a of \f$ F_{p}(\alpha ) \f$ , 444 444 /// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns 445 /// fail if there are no field elements left which have not been used before 446 static inline CanonicalForm 445 /// fail if there are no field elements left which have not been used before 446 static inline CanonicalForm 447 447 randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list, 448 bool & fail) 448 bool & fail) 449 449 { 450 450 fail= false; … … 457 457 int d= degree (mipo); 458 458 int bound= ipower (p, d); 459 do 459 do 460 460 { 461 461 if (list.length() == bound) … … 464 464 break; 465 465 } 466 if (list.length() < p) 466 if (list.length() < p) 467 467 { 468 468 random= genFF.generate(); … … 470 470 random= genFF.generate(); 471 471 } 472 else 472 else 473 473 { 474 474 random= genAlgExt.generate(); … … 476 476 random= genAlgExt.generate(); 477 477 } 478 if (F (random, x) == 0) 478 if (F (random, x) == 0) 479 479 { 480 480 list.append (random); … … 485 485 } 486 486 487 /// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$ 487 /// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$ 488 488 /// we do not extend \f$ F_{p}(\alpha ) \f$ itself but extend \f$ F_{p} \f$ , 489 /// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$ 489 /// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$ 490 490 static inline 491 void choose_extension (const int& d, const int& num_vars, Variable& beta) 491 void choose_extension (const int& d, const int& num_vars, Variable& beta) 492 492 { 493 493 int p= getCharacteristic(); … … 496 496 ZZ_pX NTLirredpoly; 497 497 //TODO: replace d by max_{i} (deg_x{i}(f)) 498 int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); 498 int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); 499 499 int m= degree (getMipo (beta)); 500 500 if (i <= 1) 501 501 i= 2; 502 BuildIrred (NTLirredpoly, i*m); 503 CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); 504 beta= rootOf (mipo); 505 } 506 507 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 502 BuildIrred (NTLirredpoly, i*m); 503 CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); 504 beta= rootOf (mipo); 505 } 506 507 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 508 508 /// l and topLevel are only used internally, output is monic 509 509 /// based on Alg. 7.2. as described in "Algorithms for 510 510 /// Computer Algebra" by Geddes, Czapor, Labahn 511 CanonicalForm 512 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 513 Variable & alpha, CFList& l, bool& topLevel) 514 { 511 CanonicalForm 512 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 513 Variable & alpha, CFList& l, bool& topLevel) 514 { 515 515 CanonicalForm A= F; 516 516 CanonicalForm B= G; … … 522 522 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); 523 523 if (F == G) return F/Lc(F); 524 524 525 525 CFMap M,N; 526 526 int best_level= myCompress (A, B, M, N, topLevel); … … 533 533 Variable x= Variable(1); 534 534 535 //univariate case 536 if (A.isUnivariate() && B.isUnivariate()) 537 return N (gcd(A,B)); 538 535 //univariate case 536 if (A.isUnivariate() && B.isUnivariate()) 537 return N (gcd(A,B)); 538 539 539 int substitute= substituteCheck (A, B); 540 540 541 541 if (substitute > 1) 542 542 subst (A, B, A, B, substitute); … … 550 550 if (best_level <= 2) 551 551 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 552 else 553 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 552 else 553 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 554 554 } 555 555 else 556 556 { 557 557 cA = uni_content (A); 558 cB = uni_content (B); 558 cB = uni_content (B); 559 559 gcdcAcB= gcd (cA, cB); 560 560 ppA= A/cA; … … 563 563 564 564 CanonicalForm lcA, lcB; // leading coefficients of A and B 565 CanonicalForm gcdlcAlcB; 565 CanonicalForm gcdlcAlcB; 566 566 567 567 lcA= uni_lcoeff (ppA); 568 568 lcB= uni_lcoeff (ppB); 569 570 if (fdivides (lcA, lcB)) 571 { 569 570 if (fdivides (lcA, lcB)) 571 { 572 572 if (fdivides (A, B)) 573 return F/Lc(F); 573 return F/Lc(F); 574 574 } 575 575 if (fdivides (lcB, lcA)) 576 { 577 if (fdivides (B, A)) 576 { 577 if (fdivides (B, A)) 578 578 return G/Lc(G); 579 579 } … … 587 587 if (substitute > 1) 588 588 return N (reverseSubst (gcdcAcB, substitute)); 589 else 589 else 590 590 return N(gcdcAcB); 591 591 } 592 592 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 593 if (d0 < d) 594 d= d0; 593 if (d0 < d) 594 d= d0; 595 595 if (d == 0) 596 596 { 597 597 if (substitute > 1) 598 598 return N (reverseSubst (gcdcAcB, substitute)); 599 else 599 else 600 600 return N(gcdcAcB); 601 601 } … … 614 614 CanonicalForm prim_elem, im_prim_elem; 615 615 CFList source, dest; 616 do 616 do 617 617 { 618 618 random_element= randomElement (m, V_buf, l, fail); 619 if (fail) 619 if (fail) 620 620 { 621 621 source= CFList(); … … 638 638 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2)); 639 639 inextension= true; 640 for (CFListIterator i= l; i.hasItem(); i++) 641 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 640 for (CFListIterator i= l; i.hasItem(); i++) 641 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 642 642 im_prim_elem, source, dest); 643 643 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 644 644 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 645 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 645 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 646 646 source, dest); 647 647 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 648 648 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 649 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 649 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 650 650 source, dest); 651 651 … … 655 655 CFList list; 656 656 TIMING_START (gcd_recursion); 657 G_random_element= 657 G_random_element= 658 658 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 659 659 list, topLevel); 660 TIMING_END_AND_PRINT (gcd_recursion, 660 TIMING_END_AND_PRINT (gcd_recursion, 661 661 "time for recursive call: "); 662 662 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 663 663 } 664 else 664 else 665 665 { 666 666 CFList list; 667 667 TIMING_START (gcd_recursion); 668 G_random_element= 668 G_random_element= 669 669 GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf, 670 670 list, topLevel); 671 TIMING_END_AND_PRINT (gcd_recursion, 671 TIMING_END_AND_PRINT (gcd_recursion, 672 672 "time for recursive call: "); 673 673 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 674 674 } 675 675 676 d0= totaldegree (G_random_element, Variable(2), 676 d0= totaldegree (G_random_element, Variable(2), 677 677 Variable (G_random_element.level())); 678 678 if (d0 == 0) … … 680 680 if (substitute > 1) 681 681 return N (reverseSubst (gcdcAcB, substitute)); 682 else 682 else 683 683 return N(gcdcAcB); 684 684 } 685 if (d0 > d) 685 if (d0 > d) 686 686 { 687 687 if (!find (l, random_element)) … … 690 690 } 691 691 692 G_random_element= 692 G_random_element= 693 693 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element)) 694 694 * G_random_element; 695 695 696 d0= totaldegree (G_random_element, Variable(2), 696 d0= totaldegree (G_random_element, Variable(2), 697 697 Variable(G_random_element.level())); 698 if (d0 < d) 698 if (d0 < d) 699 699 { 700 700 m= gcdlcAlcB; … … 706 706 TIMING_START (newton_interpolation); 707 707 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); 708 TIMING_END_AND_PRINT (newton_interpolation, 708 TIMING_END_AND_PRINT (newton_interpolation, 709 709 "time for newton interpolation: "); 710 710 711 //termination test 712 if (uni_lcoeff (H) == gcdlcAlcB) 711 //termination test 712 if (uni_lcoeff (H) == gcdlcAlcB) 713 713 { 714 714 cH= uni_content (H); 715 715 ppH= H/cH; 716 if (inextension) 717 { 718 CFList u, v; 716 if (inextension) 717 { 718 CFList u, v; 719 719 //maybe it's better to test if ppH is an element of F(\alpha) before 720 720 //mapping down 721 721 DEBOUTLN (cerr, "ppH before mapDown= " << ppH); 722 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 722 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 723 723 ppH /= Lc(ppH); 724 724 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 725 if (fdivides (ppH, A) && fdivides (ppH, B)) 725 if (fdivides (ppH, A) && fdivides (ppH, B)) 726 726 { 727 727 if (substitute > 1) … … 733 733 } 734 734 } 735 else if (fdivides (ppH, A) && fdivides (ppH, B)) 735 else if (fdivides (ppH, A) && fdivides (ppH, B)) 736 736 { 737 737 if (substitute > 1) … … 752 752 } 753 753 754 /// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a 754 /// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a 755 755 /// univariate polynomial, returns fail if there are no field elements left 756 756 /// which have not been used before … … 766 766 int d= getGFDegree(); 767 767 int bound= ipower (p, d); 768 do 768 do 769 769 { 770 770 if (list.length() == bound) … … 775 775 if (list.length() < 1) 776 776 random= 0; 777 else 777 else 778 778 { 779 779 random= genGF.generate(); … … 781 781 random= genGF.generate(); 782 782 } 783 if (F (random, x) == 0) 783 if (F (random, x) == 0) 784 784 { 785 785 list.append (random); … … 796 796 /// the size of the base field is bounded by 2^16, output is monic 797 797 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 798 CFList& l, bool& topLevel) 799 { 798 CFList& l, bool& topLevel) 799 { 800 800 CanonicalForm A= F; 801 801 CanonicalForm B= G; … … 807 807 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); 808 808 if (F == G) return F/Lc(F); 809 809 810 810 CFMap M,N; 811 811 int best_level= myCompress (A, B, M, N, topLevel); … … 818 818 Variable x= Variable(1); 819 819 820 //univariate case 821 if (A.isUnivariate() && B.isUnivariate()) 822 return N (gcd (A, B)); 820 //univariate case 821 if (A.isUnivariate() && B.isUnivariate()) 822 return N (gcd (A, B)); 823 823 824 824 int substitute= substituteCheck (A, B); 825 825 826 826 if (substitute > 1) 827 827 subst (A, B, A, B, substitute); … … 835 835 if (best_level <= 2) 836 836 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 837 else 838 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 837 else 838 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 839 839 } 840 840 else 841 841 { 842 842 cA = uni_content (A); 843 cB = uni_content (B); 843 cB = uni_content (B); 844 844 gcdcAcB= gcd (cA, cB); 845 845 ppA= A/cA; … … 848 848 849 849 CanonicalForm lcA, lcB; // leading coefficients of A and B 850 CanonicalForm gcdlcAlcB; 850 CanonicalForm gcdlcAlcB; 851 851 852 852 lcA= uni_lcoeff (ppA); 853 853 lcB= uni_lcoeff (ppB); 854 854 855 if (fdivides (lcA, lcB)) 856 { 855 if (fdivides (lcA, lcB)) 856 { 857 857 if (fdivides (A, B)) 858 return F; 859 } 860 if (fdivides (lcB, lcA)) 861 { 862 if (fdivides (B, A)) 858 return F; 859 } 860 if (fdivides (lcB, lcA)) 861 { 862 if (fdivides (B, A)) 863 863 return G; 864 864 } … … 871 871 if (substitute > 1) 872 872 return N (reverseSubst (gcdcAcB, substitute)); 873 else 873 else 874 874 return N(gcdcAcB); 875 875 } 876 876 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 877 if (d0 < d) 878 d= d0; 877 if (d0 < d) 878 d= d0; 879 879 if (d == 0) 880 880 { 881 881 if (substitute > 1) 882 882 return N (reverseSubst (gcdcAcB, substitute)); 883 else 883 else 884 884 return N(gcdcAcB); 885 885 } … … 900 900 int expon; 901 901 char gf_name_buf= gf_name; 902 do 902 do 903 903 { 904 904 random_element= GFRandomElement (m, l, fail); 905 if (fail) 906 { 905 if (fail) 906 { 907 907 int num_vars= tmin (getNumVars(A), getNumVars(B)); 908 908 num_vars--; … … 911 911 if (expon < 2) 912 912 expon= 2; 913 kk= getGFDegree(); 914 if (ipower (p, kk*expon) < (1 << 16)) 913 kk= getGFDegree(); 914 if (ipower (p, kk*expon) < (1 << 16)) 915 915 setCharacteristic (p, kk*(int)expon, 'b'); 916 else 916 else 917 917 { 918 918 expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk)); … … 922 922 inextension= true; 923 923 fail= false; 924 for (CFListIterator i= l; i.hasItem(); i++) 924 for (CFListIterator i= l; i.hasItem(); i++) 925 925 i.getItem()= GFMapUp (i.getItem(), kk); 926 926 m= GFMapUp (m, kk); … … 934 934 CFList list; 935 935 TIMING_START (gcd_recursion); 936 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 936 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 937 937 list, topLevel); 938 TIMING_END_AND_PRINT (gcd_recursion, 938 TIMING_END_AND_PRINT (gcd_recursion, 939 939 "time for recursive call: "); 940 940 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 941 941 } 942 else 942 else 943 943 { 944 944 CFList list; 945 945 TIMING_START (gcd_recursion); 946 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 946 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 947 947 list, topLevel); 948 TIMING_END_AND_PRINT (gcd_recursion, 948 TIMING_END_AND_PRINT (gcd_recursion, 949 949 "time for recursive call: "); 950 950 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 951 951 } 952 952 953 d0= totaldegree (G_random_element, Variable(2), 953 d0= totaldegree (G_random_element, Variable(2), 954 954 Variable (G_random_element.level())); 955 if (d0 == 0) 956 { 957 if (inextension) 955 if (d0 == 0) 956 { 957 if (inextension) 958 958 { 959 959 setCharacteristic (p, k, gf_name_buf); … … 962 962 return N (reverseSubst (gcdcAcB, substitute)); 963 963 else 964 return N(gcdcAcB); 965 } 964 return N(gcdcAcB); 965 } 966 966 } 967 967 else … … 970 970 return N (reverseSubst (gcdcAcB, substitute)); 971 971 else 972 return N(gcdcAcB); 973 } 974 } 975 if (d0 > d) 972 return N(gcdcAcB); 973 } 974 } 975 if (d0 > d) 976 976 { 977 977 if (!find (l, random_element)) … … 980 980 } 981 981 982 G_random_element= 982 G_random_element= 983 983 (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) * 984 984 G_random_element; 985 d0= totaldegree (G_random_element, Variable(2), 985 d0= totaldegree (G_random_element, Variable(2), 986 986 Variable (G_random_element.level())); 987 987 988 if (d0 < d) 988 if (d0 < d) 989 989 { 990 990 m= gcdlcAlcB; … … 998 998 TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: "); 999 999 1000 //termination test 1001 if (uni_lcoeff (H) == gcdlcAlcB) 1000 //termination test 1001 if (uni_lcoeff (H) == gcdlcAlcB) 1002 1002 { 1003 1003 cH= uni_content (H); 1004 1004 ppH= H/cH; 1005 if (inextension) 1006 { 1007 if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) 1005 if (inextension) 1006 { 1007 if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) 1008 1008 { 1009 1009 DEBOUTLN (cerr, "ppH before mapDown= " << ppH); … … 1019 1019 } 1020 1020 } 1021 else 1022 { 1023 if (fdivides (ppH, A) && fdivides (ppH, B)) 1021 else 1022 { 1023 if (fdivides (ppH, A) && fdivides (ppH, B)) 1024 1024 { 1025 1025 if (substitute > 1) … … 1042 1042 1043 1043 /// F is assumed to be an univariate polynomial in x, 1044 /// computes a random monic irreducible univariate polynomial of random 1044 /// computes a random monic irreducible univariate polynomial of random 1045 1045 /// degree < i in x which does not divide F 1046 CanonicalForm 1047 randomIrredpoly (int i, const Variable & x) 1046 CanonicalForm 1047 randomIrredpoly (int i, const Variable & x) 1048 1048 { 1049 1049 int p= getCharacteristic(); 1050 1050 ZZ NTLp= to_ZZ (p); 1051 1051 ZZ_p::init (NTLp); 1052 ZZ_pX NTLirredpoly; 1052 ZZ_pX NTLirredpoly; 1053 1053 CanonicalForm CFirredpoly; 1054 1054 BuildIrred (NTLirredpoly, i + 1); 1055 1055 CFirredpoly= convertNTLZZpX2CF (NTLirredpoly, x); 1056 1056 return CFirredpoly; 1057 } 1057 } 1058 1058 1059 1059 static inline … … 1067 1067 int p= getCharacteristic(); 1068 1068 int bound= p; 1069 do 1069 do 1070 1070 { 1071 1071 if (list.length() == bound) … … 1076 1076 if (list.length() < 1) 1077 1077 random= 0; 1078 else 1078 else 1079 1079 { 1080 1080 random= genFF.generate(); … … 1082 1082 random= genFF.generate(); 1083 1083 } 1084 if (F (random, x) == 0) 1084 if (F (random, x) == 0) 1085 1085 { 1086 1086 list.append (random); … … 1091 1091 } 1092 1092 1093 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1094 bool& topLevel, CFList& l) 1093 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1094 bool& topLevel, CFList& l) 1095 1095 { 1096 1096 CanonicalForm A= F; … … 1114 1114 Variable x= Variable (1); 1115 1115 1116 //univariate case 1117 if (A.isUnivariate() && B.isUnivariate()) 1116 //univariate case 1117 if (A.isUnivariate() && B.isUnivariate()) 1118 1118 return N (gcd (A, B)); 1119 1119 1120 1120 int substitute= substituteCheck (A, B); 1121 1121 1122 1122 if (substitute > 1) 1123 1123 subst (A, B, A, B, substitute); … … 1131 1131 if (best_level <= 2) 1132 1132 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 1133 else 1134 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 1133 else 1134 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 1135 1135 } 1136 1136 else 1137 1137 { 1138 1138 cA = uni_content (A); 1139 cB = uni_content (B); 1139 cB = uni_content (B); 1140 1140 gcdcAcB= gcd (cA, cB); 1141 1141 ppA= A/cA; … … 1144 1144 1145 1145 CanonicalForm lcA, lcB; // leading coefficients of A and B 1146 CanonicalForm gcdlcAlcB; 1146 CanonicalForm gcdlcAlcB; 1147 1147 lcA= uni_lcoeff (ppA); 1148 1148 lcB= uni_lcoeff (ppB); 1149 1149 1150 if (fdivides (lcA, lcB)) 1151 { 1150 if (fdivides (lcA, lcB)) 1151 { 1152 1152 if (fdivides (A, B)) 1153 1153 return F/Lc(F); 1154 } 1154 } 1155 1155 if (fdivides (lcB, lcA)) 1156 { 1157 if (fdivides (B, A)) 1156 { 1157 if (fdivides (B, A)) 1158 1158 return G/Lc(G); 1159 1159 } 1160 1160 1161 gcdlcAlcB= gcd (lcA, lcB); 1162 1161 gcdlcAlcB= gcd (lcA, lcB); 1162 1163 1163 int d= totaldegree (ppA, Variable (2), Variable (ppA.level())); 1164 1164 int d0; … … 1168 1168 if (substitute > 1) 1169 1169 return N (reverseSubst (gcdcAcB, substitute)); 1170 else 1170 else 1171 1171 return N(gcdcAcB); 1172 1172 } 1173 1173 d0= totaldegree (ppB, Variable (2), Variable (ppB.level())); 1174 1174 1175 if (d0 < d) 1175 if (d0 < d) 1176 1176 d= d0; 1177 1177 1178 if (d == 0) 1178 if (d == 0) 1179 1179 { 1180 1180 if (substitute > 1) 1181 1181 return N (reverseSubst (gcdcAcB, substitute)); 1182 else 1182 else 1183 1183 return N(gcdcAcB); 1184 1184 } … … 1195 1195 topLevel= false; 1196 1196 CFList source, dest; 1197 do 1197 do 1198 1198 { 1199 1199 if (inextension) … … 1206 1206 CFList list; 1207 1207 TIMING_START (gcd_recursion); 1208 G_random_element= 1209 GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 1208 G_random_element= 1209 GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 1210 1210 list); 1211 TIMING_END_AND_PRINT (gcd_recursion, 1211 TIMING_END_AND_PRINT (gcd_recursion, 1212 1212 "time for recursive call: "); 1213 1213 DEBOUTLN (cerr, "G_random_element= " << G_random_element); … … 1215 1215 else if (!fail && inextension) 1216 1216 { 1217 CFList list; 1217 CFList list; 1218 1218 TIMING_START (gcd_recursion); 1219 G_random_element= 1220 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1219 G_random_element= 1220 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1221 1221 list, topLevel); 1222 TIMING_END_AND_PRINT (gcd_recursion, 1222 TIMING_END_AND_PRINT (gcd_recursion, 1223 1223 "time for recursive call: "); 1224 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 1224 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 1225 1225 } 1226 1226 else if (fail && !inextension) … … 1232 1232 int deg= 2; 1233 1233 do { 1234 mipo= randomIrredpoly (deg, x); 1234 mipo= randomIrredpoly (deg, x); 1235 1235 alpha= rootOf (mipo); 1236 1236 inextension= true; 1237 1237 fail= false; 1238 random_element= randomElement (m, alpha, l, fail); 1238 random_element= randomElement (m, alpha, l, fail); 1239 1239 deg++; 1240 } while (fail); 1240 } while (fail); 1241 1241 list= CFList(); 1242 1242 TIMING_START (gcd_recursion); 1243 G_random_element= 1244 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1243 G_random_element= 1244 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1245 1245 list, topLevel); 1246 TIMING_END_AND_PRINT (gcd_recursion, 1246 TIMING_END_AND_PRINT (gcd_recursion, 1247 1247 "time for recursive call: "); 1248 1248 DEBOUTLN (cerr, "G_random_element= " << G_random_element); … … 1260 1260 CanonicalForm prim_elem, im_prim_elem; 1261 1261 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 1262 1262 1263 1263 ASSERT (!prim_fail, "failure in integer factorizer"); 1264 1264 if (prim_fail) … … 1271 1271 1272 1272 inextensionextension= true; 1273 for (CFListIterator i= l; i.hasItem(); i++) 1274 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 1273 for (CFListIterator i= l; i.hasItem(); i++) 1274 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 1275 1275 im_prim_elem, source, dest); 1276 1276 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1277 1277 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1278 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 1278 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 1279 1279 source, dest); 1280 1280 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1281 1281 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1282 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 1282 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 1283 1283 source, dest); 1284 1284 fail= false; … … 1287 1287 CFList list; 1288 1288 TIMING_START (gcd_recursion); 1289 G_random_element= 1289 G_random_element= 1290 1290 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 1291 1291 list, topLevel); 1292 TIMING_END_AND_PRINT (gcd_recursion, 1292 TIMING_END_AND_PRINT (gcd_recursion, 1293 1293 "time for recursive call: "); 1294 1294 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 1295 1295 alpha= V_buf; 1296 } 1297 1298 d0= totaldegree (G_random_element, Variable(2), 1296 } 1297 1298 d0= totaldegree (G_random_element, Variable(2), 1299 1299 Variable (G_random_element.level())); 1300 1300 … … 1304 1304 return N (reverseSubst (gcdcAcB, substitute)); 1305 1305 else 1306 return N(gcdcAcB); 1307 } 1308 if (d0 > d) 1309 { 1306 return N(gcdcAcB); 1307 } 1308 if (d0 > d) 1309 { 1310 1310 if (!find (l, random_element)) 1311 1311 l.append (random_element); … … 1314 1314 1315 1315 G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element)) 1316 *G_random_element; 1317 1318 1319 d0= totaldegree (G_random_element, Variable(2), 1316 *G_random_element; 1317 1318 1319 d0= totaldegree (G_random_element, Variable(2), 1320 1320 Variable(G_random_element.level())); 1321 1321 1322 if (d0 < d) 1322 if (d0 < d) 1323 1323 { 1324 1324 m= gcdlcAlcB; … … 1327 1327 d= d0; 1328 1328 } 1329 1329 1330 1330 TIMING_START (newton_interpolation); 1331 1331 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); 1332 TIMING_END_AND_PRINT (newton_interpolation, 1332 TIMING_END_AND_PRINT (newton_interpolation, 1333 1333 "time for newton_interpolation: "); 1334 1334 1335 1335 //termination test 1336 if (uni_lcoeff (H) == gcdlcAlcB) 1336 if (uni_lcoeff (H) == gcdlcAlcB) 1337 1337 { 1338 1338 cH= uni_content (H); … … 1340 1340 ppH /= Lc (ppH); 1341 1341 DEBOUTLN (cerr, "ppH= " << ppH); 1342 if (fdivides (ppH, A) && fdivides (ppH, B)) 1342 if (fdivides (ppH, A) && fdivides (ppH, B)) 1343 1343 { 1344 1344 if (substitute > 1) -
factory/cf_gcd_smallp.h
rc840d9 r806c18 5 5 /** @file cf_gcd_smallp.h 6 6 * 7 * @author Martin Lee 7 * @author Martin Lee 8 8 * @date 22.10.2009 9 9 * 10 * This file defines the functions GCD_Fp_extension which computes the GCD of 11 * two polynomials over \f$ F_{p}(\alpha ) \f$ , GCD_small_p which computes the 12 * GCD of two polynomials over \f$ F_{p} \f$ , and GCD_GF which computes the 13 * GCD of two polynomials over GF. Algorithms are based on "On the Genericity of 10 * This file defines the functions GCD_Fp_extension which computes the GCD of 11 * two polynomials over \f$ F_{p}(\alpha ) \f$ , GCD_small_p which computes the 12 * GCD of two polynomials over \f$ F_{p} \f$ , and GCD_GF which computes the 13 * GCD of two polynomials over GF. Algorithms are based on "On the Genericity of 14 14 * the Modular Polynomial GCD Algorithm" by E. Kaltofen & M. Monagan 15 15 * … … 17 17 * (c) by The SINGULAR Team, see LICENSE file 18 18 * 19 * @internal 19 * @internal 20 20 * @version \$Id$ 21 21 * … … 24 24 25 25 #include <config.h> 26 #include "assert.h" 26 #include "assert.h" 27 27 28 28 CanonicalForm GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 29 29 Variable & alpha, CFList& l, bool& top_level); 30 30 31 /// GCD of A and B over \f$ F_{p}(\alpha ) \f$ 32 static inline CanonicalForm GCD_Fp_extension (const CanonicalForm& A, const CanonicalForm& B, 33 Variable & alpha) 31 /// GCD of A and B over \f$ F_{p}(\alpha ) \f$ 32 static inline CanonicalForm GCD_Fp_extension (const CanonicalForm& A, const CanonicalForm& B, 33 Variable & alpha) 34 34 { 35 35 CFList list; … … 54 54 55 55 /// GCD of A and B over GF 56 static inline CanonicalForm GCD_GF (const CanonicalForm& A, const CanonicalForm& B) 56 static inline CanonicalForm GCD_GF (const CanonicalForm& A, const CanonicalForm& B) 57 57 { 58 ASSERT (CFFactory::gettype() == GaloisFieldDomain, 58 ASSERT (CFFactory::gettype() == GaloisFieldDomain, 59 59 "GF as base field expected"); 60 60 CFList list; … … 63 63 } 64 64 65 CanonicalForm 65 CanonicalForm 66 66 randomIrredpoly (int i, const Variable & x) ; 67 67 #endif -
factory/cf_generator.cc
rc840d9 r806c18 55 55 ASSERT( current != gf_q + 1, "no more items" ); 56 56 if ( gf_iszero( current ) ) 57 57 current = 0; 58 58 else if ( current == gf_q1 - 1 ) 59 59 current = gf_q + 1; 60 60 else 61 61 current++; 62 62 } 63 63 … … 137 137 { 138 138 for ( int i = 0; i < n; i++ ) 139 139 result += power( algext, i ) * gensg[i]->item(); 140 140 } 141 141 else 142 142 { 143 143 for ( int i = 0; i < n; i++ ) 144 144 result += power( algext, i ) * gensf[i]->item(); 145 145 } 146 146 return result; … … 156 156 while ( ! stop && i < n ) 157 157 { 158 159 158 gensg[i]->next(); 159 if ( ! gensg[i]->hasItems() ) 160 160 { 161 162 163 164 165 161 gensg[i]->reset(); 162 i++; 163 } 164 else 165 stop = true; 166 166 } 167 167 } … … 170 170 while ( ! stop && i < n ) 171 171 { 172 173 172 gensf[i]->next(); 173 if ( ! gensf[i]->hasItems() ) 174 174 { 175 176 177 178 179 175 gensf[i]->reset(); 176 i++; 177 } 178 else 179 stop = true; 180 180 } 181 181 } 182 182 if ( ! stop ) 183 183 nomoreitems = true; 184 184 } 185 185 … … 188 188 ASSERT( getCharacteristic() > 0, "not a finite field" ); 189 189 if ( getGFDegree() > 1 ) 190 190 return new GFGenerator(); 191 191 else 192 192 return new FFGenerator(); 193 193 } -
factory/cf_inline.cc
rc840d9 r806c18 253 253 { 254 254 if ( (! is_imm( value )) && value->deleteObject() ) 255 255 delete value; 256 256 } 257 257 //}}} … … 273 273 { 274 274 if ( this != &cf ) { 275 276 277 275 if ( (! is_imm( value )) && value->deleteObject() ) 276 delete value; 277 value = (is_imm( cf.value )) ? cf.value : cf.value->copyObject(); 278 278 } 279 279 return *this; … … 306 306 { 307 307 if ( (! is_imm( value )) && value->deleteObject() ) 308 308 delete value; 309 309 value = CFFactory::basic( cf ); 310 310 return *this; … … 351 351 // Use `mpz_cpm_ui()' resp. `mpz_sgn()' to check the underlying 352 352 // mpi. 353 // 353 // 354 354 //}}} 355 355 CF_INLINE bool … … 359 359 360 360 if ( ! what ) 361 361 return value->isOne(); 362 362 else if ( what == INTMARK ) 363 363 return imm_isone( value ); 364 364 else if ( what == FFMARK ) 365 365 return imm_isone_p( value ); 366 366 else 367 367 return imm_isone_gf( value ); 368 368 } 369 369 … … 374 374 375 375 if ( what == 0 ) 376 376 return value->isZero(); 377 377 else if ( what == INTMARK ) 378 378 return imm_iszero( value ); 379 379 else if ( what == FFMARK ) 380 380 return imm_iszero_p( value ); 381 381 else 382 382 return imm_iszero_gf( value ); 383 383 } 384 384 //}}} … … 396 396 // has to be created. 397 397 // 398 // Type info: 398 // Type info: 399 399 // ---------- 400 400 // cf: CurrentPP … … 441 441 442 442 if ( ! what ) 443 443 result.value = result.value->neg(); 444 444 else if ( what == INTMARK ) 445 445 result.value = imm_neg( result.value ); 446 446 else if ( what == FFMARK ) 447 447 result.value = imm_neg_p( result.value ); 448 448 else 449 449 result.value = imm_neg_gf( result.value ); 450 450 451 451 return result; -
factory/cf_irred.cc
rc840d9 r806c18 24 24 int i; 25 25 do { 26 27 28 26 result = power( x, deg ); 27 for ( i = deg-1; i >= 0; i-- ) 28 result += gen.generate() * power( x, i ); 29 29 } while ( ! is_irreducible( result ) ); 30 30 return result; -
factory/cf_iter_inline.cc
rc840d9 r806c18 90 90 ASSERT( hasterms, "lib error: iterator out of terms" ); 91 91 if ( ispoly ) 92 92 return cursor->coeff; 93 93 else 94 94 return data; 95 95 } 96 96 //}}} … … 109 109 ASSERT( hasterms, "lib error: iterator out of terms" ); 110 110 if ( ispoly ) 111 111 return cursor->exp; 112 112 else 113 113 return 0; 114 114 } 115 115 //}}} … … 132 132 ASSERT( hasterms, "lib error: iterator out of terms" ); 133 133 if ( ispoly ) { 134 135 134 cursor = cursor->next; 135 hasterms = cursor != 0; 136 136 } else 137 137 hasterms = false; 138 138 139 139 return *this; … … 145 145 ASSERT( hasterms, "lib error: iterator out of terms" ); 146 146 if ( ispoly ) { 147 148 147 cursor = cursor->next; 148 hasterms = cursor != 0; 149 149 } else 150 150 hasterms = false; 151 151 152 152 return *this; -
factory/cf_map.cc
rc840d9 r806c18 381 381 degsg = degrees( g, degsg ); 382 382 optvalues( degsf, degsg, n, p1, pe ); 383 383 384 384 i = 1; k = 1; 385 385 if ( pe > 1 ) -
factory/cf_map_ext.cc
rc840d9 r806c18 3 3 /** @file cf_map_ext.cc 4 4 * 5 * @author Martin Lee 5 * @author Martin Lee 6 6 * @date 16.11.2009 7 7 * … … 11 11 * (c) by The SINGULAR Team, see LICENSE file 12 12 * 13 * @internal 13 * @internal 14 14 * @version \$Id$ 15 15 * … … 35 35 #include "cf_map_ext.h" 36 36 37 #ifdef HAVE_NTL 37 #ifdef HAVE_NTL 38 38 39 39 /// helper function 40 40 static inline 41 int findItem (const CFList& list, const CanonicalForm& item) 42 { 43 int result= 1; 44 for (CFListIterator i= list; i.hasItem(); i++, result++) 41 int findItem (const CFList& list, const CanonicalForm& item) 42 { 43 int result= 1; 44 for (CFListIterator i= list; i.hasItem(); i++, result++) 45 45 { 46 46 if (i.getItem() == item) … … 52 52 /// helper function 53 53 static inline 54 CanonicalForm getItem (const CFList& list, const int& pos) 54 CanonicalForm getItem (const CFList& list, const int& pos) 55 55 { 56 56 int j= 1; 57 57 if (pos > list.length() || pos < 1) return 0; 58 for (CFListIterator i= list; j <= pos; i++, j++) 59 { 60 if (j == pos) 58 for (CFListIterator i= list; j <= pos; i++, j++) 59 { 60 if (j == pos) 61 61 return i.getItem(); 62 62 } 63 } 64 65 66 /// \f$ F_{p} (\alpha ) \subset F_{p}(\beta ) \f$ and \f$ \alpha \f$ is a 67 /// primitive element, returns the image of \f$ \alpha \f$ 68 static inline 69 CanonicalForm mapUp (const Variable& alpha, const Variable& beta) 63 } 64 65 66 /// \f$ F_{p} (\alpha ) \subset F_{p}(\beta ) \f$ and \f$ \alpha \f$ is a 67 /// primitive element, returns the image of \f$ \alpha \f$ 68 static inline 69 CanonicalForm mapUp (const Variable& alpha, const Variable& beta) 70 70 { 71 71 int p= getCharacteristic (); 72 72 zz_p::init (p); 73 zz_pX NTL_mipo= convertFacCF2NTLzzpX (getMipo (beta)); 73 zz_pX NTL_mipo= convertFacCF2NTLzzpX (getMipo (beta)); 74 74 zz_pE::init (NTL_mipo); 75 75 zz_pEX NTL_alpha_mipo= convertFacCF2NTLzz_pEX (getMipo(alpha), NTL_mipo); … … 78 78 } 79 79 80 /// the CanonicalForm G is the output of map_up, returns F considered as an 80 /// the CanonicalForm G is the output of map_up, returns F considered as an 81 81 /// element over \f$ F_{p}(\alpha ) \f$, WARNING: make sure coefficients of F 82 82 /// are really elements of a subfield of \f$ F_{p}(\beta ) \f$ which is 83 /// isomorphic to \f$ F_{p}(\alpha ) \f$ 84 static inline 85 CanonicalForm 83 /// isomorphic to \f$ F_{p}(\alpha ) \f$ 84 static inline 85 CanonicalForm 86 86 mapDown (const CanonicalForm& F, const Variable& alpha, const 87 CanonicalForm& G, CFList& source, CFList& dest) 88 { 87 CanonicalForm& G, CFList& source, CFList& dest) 88 { 89 89 CanonicalForm buf, buf2; 90 90 int counter= 0; … … 97 97 CanonicalForm alpha_power; 98 98 if (degree(F) == 0) return F; 99 if (F.level() < 0 && F.isUnivariate()) 99 if (F.level() < 0 && F.isUnivariate()) 100 100 { 101 101 buf= F; 102 102 remainder= mod (buf, G); 103 ASSERT (remainder.isZero(), "alpha is not primitive"); 103 ASSERT (remainder.isZero(), "alpha is not primitive"); 104 104 pos= findItem (source, buf); 105 105 if (pos == 0) 106 106 source.append (buf); 107 107 buf2= buf; 108 while (degree (buf) != 0 && counter < bound) 108 while (degree (buf) != 0 && counter < bound) 109 109 { 110 110 buf /= G; 111 111 counter++; 112 112 if (buf == buf2) break; 113 } 114 ASSERT (counter >= bound, "alpha is not primitive"); 115 if (pos == 0) 113 } 114 ASSERT (counter >= bound, "alpha is not primitive"); 115 if (pos == 0) 116 116 { 117 117 alpha_power= power (alpha, counter); … … 119 119 } 120 120 else 121 alpha_power= getItem (dest, pos); 121 alpha_power= getItem (dest, pos); 122 122 result = alpha_power*buf; 123 123 return result; 124 124 } 125 else 126 { 127 for (CFIterator i= F; i.hasTerms(); i++) 125 else 126 { 127 for (CFIterator i= F; i.hasTerms(); i++) 128 128 { 129 129 buf= mapDown (i.coeff(), alpha, G, source, dest); … … 142 142 char gf_name_buf= gf_name; 143 143 InternalCF* buf; 144 if (F.inBaseDomain()) 144 if (F.inBaseDomain()) 145 145 { 146 146 if (F.isOne()) return 1; … … 149 149 result= power (alpha, exp).mapinto(); 150 150 return result; 151 } 151 } 152 152 for (CFIterator i= F; i.hasTerms(); i++) 153 153 result += GF2FalphaHelper (i.coeff(), alpha)*power (F.mvar(), i.exp()); … … 155 155 } 156 156 157 /// changes representation by primitive element to representation by residue 158 /// classes modulo a Conway polynomial 159 CanonicalForm GF2FalphaRep (const CanonicalForm& F, const Variable& alpha) 157 /// changes representation by primitive element to representation by residue 158 /// classes modulo a Conway polynomial 159 CanonicalForm GF2FalphaRep (const CanonicalForm& F, const Variable& alpha) 160 160 { 161 161 Variable beta= rootOf (gf_mipo); … … 163 163 } 164 164 165 /// change representation by residue classes modulo a Conway polynomial 165 /// change representation by residue classes modulo a Conway polynomial 166 166 /// to representation by primitive element 167 CanonicalForm Falpha2GFRep (const CanonicalForm& F) 167 CanonicalForm Falpha2GFRep (const CanonicalForm& F) 168 168 { 169 169 CanonicalForm result= 0; 170 170 InternalCF* buf; 171 171 172 if (F.inCoeffDomain()) 172 if (F.inCoeffDomain()) 173 173 { 174 174 if (F.inBaseDomain()) 175 175 return F.mapinto(); 176 else 177 { 178 for (CFIterator i= F; i.hasTerms(); i++) 176 else 177 { 178 for (CFIterator i= F; i.hasTerms(); i++) 179 179 { 180 180 buf= int2imm_gf (i.exp()); … … 183 183 } 184 184 return result; 185 } 186 for (CFIterator i= F; i.hasTerms(); i++) 185 } 186 for (CFIterator i= F; i.hasTerms(); i++) 187 187 result += Falpha2GFRep (i.coeff())*power (F.mvar(), i.exp()); 188 188 return result; … … 191 191 /// GF_map_up helper 192 192 static inline 193 CanonicalForm GFPowUp (const CanonicalForm & F, int k) 193 CanonicalForm GFPowUp (const CanonicalForm & F, int k) 194 194 { 195 195 if (F.isOne()) return F; 196 196 CanonicalForm result= 0; 197 if (F.inBaseDomain()) 197 if (F.inBaseDomain()) 198 198 return power(F, k); 199 for (CFIterator i= F; i.hasTerms(); i++) 199 for (CFIterator i= F; i.hasTerms(); i++) 200 200 result += GFPowUp (i.coeff(), k)*power (F.mvar(), i.exp()); 201 201 return result; 202 202 } 203 203 204 /// maps a polynomial over \f$ GF(p^{k}) \f$ to a polynomial over 204 /// maps a polynomial over \f$ GF(p^{k}) \f$ to a polynomial over 205 205 /// \f$ GF(p^{d}) \f$ , d needs to be a multiple of k 206 CanonicalForm GFMapUp (const CanonicalForm & F, int k) 207 { 206 CanonicalForm GFMapUp (const CanonicalForm & F, int k) 207 { 208 208 int d= getGFDegree(); 209 209 ASSERT (d%k == 0, "multiple of GF degree expected"); … … 217 217 /// GFMapDown helper 218 218 static inline 219 CanonicalForm GFPowDown (const CanonicalForm & F, int k) 219 CanonicalForm GFPowDown (const CanonicalForm & F, int k) 220 220 { 221 221 if (F.isOne()) return F; … … 223 223 int exp; 224 224 InternalCF* buf; 225 if (F.inBaseDomain()) 225 if (F.inBaseDomain()) 226 226 { 227 227 buf= F.getval(); … … 229 229 if ((exp % k) == 0) 230 230 exp= exp/k; 231 else 231 else 232 232 return -1; 233 233 234 234 buf= int2imm_gf (exp); 235 235 return CanonicalForm (buf); 236 } 237 for (CFIterator i= F; i.hasTerms(); i++) 236 } 237 for (CFIterator i= F; i.hasTerms(); i++) 238 238 result += GFPowDown (i.coeff(), k)*power (F.mvar(), i.exp()); 239 239 return result; 240 240 } 241 241 242 /// maps a polynomial over \f$ GF(p^{d}) \f$ to a polynomial over 242 /// maps a polynomial over \f$ GF(p^{d}) \f$ to a polynomial over 243 243 /// \f$ GF(p^{k})\f$ , d needs to be a multiple of k 244 CanonicalForm GFMapDown (const CanonicalForm & F, int k) 244 CanonicalForm GFMapDown (const CanonicalForm & F, int k) 245 245 { 246 246 int d= getGFDegree(); … … 253 253 } 254 254 255 static inline 256 CanonicalForm mapUp (const CanonicalForm& F, const CanonicalForm& G, 257 const Variable& alpha, const CanonicalForm& H, 255 static inline 256 CanonicalForm mapUp (const CanonicalForm& F, const CanonicalForm& G, 257 const Variable& alpha, const CanonicalForm& H, 258 258 CFList& source, CFList& dest) 259 { 259 { 260 260 CanonicalForm buf, buf2; 261 261 int counter= 0; … … 268 268 CanonicalForm H_power; 269 269 if (degree(F) <= 0) return F; 270 if (F.level() < 0 && F.isUnivariate()) 270 if (F.level() < 0 && F.isUnivariate()) 271 271 { 272 272 buf= F; 273 273 remainder= mod (buf, G); 274 ASSERT (remainder.isZero(), "alpha is not primitive"); 274 ASSERT (remainder.isZero(), "alpha is not primitive"); 275 275 pos= findItem (source, buf); 276 276 if (pos == 0) 277 277 source.append (buf); 278 278 buf2= buf; 279 while (degree (buf) != 0 && counter < bound) 279 while (degree (buf) != 0 && counter < bound) 280 280 { 281 281 buf /= G; 282 282 counter++; 283 283 if (buf == buf2) break; 284 } 285 ASSERT (counter >= bound, "alpha is not primitive"); 286 if (pos == 0) 284 } 285 ASSERT (counter >= bound, "alpha is not primitive"); 286 if (pos == 0) 287 287 { 288 288 H_power= power (H, counter); … … 290 290 } 291 291 else 292 H_power= getItem (dest, pos); 292 H_power= getItem (dest, pos); 293 293 result = H_power*buf; 294 294 return result; 295 295 } 296 else 297 { 298 for (CFIterator i= F; i.hasTerms(); i++) 296 else 297 { 298 for (CFIterator i= F; i.hasTerms(); i++) 299 299 { 300 300 buf= mapUp (i.coeff(), G, alpha, H, source, dest); … … 305 305 } 306 306 307 /// determine a primitive element of \f$ F_{p} (\alpha ) \f$, 308 /// \f$ \beta \f$ is a primitive element of a field which is isomorphic to 307 /// determine a primitive element of \f$ F_{p} (\alpha ) \f$, 308 /// \f$ \beta \f$ is a primitive element of a field which is isomorphic to 309 309 /// \f$ F_{p}(\alpha ) \f$ 310 CanonicalForm 311 primitiveElement (const Variable& alpha, Variable& beta, bool fail) 310 CanonicalForm 311 primitiveElement (const Variable& alpha, Variable& beta, bool fail) 312 312 { 313 313 bool primitive= false; … … 316 316 if (fail) 317 317 return 0; 318 if (primitive) 318 if (primitive) 319 319 { 320 320 beta= alpha; 321 321 return alpha; 322 322 } 323 CanonicalForm mipo= getMipo (alpha); 323 CanonicalForm mipo= getMipo (alpha); 324 324 int d= degree (mipo); 325 325 int p= getCharacteristic (); … … 331 331 do 332 332 { 333 BuildIrred (NTL_mipo, d); 333 BuildIrred (NTL_mipo, d); 334 334 mipo2= convertNTLzzpX2CF (NTL_mipo, Variable (1)); 335 335 beta= rootOf (mipo2); 336 336 primitive= isPrimitive (beta, fail); 337 337 if (primitive) 338 break; 338 break; 339 339 if (fail) 340 340 return 0; … … 348 348 CanonicalForm 349 349 mapDown (const CanonicalForm& F, const CanonicalForm& prim_elem, const 350 CanonicalForm& im_prim_elem, const Variable& alpha, CFList& source, 351 CFList& dest) 352 { 350 CanonicalForm& im_prim_elem, const Variable& alpha, CFList& source, 351 CFList& dest) 352 { 353 353 return mapUp (F, im_prim_elem, alpha, prim_elem, dest, source); 354 354 } 355 355 356 356 CanonicalForm 357 mapUp (const CanonicalForm& F, const Variable& alpha, const Variable& beta, 358 const CanonicalForm& prim_elem, const CanonicalForm& im_prim_elem, 359 CFList& source, CFList& dest) 360 { 361 if (prim_elem == alpha) 357 mapUp (const CanonicalForm& F, const Variable& alpha, const Variable& beta, 358 const CanonicalForm& prim_elem, const CanonicalForm& im_prim_elem, 359 CFList& source, CFList& dest) 360 { 361 if (prim_elem == alpha) 362 362 return F (im_prim_elem, alpha); 363 363 return mapUp (F, prim_elem, alpha, im_prim_elem, source, dest); 364 364 } 365 365 366 CanonicalForm 367 mapPrimElem (const CanonicalForm& prim_elem, const Variable& alpha, 366 CanonicalForm 367 mapPrimElem (const CanonicalForm& prim_elem, const Variable& alpha, 368 368 const Variable& beta) 369 369 { … … 377 377 result += power (im_alpha, i.exp())*i.coeff(); 378 378 return result; 379 } 379 } 380 380 } 381 381 -
factory/cf_map_ext.h
rc840d9 r806c18 5 5 /** @file cf_map_ext.h 6 6 * 7 * @author Martin Lee 7 * @author Martin Lee 8 8 * @date 16.11.2009 9 9 * … … 13 13 * (c) by The SINGULAR Team, see LICENSE file 14 14 * 15 * @internal 15 * @internal 16 16 * @version \$Id$ 17 17 * … … 26 26 27 27 CanonicalForm 28 mapUp (const CanonicalForm& F, const Variable& alpha, const Variable& beta, 29 const CanonicalForm& prim_elem, const CanonicalForm& im_prim_elem, 28 mapUp (const CanonicalForm& F, const Variable& alpha, const Variable& beta, 29 const CanonicalForm& prim_elem, const CanonicalForm& im_prim_elem, 30 30 CFList& source, CFList& dest); 31 31 32 32 CanonicalForm 33 33 mapDown (const CanonicalForm& F, const CanonicalForm& prim_elem, const 34 CanonicalForm& im_prim_elem, const Variable& alpha, CFList& source, 35 CFList& dest); 34 CanonicalForm& im_prim_elem, const Variable& alpha, CFList& source, 35 CFList& dest); 36 36 37 CanonicalForm 37 CanonicalForm 38 38 primitiveElement (const Variable& alpha, Variable& beta, bool fail); 39 39 40 CanonicalForm 41 mapPrimElem (const CanonicalForm& prim_elem, const Variable& alpha, 40 CanonicalForm 41 mapPrimElem (const CanonicalForm& prim_elem, const Variable& alpha, 42 42 const Variable& beta); 43 43 -
factory/cf_primes.cc
rc840d9 r806c18 15 15 ASSERT( i >= 0 && i < NUMPRIMES, "index to primes too high" ); 16 16 if ( i >= NUMSMALLPRIMES ) 17 17 return bigprimes[i-NUMSMALLPRIMES]; 18 18 else 19 19 return smallprimes[i]; 20 20 } 21 21 -
factory/cf_random.cc
rc840d9 r806c18 144 144 CanonicalForm result; 145 145 for ( int i = 0; i < n; i++ ) 146 146 result += power( algext, i ) * gen->generate(); 147 147 return result; 148 148 }; … … 156 156 { 157 157 if ( getCharacteristic() == 0 ) 158 158 return new IntRandom(); 159 159 if ( getGFDegree() > 1 ) 160 160 return new GFRandom(); 161 161 else 162 162 return new FFRandom(); 163 163 } 164 164 … … 166 166 { 167 167 if ( n == 0 ) 168 168 return (int)ranGen.generate(); 169 169 else 170 170 return ranGen.generate() % n; 171 171 } 172 172 -
factory/cf_resultant.cc
rc840d9 r806c18 48 48 // some checks on triviality 49 49 if ( f.isZero() || g.isZero() ) { 50 51 50 trivialResult[0] = 0; 51 return trivialResult; 52 52 } 53 53 54 54 // make x main variable 55 55 if ( f.mvar() > x || g.mvar() > x ) { 56 57 58 59 60 61 56 if ( f.mvar() > g.mvar() ) 57 X = f.mvar(); 58 else 59 X = g.mvar(); 60 F = swapvar( f, X, x ); 61 G = swapvar( g, X, x ); 62 62 } 63 63 else { 64 65 66 64 X = x; 65 F = f; 66 G = g; 67 67 } 68 68 // at this point, we have to calculate the sequence of F and … … 83 83 // make sure that S[j+1] is regular and j < n 84 84 if ( m == n && j > 0 ) { 85 86 85 S[j-1] = LC( S[j], X ) * psr( S[j+1], S[j], X ); 86 j--; 87 87 } else if ( m < n ) { 88 89 88 S[j-1] = LC( S[j], X ) * LC( S[j], X ) * S[j+1]; 89 j--; 90 90 } else if ( m > n && j > 0 ) { 91 92 93 94 95 96 97 98 99 100 101 102 103 91 // calculate first step 92 r = degree( S[j], X ); 93 R = LC( S[j+1], X ); 94 95 // if there was a gap calculate similar polynomial 96 if ( j > r && r >= 0 ) 97 S[r] = power( LC( S[j], X ), j - r ) * S[j] * power( R, j - r ); 98 99 if ( r > 0 ) { 100 // calculate remainder 101 S[r-1] = psr( S[j+1], S[j], X ) * power( -R, j - r ); 102 j = r-1; 103 } 104 104 } 105 105 106 106 while ( j > 0 ) { 107 108 109 110 111 112 113 114 115 116 117 118 119 120 107 // at this point, 0 < j < n and S[j+1] is regular 108 r = degree( S[j], X ); 109 R = LC( S[j+1], X ); 110 111 // if there was a gap calculate similar polynomial 112 if ( j > r && r >= 0 ) 113 S[r] = (power( LC( S[j], X ), j - r ) * S[j]) / power( R, j - r ); 114 115 if ( r <= 0 ) break; 116 // calculate remainder 117 S[r-1] = psr( S[j+1], S[j], X ) / power( -R, j - r + 2 ); 118 119 j = r-1; 120 // again 0 <= j < r <= jOld and S[j+1] is regular 121 121 } 122 122 123 123 for ( j = 0; j <= S.max(); j++ ) { 124 125 126 127 124 // reswap variables if necessary 125 if ( X != x ) { 126 S[j] = swapvar( S[j], X, x ); 127 } 128 128 } 129 129 … … 148 148 // f or g in R 149 149 if ( degree( f, x ) == 0 ) 150 150 return power( f, degree( g, x ) ); 151 151 if ( degree( g, x ) == 0 ) 152 152 return power( g, degree( f, x ) ); 153 153 154 154 // f and g are linear polynomials … … 180 180 // here because this may involve variable swapping. 181 181 if ( f.isZero() || g.isZero() ) 182 182 return 0; 183 183 if ( f.mvar() < x ) 184 184 return power( f, g.degree( x ) ); 185 185 if ( g.mvar() < x ) 186 186 return power( g, f.degree( x ) ); 187 187 188 188 // make x main variale … … 190 190 Variable X; 191 191 if ( f.mvar() > x || g.mvar() > x ) { 192 193 194 195 196 197 192 if ( f.mvar() > g.mvar() ) 193 X = f.mvar(); 194 else 195 X = g.mvar(); 196 F = swapvar( f, X, x ); 197 G = swapvar( g, X, x ); 198 198 } 199 199 else { 200 201 202 200 X = x; 201 F = f; 202 G = g; 203 203 } 204 204 // at this point, we have to calculate resultant( F, G, X ) … … 210 210 // catch trivial cases 211 211 if ( m+n <= 2 || m == 0 || n == 0 ) 212 212 return swapvar( trivialResultant( F, G, X ), X, x ); 213 213 214 214 // exchange F and G if necessary 215 215 int flipFactor; 216 216 if ( m < n ) { 217 218 219 220 221 222 223 224 217 CanonicalForm swap = F; 218 F = G; G = swap; 219 int degswap = m; 220 m = n; n = degswap; 221 if ( m & 1 && n & 1 ) 222 flipFactor = -1; 223 else 224 flipFactor = 1; 225 225 } else 226 226 flipFactor = 1; 227 227 228 228 // this is not an effective way to calculate the resultant! 229 229 CanonicalForm extFactor; 230 230 if ( m == n ) { 231 232 233 234 231 if ( n & 1 ) 232 extFactor = -LC( G, X ); 233 else 234 extFactor = LC( G, X ); 235 235 } else 236 236 extFactor = power( LC( F, X ), m-n-1 ); 237 237 238 238 CanonicalForm result; -
factory/cf_reval.cc
rc840d9 r806c18 15 15 { 16 16 if ( e.gen == 0 ) 17 17 gen = 0; 18 18 else 19 19 gen = e.gen->clone(); 20 20 values = e.values; 21 21 } … … 24 24 { 25 25 if ( gen != 0 ) 26 26 delete gen; 27 27 } 28 28 … … 32 32 if ( this != &e ) { 33 33 if (gen!=0) 34 35 36 37 38 39 34 delete gen; 35 values = e.values; 36 if ( e.gen == 0 ) 37 gen = 0; 38 else 39 gen = e.gen->clone(); 40 40 } 41 41 return *this; … … 47 47 int n = values.max(); 48 48 for ( int i = values.min(); i <= n; i++ ) 49 49 values[i] = gen->generate(); 50 50 } -
factory/cf_switches.cc
rc840d9 r806c18 26 26 { 27 27 for ( int i = 0; i < CFSwitchesMax; i++ ) 28 28 switches[i] = false; 29 29 // and set the default (recommended) On-values: 30 30 #ifdef HAVE_NTL -
factory/debug.cc
rc840d9 r806c18 11 11 // deb_level == -1 iff we enter this function for the first time 12 12 if ( deb_level == -1 ) 13 13 deb_level = 0; 14 14 else 15 15 delete [] deb_level_msg; 16 16 17 17 deb_level++; 18 18 deb_level_msg = new char[3*deb_level+1]; 19 19 for ( i = 0; i < 3*deb_level; i++ ) 20 20 deb_level_msg[i] = ' '; 21 21 deb_level_msg[3*deb_level] = '\0'; 22 22 } … … 25 25 { 26 26 if ( deb_level > 0 ) { 27 28 29 30 31 32 33 27 int i; 28 deb_level--; 29 delete [] deb_level_msg; 30 deb_level_msg = new char[3*deb_level+1]; 31 for ( i = 0; i < 3*deb_level; i++ ) 32 deb_level_msg[i] = ' '; 33 deb_level_msg[3*deb_level] = '\0'; 34 34 } 35 35 } -
factory/facAlgExt.cc
rc840d9 r806c18 4 4 * 5 5 * @author Martin Lee 6 * @date 6 * @date 7 7 * 8 8 * Univariate factorization over algebraic extension of Q using Trager's … … 12 12 * (c) by The SINGULAR Team, see LICENSE file 13 13 * 14 * @internal 14 * @internal 15 15 * @version \$Id$ 16 16 * … … 38 38 CanonicalForm G= deriv (F, F.mvar()); 39 39 G= gcd (F, G); 40 return F/G; 40 return F/G; 41 41 } 42 42 … … 49 49 CanonicalForm mipo= getMipo (alpha); 50 50 mipo= mipo (x, alpha); 51 51 52 52 CanonicalForm norm= resultant (g, mipo, x); 53 53 i= 0; … … 88 88 ASSERT (F.isUnivariate(), "univariate input expected"); 89 89 ASSERT (getCharacteristic() == 0, "characteristic 0 expected"); 90 90 91 91 On (SW_RATIONAL); 92 CanonicalForm f= F; 92 CanonicalForm f= F; 93 93 int shift; 94 94 CanonicalForm norm= sqrfNorm (f, alpha, shift); … … 98 98 if (normFactors.length() <= 2) 99 99 return CFList (f); 100 100 101 101 normFactors.removeFirst(); 102 102 CanonicalForm buf; … … 104 104 buf= f (f.mvar() - shift*alpha, f.mvar()); 105 105 else 106 buf= f; 106 buf= f; 107 107 CanonicalForm factor; 108 108 for (CFFListIterator i= normFactors; i.hasItem(); i++) … … 112 112 buf /= factor; 113 113 if (shift != 0) 114 factor= factor (f.mvar() + shift*alpha, f.mvar()); 114 factor= factor (f.mvar() + shift*alpha, f.mvar()); 115 115 factors.append (factor); 116 } 116 } 117 117 ASSERT (degree (buf) <= 0, "incomplete factorization"); 118 118 Off (SW_RATIONAL); … … 120 120 } 121 121 122 CFFList 122 CFFList 123 123 AlgExtFactorize (const CanonicalForm& F, const Variable& alpha) 124 124 { 125 125 ASSERT (F.isUnivariate(), "univariate input expected"); 126 126 ASSERT (getCharacteristic() == 0, "characteristic 0 expected"); 127 127 128 128 if (F.inCoeffDomain()) 129 129 return CFFList (CFFactor (F, 1)); … … 135 135 CanonicalForm buf= F/Lc (F); 136 136 CFFList factors; 137 int multi; 137 int multi; 138 138 for (CFListIterator i= sqrfFactors; i.hasItem(); i++) 139 139 { … … 146 146 } 147 147 factors.append (CFFactor (i.getItem(), multi)); 148 } 148 } 149 149 Off (SW_RATIONAL); 150 150 factors.insert (CFFactor (Lc(F), 1)); 151 151 ASSERT (degree (buf) <= 0, "bug in AlgExtFactorize"); 152 return factors; 152 return factors; 153 153 } 154 154 -
factory/facAlgExt.h
rc840d9 r806c18 4 4 * 5 5 * @author Martin Lee 6 * @date 6 * @date 7 7 * 8 8 * Univariate factorization over algebraic extension of Q using Trager's … … 12 12 * (c) by The SINGULAR Team, see LICENSE file 13 13 * 14 * @internal 14 * @internal 15 15 * @version \$Id$ 16 16 * … … 24 24 ///factorize a univariate squarefree polynomial over algebraic extension of Q 25 25 /// 26 /// @return @a AlgExtSqrfFactorize returns a list of factors of F 27 CFList 26 /// @return @a AlgExtSqrfFactorize returns a list of factors of F 27 CFList 28 28 AlgExtSqrfFactorize (const CanonicalForm& F, ///<[in] a univariate squarefree 29 29 ///< polynomial … … 33 33 /// factorize a univariate polynomial over algebraic extension of Q 34 34 /// 35 /// @return @a AlgExtFactorize returns a list of factors of F with multiplicity 36 CFFList 35 /// @return @a AlgExtFactorize returns a list of factors of F with multiplicity 36 CFFList 37 37 AlgExtFactorize (const CanonicalForm& F, ///<[in] a univariate polynomial 38 38 const Variable& alpha ///<[in] an algebraic variable -
factory/facFqBivar.cc
rc840d9 r806c18 1 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 2 * Computer Algebra System SINGULAR 3 3 \*****************************************************************************/ 4 4 /** @file facFqBivar.cc 5 * 5 * 6 6 * This file provides functions for factorizing a bivariate polynomial over 7 * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer 8 * Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring 9 * multivariate polynomials over a finite field" by L. Bernardin. 7 * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer 8 * Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring 9 * multivariate polynomials over a finite field" by L. Bernardin. 10 10 * 11 * ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and 12 * factorize the polynomial in both variables. So far factor recombination is 13 * done naive! 11 * ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and 12 * factorize the polynomial in both variables. So far factor recombination is 13 * done naive! 14 14 * 15 15 * @author Martin Lee … … 43 43 TIMING_DEFINE_PRINT(fac_factor_recombination); 44 44 45 CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M) 45 CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M) 46 46 { 47 47 if (L.isEmpty()) … … 57 57 CFList tmp1, tmp2; 58 58 CanonicalForm buf1, buf2; 59 for (int j= 1; j <= l; j++, i++) 59 for (int j= 1; j <= l; j++, i++) 60 60 tmp1.append (i.getItem()); 61 61 tmp2= Difference (L, tmp1); … … 66 66 } 67 67