1 // -*- mode:C++ ; compile-command: "g++ -I.. -g -c modpoly.cc" -*- 2 /* 3 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 19 #ifndef _GIAC_MODPOLY_H_ 20 #define _GIAC_MODPOLY_H_ 21 #include "first.h" 22 #include "global.h" 23 #include "fraction.h" 24 #ifdef HAVE_LIBNTL 25 #include <NTL/ZZXFactoring.h> 26 #include <NTL/ZZ.h> 27 #include <NTL/ZZ_p.h> 28 #include <NTL/GF2X.h> 29 #include <NTL/pair_GF2X_long.h> 30 #include <NTL/GF2XFactoring.h> 31 #include <NTL/ZZ_pX.h> 32 #include <NTL/ZZX.h> 33 #include <NTL/pair_ZZ_pX_long.h> 34 #include <NTL/ZZ_pXFactoring.h> 35 #ifdef HAVE_PTHREAD_H 36 #include <pthread.h> 37 #endif 38 #endif // HAVE_LIBNTL 39 40 #ifndef NO_NAMESPACE_GIAC 41 namespace giac { 42 #endif // ndef NO_NAMESPACE_GIAC 43 // Previous Fourier prime if fourier_for_n!=0, or previous prime <= maxp 44 int prevprimep1p2p3(int p,int maxp,int fourier_for_n=0); 45 template<class T> class tensor; 46 typedef tensor<gen> polynome; 47 typedef vecteur modpoly; 48 typedef vecteur dense_POLY1; // same internal rep but assumes non modular op 49 50 // set env->pn, returns true if size of field <=2^31 51 bool normalize_env(environment * env); 52 // modular or GF inverse 53 gen invenv(const gen & g,environment * env); 54 55 // *********************** 56 // Building and converting 57 // *********************** 58 // Conversion from univariate to multivariate polynomials 59 modpoly modularize(const polynome & p,const gen & n,environment * env); 60 modpoly modularize(const dense_POLY1 & p,const gen & n,environment * env); 61 polynome unmodularize(const modpoly & a); 62 bool is_one(const modpoly & p); 63 int is_cyclotomic(const modpoly & p,GIAC_CONTEXT); 64 int is_cyclotomic(const modpoly & p,double eps); 65 modpoly one(); // 1 66 modpoly xpower1(); // x 67 modpoly xpowerpn(environment * env); // x^modulo 68 vecteur x_to_xp(const vecteur & v, int p); // x -> x^p non modular 69 gen nrandom(environment * env); // random modular integer 70 modpoly random(int i,environment * env); // random univariate polynomial of degree i 71 void shiftmodpoly(modpoly & a,int n); // multiply by x^n 72 // high = high*x^n + low, size of low must be < n 73 void mergemodpoly(modpoly & high,const modpoly & low,int n); 74 modpoly trim(const modpoly & p,environment * env); 75 void trim_inplace(modpoly & p); 76 void fast_trim_inplace(std::vector<int> & p,int modulo,int maxsize=-1); 77 bool trim(modpoly & v); // true if v is empty after trimming 78 void rrdm(modpoly & p, int n); // right redimension poly to degree n 79 80 // *************** 81 // Arithmetic 82 // *************** 83 // !! Do not call Addmodpoly with modpoly slices if new_coord and th/other overlapp 84 void Addmodpoly(modpoly::const_iterator th_it,modpoly::const_iterator th_itend,modpoly::const_iterator other_it,modpoly::const_iterator other_itend,environment * env, modpoly & new_coord); 85 void addmodpoly(const modpoly & th, const modpoly & other, modpoly & new_coord); 86 void addmodpoly(const modpoly & th, const modpoly & other, environment * env,modpoly & new_coord); 87 modpoly operator_plus (const modpoly & th,const modpoly & other,environment * env); 88 modpoly operator + (const modpoly & th,const modpoly & other); 89 void Submodpoly(modpoly::const_iterator th_it,modpoly::const_iterator th_itend,modpoly::const_iterator other_it,modpoly::const_iterator other_itend,environment * env,modpoly & new_coord); 90 void submodpoly(const modpoly & th, const modpoly & other, environment * env,modpoly & new_coord); 91 void submodpoly(const modpoly & th, const modpoly & other, modpoly & new_coord); 92 modpoly operator_minus (const modpoly & th,const modpoly & other,environment * env); 93 modpoly operator - (const modpoly & th,const modpoly & other); 94 void negmodpoly(const modpoly & th, modpoly & new_coord); 95 modpoly operator - (const modpoly & th) ; 96 97 // Warning: mulmodpoly assumes that coeff are integers. 98 // Use operator_times unless you know what you do 99 void mulmodpoly(const modpoly & th, const gen & fact, modpoly & new_coord); 100 void mulmodpoly(const modpoly & th, const gen & fact, environment * env, modpoly & new_coord); 101 void mulmodpoly_naive(modpoly::const_iterator ita,modpoly::const_iterator ita_end,modpoly::const_iterator itb,modpoly::const_iterator itb_end,environment * env,modpoly & new_coord); 102 void mulmodpoly_kara_naive(const modpoly & a, const modpoly & b,environment * env,modpoly & new_coord,int seuil_kara); 103 // new_coord += a*b in place 104 void add_mulmodpoly(const modpoly::const_iterator & ita0,const modpoly::const_iterator & ita_end,const modpoly::const_iterator & itb0,const modpoly::const_iterator & itb_end,environment * env,modpoly & new_coord); 105 // modpoly operator * (const modpoly & th, const gen & fact); 106 modpoly operator_times (const modpoly & th, const gen & fact,environment * env); 107 // commented otherwise int * gen might be interpreted as 108 // make a modpoly of size the int and multiply 109 modpoly operator_times (const gen & fact,const modpoly & th,environment * env); 110 modpoly operator * (const gen & fact,const modpoly & th); 111 void mulmodpoly(const modpoly & a, const modpoly & b, environment * env,modpoly & new_coord,int maxdeg=RAND_MAX); 112 modpoly operator * (const modpoly & th, const modpoly & other) ; 113 modpoly operator_times (const modpoly & th, const modpoly & other,environment * env) ; 114 void operator_times (const modpoly & a, const modpoly & b,environment * env,modpoly & new_coord,int maxdeg=RAND_MAX); 115 // res=(*it) * ... (*(it_end-1)) 116 void mulmodpoly(std::vector<modpoly>::const_iterator it,std::vector<modpoly>::const_iterator it_end,environment * env,modpoly & new_coord); 117 void mulmodpoly(std::vector<modpoly>::const_iterator * it,int debut,int fin,environment * env,modpoly & pi); 118 119 void divmodpoly(const modpoly & th, const gen & fact, modpoly & new_coord); 120 void divmodpoly(const modpoly & th, const gen & fact, environment * env,modpoly & new_coord); 121 modpoly operator / (const modpoly & th,const modpoly & other) ; 122 int coefftype(const modpoly & v,gen & coefft); 123 bool DivRem(const modpoly & th, const modpoly & other, environment * env,modpoly & quo, modpoly & rem,bool allowrational=true); 124 bool DenseDivRem(const modpoly & th, const modpoly & other,modpoly & quo, modpoly & rem,bool fastfalsetest=false); 125 // Pseudo division a*th = other*quo + rem 126 void PseudoDivRem(const dense_POLY1 & th, const dense_POLY1 & other, dense_POLY1 & quo, dense_POLY1 & rem, gen & a); 127 modpoly operator / (const modpoly & th,const gen & fact ) ; 128 modpoly operator_div (const modpoly & th,const modpoly & other,environment * env) ; 129 modpoly operator % (const modpoly & th,const modpoly & other) ; 130 modpoly operator_mod (const modpoly & th,const modpoly & other,environment * env) ; 131 // division by ascending power (used for power series) 132 dense_POLY1 AscPowDivRem(const dense_POLY1 & num, const dense_POLY1 & den,int order); 133 134 modpoly powmod(const modpoly & p,const gen & n,const modpoly & pmod,environment * env); 135 gen cstcoeff(const modpoly & q); 136 // p=(X-x)q+p(x), first horner returns p(x), second one compute q as well 137 gen horner(const gen & g,const gen & x); 138 gen horner(const gen & g,const gen & x); 139 gen hornermod(const vecteur & v,const gen & alpha,const gen & modulo); 140 int hornermod(const vecteur & v,int alpha,int modulo); 141 142 extern const unary_function_ptr * const at_horner ; 143 gen horner(const modpoly & p,const gen & x,environment * env,bool simp=true); 144 gen horner(const modpoly & p,const gen & x); 145 gen horner(const modpoly & p,const gen & x,environment * env,modpoly & q); 146 gen horner(const modpoly & p,const fraction & f,bool simp); 147 gen _horner(const gen & args,GIAC_CONTEXT); 148 std::complex<double> horner_newton(const vecteur & p,const std::complex<double> &x,GIAC_CONTEXT); // x-p(x)/p'(x) 149 150 void hornerfrac(const modpoly & p,const gen &num, const gen &den,gen & res,gen & d); // res/d=p(num/den) 151 // find bounds for p(interval[l,r]) with p, l and r real and exact 152 vecteur horner_interval(const modpoly & p,const gen & l,const gen & r); 153 154 gen symb_horner(const modpoly & p,const gen & x,int d); 155 gen symb_horner(const modpoly & p,const gen & x); 156 // shift polynomial 157 modpoly taylor(const modpoly & p,const gen & x,environment * env=0); 158 // split P=Pp-Pn in two parts, Pp positive coeffs and Pn negative coeffs 159 void splitP(const vecteur &P,vecteur &Pp,vecteur &Pn); 160 161 gen norm(const dense_POLY1 & q,GIAC_CONTEXT); // max of |coeff| 162 modpoly derivative(const modpoly & p); 163 modpoly derivative(const modpoly & p,environment * env); 164 // integration of p with coeff in *ascending* order, does not add 0 165 // if you use it with usual modpoly, you must reverse(p.begin(),p.end()) 166 // before and after usage, set shift_coeff to 1, and push_back a zero 167 modpoly integrate(const modpoly & p,const gen & shift_coeff); 168 169 // *************** 170 // GCD and related. Works currently only if modular aritmetic 171 // *************** 172 modpoly simplify(modpoly & a, modpoly & b,environment * env); 173 gen ppz(dense_POLY1 & q); 174 gen lgcd(const dense_POLY1 & p); // gcd of coeffs 175 gen lgcd(const dense_POLY1 & p,const gen & g); 176 // modular gcd 177 void modpoly2smallmodpoly(const modpoly & p,std::vector<int> & v,int m); 178 179 bool gcdmodpoly(const modpoly &p,const modpoly & q,environment * env,modpoly &a); 180 bool gcd_modular_algo(const modpoly &p,const modpoly &q,modpoly &d,modpoly * p_simp,modpoly * q_simp); // p and q must have coeffs in Z or Z[i] 181 182 // half-gcd: a0.size() must be > a1.size(), returns [[A,B],[C,D]] 183 bool hgcd(const modpoly & a0,const modpoly & a1,const gen & modulo,modpoly &A,modpoly &B,modpoly &C,modpoly &D); // a0 is A in Yap, a1 is B 184 // fast modular inverse: f*g=1 mod x^l, f must be invertible (f.back()!=0) 185 bool invmod(const modpoly & f,int l,environment * env,modpoly & g); 186 // for p prime such that p-1 is divisible by 2^N, compute a 2^N-th root of 1 187 // otherwise return 0 188 unsigned nthroot(unsigned p,unsigned N); 189 // euclidean quotient using modular inverse 190 // returns 0 on failure, 1 on success, 2 if division is exact 191 int DivQuo(const modpoly & a, const modpoly & b, environment * env,modpoly & q); 192 // 1-d modular for small modulus<sqrt(RAND_MAX) 193 bool gcdsmallmodpoly(const polynome &p,const polynome & q,int m,polynome & d,polynome & dp,polynome & dq,bool compute_cof); 194 void smallmodpoly2modpoly(const std::vector<int> & v,modpoly & p,int m); 195 void gcdsmallmodpoly(const std::vector<int> &p,const std::vector<int> & q,int m,std::vector<int> & d); 196 void gcdsmallmodpoly(const std::vector<int> &p,const std::vector<int> & q,int m,std::vector<int> & d,std::vector<int> * pcof,std::vector<int> * qcof); 197 void gcdsmallmodpoly(const std::vector<int> &p,const std::vector<int> & q,int m,std::vector<int> & d); 198 void gcdsmallmodpoly(const modpoly &p,const modpoly & q,int m,modpoly & d); 199 void DivRem(const std::vector<int> & th, const std::vector<int> & other,int m,std::vector<int> & quo, std::vector<int> & rem,bool ck_exactquo=false); 200 modpoly gcd(const modpoly & a,const modpoly &b,environment * env,bool call_ntl=false); 201 // n-var modular gcd 202 bool gcd_modular(const polynome &p_orig, const polynome & q_orig, polynome & pgcd,polynome & pcofactor,polynome & qcofactor,bool compute_cofactors); 203 204 bool convert(const polynome &p_orig, const polynome & q_orig,index_t & d,std::vector<hashgcd_U> & vars,std::vector< T_unsigned<gen,hashgcd_U> > & p,std::vector< T_unsigned<gen,hashgcd_U> > & q); 205 bool convert(const polynome &p_orig, const polynome & q_orig,index_t & d,std::vector<hashgcd_U> & vars,std::vector< T_unsigned<int,hashgcd_U> > & p,std::vector< T_unsigned<int,hashgcd_U> > & q,int modulo); 206 bool modgcd(const polynome &p_orig, const polynome & q_orig, const gen & modulo, polynome & d,polynome & pcofactor,polynome & qcofactor,bool compute_cofactors); 207 bool mod_gcd_c(const polynome &p_orig, const polynome & q_orig, const gen & modulo, polynome & d,polynome & pcofactor,polynome & qcofactor,bool compute_cofactors); 208 bool mod_gcd(const polynome &p_orig, const polynome & q_orig, const gen & modulo, polynome & pgcd,polynome & pcofactor,polynome & qcofactor,bool compute_cofactors); 209 modpoly lcm(const modpoly & a,const modpoly &b,environment * env); 210 bool gcd_modular_algo1(polynome &p,polynome &q,polynome &d,bool compute_cof); 211 212 // check if a (normally a fraction) modulo m is equal to p 213 bool chk_equal_mod(const gen & a,longlong p,int m); 214 // p1*u+p2*v=d 215 void egcd(const modpoly &p1, const modpoly & p2, environment * env,modpoly & u,modpoly & v,modpoly & d,bool deterministic=true); 216 // a*u+b*v=d optimized for GMP, if u_ptr and v_ptr are non 0 compute previous line of Euclide algorithm, stops as soon as degree<degstop 217 bool egcd_mpz(const modpoly & a,const modpoly &b,int degstop,const gen & m,modpoly & u,modpoly &v,modpoly & d,modpoly * u_ptr,modpoly * v_ptr,modpoly * r_ptr); 218 bool egcd_pade(const modpoly & n,const modpoly & x,int l,modpoly & a,modpoly &b,environment * env,bool psron=true); 219 // alg extension norme for p_y, alg extension defined by pmini 220 bool algnorme(const polynome & p_y,const polynome & pmini,polynome & n); 221 void subresultant(const modpoly & P,const modpoly & Q,gen & res); 222 int sizeinbase2(const gen & g); 223 int sizeinbase2(const vecteur & v); 224 // multinomial power by Miller pure recurrence 225 bool miller_pow(const modpoly & p_,unsigned m,modpoly & res); 226 227 // Given [v_0 ... v_(2n-1)] (begin of the recurrence sequence) 228 // return [b_n...b_0] such that b_n*v_{n+k}+...+b_0*v_k=0 229 // Example [1,-1,3,3] -> [1,-3,-6] 230 vecteur reverse_rsolve(const vecteur & v,bool psron=true); 231 // given a and c, find u such that 232 // a[0]*...a[n-1]*u[n]+a[0]*...*a[n-2]*a[n]*u[n-1]+...+a[1]*...*a[n-1]*u[0]=1 233 bool egcd(const std::vector<modpoly> & a, environment * env,std::vector<modpoly> & u); 234 // same as above 235 std::vector<modpoly> egcd(const std::vector<modpoly> & a, environment * env); 236 // assuming pmod and qmod are prime together, find r such that 237 // r = p mod pmod and r = q mod qmod 238 // hence r = p + A*pmod = q + B*qmod 239 // or A*pmod -B*qmod = q - p 240 // assuming u*pmod+v*pmod=d we get 241 // A=u*(q-p)/d 242 dense_POLY1 ichinrem(const dense_POLY1 &p,const dense_POLY1 & q,const gen & pmod,const gen & qmod); 243 modpoly chinrem(const modpoly & p,const modpoly & q, const modpoly & pmod, const modpoly & qmod,environment * env); 244 int ichinrem_inplace(dense_POLY1 &p,const std::vector<int> & q,const gen & pmod,int qmodval,int reserve_mem=0); // 0 error, 1 p changed, 2 p unchangedb 245 bool ichinrem_inplace(dense_POLY1 &p,const dense_POLY1 & q,const gen & pmod,int qmodval); 246 void divided_differences(const vecteur & x,const vecteur & y,vecteur & res,environment * env); 247 // in-place modification and exact division if divexact==true 248 void divided_differences(const vecteur & x,vecteur & res,environment * env,bool divexact); 249 void interpolate(const vecteur & x,const vecteur & y,modpoly & res,environment * env); 250 void interpolate_inplace(const vecteur & x,modpoly & res,environment * env); 251 void mulpoly_interpolate(const polynome & p,const polynome & q,polynome & res,environment * env); 252 bool dotvecteur_interp(const vecteur & a,const vecteur &b,gen & res); 253 bool mmult_interp(const matrice & a,const matrice &b,matrice & res); 254 bool poly_pcar_interp(const matrice & a,vecteur & p,bool compute_pmin,GIAC_CONTEXT); 255 void polymat2matpoly(const vecteur & R,vecteur & res); 256 257 void makepositive(int * p,int n,int modulo); 258 // Fast Fourier Transform, f the poly sum_{j<n} f_j x^j, 259 // and w=[1,omega,...,omega^[m-1]] with m a multiple of n 260 // return [f(1),f(omega),...,f(omega^[n-1]) 261 // WARNING f is given in ascending power 262 void fft(const modpoly & f,const modpoly & w ,modpoly & res,environment * env); 263 // void fft(std::vector< std::complex<double> >& f,const std::vector< std::complex<double> > & w ,std::vector<std::complex< double> > & res); 264 void fft(std::complex<double> * f,int n,const std::complex<double> * w,int m,std::complex< double> * t); 265 void fft(const std::vector<int> & f,const std::vector<int> & w ,std::vector<int> & res,int modulo); 266 // res=a .* b mod p 267 void fft_ab_p(const std::vector<int> &a,const std::vector<int> &b,std::vector<int> & res,int p); 268 // res=a ./ b mod p 269 void fft_aoverb_p(const std::vector<int> &a,const std::vector<int> &b,std::vector<int> & res,int p); 270 // reverse the table of root of unity^k for inverse fft 271 void fft_reverse(std::vector<int> & W,int p); 272 273 // res=a*b mod p 274 bool fft2mult(int ablinfnorm,const std::vector<int> & a,const std::vector<int> & b,std::vector<int> & res,int modulo,std::vector<int> & W,std::vector<int> & fftmult_p,std::vector<int> & fftmult_q,bool reverseatend,bool dividebyn,bool makeplus); 275 bool fftmultp1234(const modpoly & p,const modpoly & q,const gen &P,const gen &Q,modpoly & pq,int modulo, std::vector<int> & a,std::vector<int>&b,std::vector<int> &resp1,std::vector<int>&resp2,std::vector<int> & resp3, std::vector<int> & Wp1,std::vector<int> & Wp2, std::vector<int> & Wp3,std::vector<int> & Wp4,std::vector<int> &tmp_p,std::vector<int> &tmp_q,bool compute_pq); 276 // FFT mod 2^{r*2^l}+1, tmp1, tmp2 temporary gen must be _ZINT 277 void fft2rl(gen * f,long n,int r,int l,gen * t,bool direct,gen & tmp1, gen & tmp2,mpz_t & tmpqz); 278 // alpha[i] *= beta[i] mod 2^(expoN)+1 279 void fft2rltimes(modpoly & alpha,const modpoly & beta,unsigned long expoN,mpz_t & tmp,mpz_t & tmpqz); 280 void fft2rltimes(const modpoly & alpha,const modpoly & beta,modpoly & res,unsigned long expoN,mpz_t & tmp,mpz_t & tmpqz); 281 // pq *= -2^shift mod N=2^(expoN+1) where -2^shift is the inverse of n mod N 282 void fft2rldiv(modpoly & pq,unsigned long expoN,unsigned long shift,mpz_t & tmp,mpz_t & tmpqz); 283 gen intnorm(const dense_POLY1 & p,GIAC_CONTEXT); 284 285 // FFT representation is a triplet of FFT mod p1, p2, p3 286 struct fft_rep { 287 int modulo; 288 // modulo=0 means we make the computation for the 3 primes p1, p2 and p3 289 // modulo!=0 means we reconstruct mod modulo with p1, p2, p3 290 std::vector<int> modp1,modp2,modp3; 291 }; 292 293 // p -> f it's FFT representation, p is reversed before FFT is called 294 // a is a temporary vector = A mod modulo after call 295 // Wp1, Wp2, Wp3 is a vector of powers of the n-th root of unity mod p1,p2,p3 296 // if size is not n or Wp[0]==0, Wp is computed 297 // do not share Wp1/p2/p3 between different threads 298 void to_fft(const vecteur & A,int modulo,std::vector<int> & Wp1,std::vector<int> & Wp2,std::vector<int> & Wp3,std::vector<int> & a,int n,fft_rep & f,bool reverse,bool makeplus); 299 void to_fft(const std::vector<int> & a,int modulo,int w,std::vector<int> & Wp,int n,std::vector<int> & f,int reverse,bool makeplus,bool makemod=true); 300 void to_fft(const std::vector<int> & a,int modulo,std::vector<int> & Wp1,std::vector<int> & Wp2,std::vector<int> & Wp3,int n,fft_rep & f,bool reverse,bool makeplus,bool makemod=true); 301 // FFT representation f -> res 302 // Wp1,p2,p3 should be computed with to_fft 303 // division by n=size of f.modp1/p2/p3 is done 304 // result should normally be reversed at end 305 // tmp1/p2/p3 are temporary vectors 306 // do not share Wp1/p2/p3 between different threads 307 void from_fft(const fft_rep & f,std::vector<int> & Wp1,std::vector<int> & Wp2,std::vector<int> & Wp3,std::vector<int> & res,std::vector<int> & tmp1,std::vector<int> & tmp2,std::vector<int> & tmp3,bool reverseatend,bool revw); 308 void from_fft(const std::vector<int> & f,int p,std::vector<int> & Wp,std::vector<int> & res,bool reverseatend,bool revw); 309 310 struct multi_fft_rep { 311 gen modulo; 312 fft_rep p1p2p3; // for p1, p2 and p3 313 std::vector<fft_rep> v; // one fft_rep for each prime < p1, !=p2 and !=p3 314 }; 315 // convert A to multi fft representation, with a number of primes 316 // suitable for ichinrem reconstruction mod modulo 317 void to_multi_fft(const vecteur & A,const gen & modulo,std::vector<int> & Wp1,std::vector<int> & Wp2,std::vector<int> & Wp3,unsigned long n,multi_fft_rep & f,bool reverse,bool makeplus); 318 void from_multi_fft(const multi_fft_rep & f,std::vector<int> & Wp1,std::vector<int> & Wp2,std::vector<int> & Wp3,vecteur & res,bool reverseatend); 319 void multi_fft_ab_cd(const multi_fft_rep & a,const multi_fft_rep & b,const multi_fft_rep & c,const multi_fft_rep & d,multi_fft_rep & res); 320 321 // Convolution of p and q, omega a n-th root of unity, n=2^k 322 // WARNING p0 and q0 are given in ascending power 323 void fftconv(const modpoly & p,const modpoly & q,unsigned long k,unsigned long n,const gen & omega,modpoly & pq,environment * env); 324 // Convolution of p and q, omega a n-th root of unity, n=2^k 325 // p and q are given in descending power order 326 void fftconv(const modpoly & p0,const modpoly & q0,unsigned long k,const gen & omega,modpoly & pq,environment * env); 327 bool fftmult(const modpoly & p,const modpoly & q,modpoly & pq,int modulo,int maxdeg=RAND_MAX); 328 modpoly fftmult(const modpoly & p,const modpoly & q); 329 // input A with positive int, output fft in A 330 // w a 2^n-th root of unity mod p 331 void fft2( int *A, int n, int w, int p ,bool permute=true); 332 // float fft, theta should be +/-2*M_PI/n 333 void fft2( std::complex<double> * A, int n, double theta ); 334 // resultant on Z, if eps!=0 might be non deterministic 335 gen mod_resultant(const modpoly & P,const modpoly & Q,double eps); 336 modpoly unmod(const modpoly & a,const gen & m); 337 // resultant of P and Q modulo m, modifies P and Q, 338 int resultant_int(std::vector<int> & P,std::vector<int> & Q,std::vector<int> & tmp1,std::vector<int> & tmp2,int m,int w=0); 339 void vecteur2vector_ll(const vecteur & v,longlong m,std::vector<longlong> & res); 340 longlong resultantll(std::vector<longlong> & P,std::vector<longlong> & Q,std::vector<longlong> & tmp1,std::vector<longlong> & tmp2,longlong m); 341 342 bool ntlresultant(const modpoly &p,const modpoly &q,const gen & modulo,gen & res,bool ntl_on_check=true); 343 bool ntlxgcd(const modpoly &a,const modpoly &b,const gen & modulo,modpoly & reu,modpoly &v,modpoly & d,bool ntl_on_check=true); 344 bool ntlgcd(const modpoly &a,const modpoly &b,const gen & modulo,modpoly & d,bool ntl_on_check=true); 345 #ifdef HAVE_LIBNTL 346 #ifdef HAVE_LIBPTHREAD 347 extern pthread_mutex_t ntl_mutex; 348 #endif 349 typedef gen inttype; 350 351 bool polynome2tab(const polynome & p,int deg,inttype * tab); 352 polynome tab2polynome(const inttype * tab,int deg); 353 bool ininttype2ZZ(const inttype & temp,const inttype & step,NTL::ZZ & z,const NTL::ZZ & zzstep); 354 NTL::ZZ inttype2ZZ(const inttype & i); 355 void inZZ2inttype(const NTL::ZZ & zztemp,int pow2,inttype & temp); 356 inttype ZZ2inttype(const NTL::ZZ & z); 357 NTL::ZZX tab2ZZX(const inttype * tab,int degree); 358 void ZZX2tab(const NTL::ZZX & f,int & degree,inttype * & tab); 359 360 // Z/nZ[X] conversions 361 NTL::GF2X modpoly2GF2X(const modpoly & p); 362 modpoly GF2X2modpoly(const NTL::GF2X & f); 363 NTL::ZZ_pX modpoly2ZZ_pX(const modpoly & p); 364 modpoly ZZ_pX2modpoly(const NTL::ZZ_pX & f); 365 366 #endif 367 368 #ifndef NO_NAMESPACE_GIAC 369 } // namespace giac 370 #endif // NO_NAMESPACE_GIAC 371 372 #endif // _GIAC_MODPOLY_H 373