/* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -I../include -g -c threaded.cc -DHAVE_CONFIG_H -DIN_GIAC -D_I386_" -*- */ #include "giacPCH.h" /* Copyright (C) 2000,2007 B. Parisse, Institut Fourier, 38402 St Martin d'Heres * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ using namespace std; #include "threaded.h" #include "sym2poly.h" #include "gausspol.h" #include "usual.h" #include "monomial.h" #include "modpoly.h" #include "giacintl.h" #include "input_parser.h" #ifndef NO_NAMESPACE_GIAC namespace giac { #endif // ndef NO_NAMESPACE_GIAC #ifdef HAVE_GMPXX_H mpz_class invmod(const mpz_class & a,int reduce){ mpz_class z=(a%reduce); int tmp=z.get_si(); tmp=giac::invmod(tmp,reduce); return tmp; } mpz_class smod(const mpz_class & a,int reduce){ mpz_class z=(a%reduce); int tmp=z.get_si(); tmp=giac::smod(tmp,reduce); return tmp; } #endif double heap_mult=20000; gen _heap_mult(const gen & g0,GIAC_CONTEXT){ if ( g0.type==_STRNG && g0.subtype==-1) return g0; gen g=evalf_double(g0,1,contextptr); if (g.type!=_DOUBLE_) return heap_mult; return heap_mult=g._DOUBLE_val; } static const char _heap_mult_s []="heap_mult"; static define_unary_function_eval (__heap_mult,&_heap_mult,_heap_mult_s); define_unary_function_ptr5( at_heap_mult ,alias_at_heap_mult,&__heap_mult,0,true); double modgcd_cachesize=6291456; gen _modgcd_cachesize(const gen & g0,GIAC_CONTEXT){ if ( g0.type==_STRNG && g0.subtype==-1) return g0; gen g=evalf_double(g0,1,contextptr); if (g.type!=_DOUBLE_) return modgcd_cachesize; return modgcd_cachesize=g._DOUBLE_val; } static const char _modgcd_cachesize_s []="modgcd_cachesize"; static define_unary_function_eval (__modgcd_cachesize,&_modgcd_cachesize,_modgcd_cachesize_s); define_unary_function_ptr5( at_modgcd_cachesize ,alias_at_modgcd_cachesize,&__modgcd_cachesize,0,true); gen _debug_infolevel(const gen & g0,GIAC_CONTEXT){ if ( g0.type==_STRNG && g0.subtype==-1) return g0; gen g=evalf_double(g0,1,contextptr); if (g.type!=_DOUBLE_) return debug_infolevel; return debug_infolevel=int(g._DOUBLE_val); } static const char _debug_infolevel_s []="debug_infolevel"; static define_unary_function_eval (__debug_infolevel,&_debug_infolevel,_debug_infolevel_s); define_unary_function_ptr5( at_debug_infolevel ,alias_at_debug_infolevel,&__debug_infolevel,0,T_DIGITS); gen _step_infolevel(const gen & g0,GIAC_CONTEXT){ if ( g0.type==_STRNG && g0.subtype==-1) return g0; gen g=evalf_double(g0,1,contextptr); if (g.type!=_DOUBLE_) return step_infolevel(contextptr); return step_infolevel(contextptr)=int(g._DOUBLE_val); } static const char _step_infolevel_s []="step_infolevel"; static define_unary_function_eval (__step_infolevel,&_step_infolevel,_step_infolevel_s); define_unary_function_ptr5( at_step_infolevel ,alias_at_step_infolevel,&__step_infolevel,0,true); my_mpz operator % (const my_mpz & a,const my_mpz & b){ my_mpz tmp; mpz_fdiv_r(tmp.ptr,a.ptr,b.ptr); return tmp; } my_mpz operator %= (my_mpz & a,const my_mpz & b){ mpz_fdiv_r(a.ptr,a.ptr,b.ptr); return a; } my_mpz operator += (my_mpz & a,const my_mpz & b){ mpz_add(a.ptr,a.ptr,b.ptr); return a; } my_mpz operator -= (my_mpz & a,const my_mpz & b){ mpz_sub(a.ptr,a.ptr,b.ptr); return a; } my_mpz operator + (const my_mpz & a,const my_mpz & b){ my_mpz tmp; mpz_add(tmp.ptr,a.ptr,b.ptr); return tmp; } my_mpz operator - (const my_mpz & a,const my_mpz & b){ my_mpz tmp; mpz_sub(tmp.ptr,a.ptr,b.ptr); return tmp; } my_mpz operator * (const my_mpz & a,const my_mpz & b){ my_mpz tmp; mpz_mul(tmp.ptr,a.ptr,b.ptr); return tmp; } my_mpz operator / (const my_mpz & a,const my_mpz & b){ my_mpz tmp; mpz_fdiv_q(tmp.ptr,a.ptr,b.ptr); return tmp; } my_mpz invmod(const my_mpz & a,int reduce){ my_mpz z=(a%reduce); int tmp=mpz_get_si(z.ptr); tmp=invmod(tmp,reduce); return tmp; } my_mpz smod(const my_mpz & a,int reduce){ my_mpz z=(a%reduce); int tmp=mpz_get_si(z.ptr); tmp=smod(tmp,reduce); return tmp; } void wait_1ms(context * contextptr){ #ifdef NSPIRE sleep(1); #else usleep(1000); #endif } void background_callback(const gen & g,void * newcontextptr){ if (g.type==_VECT && g._VECTptr->size()==2){ context * cptr=(giac::context *)newcontextptr; if (cptr){ #ifdef HAVE_LIBPTHREAD pthread_mutex_lock(cptr->globalptr->_mutex_eval_status_ptr); sto(g._VECTptr->back(),g._VECTptr->front(),cptr); pthread_mutex_unlock(cptr->globalptr->_mutex_eval_status_ptr); #endif } } } gen _background(const gen & g,GIAC_CONTEXT){ if ( g.type==_STRNG && g.subtype==-1) return g; #ifdef HAVE_LIBPTHREAD if (g.type!=_VECT || g._VECTptr->size()<2) return gensizeerr(contextptr); vecteur v(*g._VECTptr); gen target=v[0]; gen toeval=v[1]; int s=v.size(); double maxsize=1e9; // in bytes double maxtime=1e9; // in microseconds int level=eval_level(contextptr); if (s>2){ gen tmp=evalf_double(v[2],level,contextptr); if (tmp.type!=_DOUBLE_) return gentypeerr(contextptr); maxsize=tmp._DOUBLE_val; } if (s>3){ gen tmp=evalf_double(v[3],level,contextptr); if (tmp.type!=_DOUBLE_) return gentypeerr(contextptr); maxtime=tmp._DOUBLE_val; } if (s>4 && v[4].type==_INT_) level=v[4].val; gen tmp; context * newcontextptr=clone_context(contextptr); newcontextptr->parent=contextptr; tmp=gen(newcontextptr,_THREAD_POINTER); sto(tmp,target,contextptr); if (!make_thread(makevecteur(symbolic(at_quote,target),toeval),level,background_callback,(void *)newcontextptr,newcontextptr)){ sto(undef,target,contextptr); return gensizeerr(gettext("Unable to make thread")); } return tmp; #else return undef; #endif } static const char _background_s []="background"; static define_unary_function_eval_quoted (__background,&_background,_background_s); define_unary_function_ptr5( at_background ,alias_at_background,&__background,_QUOTE_ARGUMENTS,true); // check if var is a power of 2 static int find_shift(const hashgcd_U & var){ hashgcd_U u(var); int shift=-1; for (;u;u = u>> 1){ ++shift; } return (int(var) == (1< & vars,index_t & shift){ shift.clear(); std::vector::const_iterator it=vars.begin(),itend=vars.end(); shift.reserve(itend-it); for (;it!=itend;++it){ shift.push_back(find_shift(*it)); if (shift.back()==-1) return false; } return true; } static bool gcdsmallmodpoly(const vector & a,const vector & b,const vector * pminptr,int modulo,vector & d){ gcdsmallmodpoly(a,b,modulo,d); return true; } static bool DivRem(const vector & a,const vector & b,const vector * pminptr,int modulo,vector & q,vector & r){ DivRem(a,b,modulo,q,r); return true; } static bool divrem_(vector< vector > & a,vector< vector > & b,const vector & pmin, int modulo, vector< vector > * qptr,bool set_q_orig_b,vector & b0,vector & b0inv,vector & tmp,vector & tmp1,vector & tmp2,vector & tmp3,vector & tmp4,vector & tmp5); static bool DivRem(const vector< vector > & a,const vector< vector > & b,const vector * pminptr,int modulo,vector< vector > & q,vector< vector > & r){ r=a; vector< vector > b0(b); vector B0,B0inv,tmp,tmp1,tmp2,tmp3,tmp4,tmp5; return divrem_(r,b0,*pminptr,modulo,&q,true,B0,B0inv,tmp,tmp1,tmp2,tmp3,tmp4,tmp5); } static bool gcdsmallmodpoly_ext(const vector< vector > & p,const vector< vector > & q,const vector & pmin,int modulo,vector< vector > & d); static bool gcdsmallmodpoly(const vector< vector > & a,const vector< vector > & b,const vector * pminptr,int modulo,vector< vector > & d){ return gcdsmallmodpoly_ext(a,b,*pminptr,modulo,d); } bool is_zero(const vector & v){ vector::const_iterator it=v.begin(),itend=v.end(); for (;it!=itend;++it){ if (*it) return false; } return true; } inline int make_unit(int ){ return 1; } inline vector make_unit(const vector){ vector v(1,1); return v; } // extract modular content of p wrt to the last variable template static bool pp_mod(vector< T_unsigned > & p,const vector * pminptr,int modulo,hashgcd_U var,hashgcd_U var2,vector & pcontxn){ typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); pcontxn.clear(); if (it==itend) return true; hashgcd_U u=(p.front().u/var)*var,curu,newu=0; if (u==p.front().u){ pcontxn.push_back(make_unit(p.front().g)); return true; } vector current,tmp,reste; for (;it!=itend;){ current.clear(); current.push_back(it->g); curu=it->u; ++it; for (;it!=itend;++it){ newu=it->u; if ( newu < u ){ break; } if (curu>newu+var2) current.insert(current.end(),(curu-newu)/var2-1,T(0)); current.push_back(it->g); curu=newu; } if (curu>u) current.insert(current.end(),(curu-u)/var2,T(0)); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || !gcdsmallmodpoly(pcontxn,current,pminptr,modulo,tmp)) return false; pcontxn=tmp; if (tmp.size()==1) return true; u=(newu/var)*var; } vector< T_unsigned > res; it=p.begin(); u=(p.front().u/var)*var; for (;it!=itend;){ current.clear(); current.push_back(it->g); curu=it->u; ++it; for (;it!=itend;++it){ newu=it->u; if ( newu < u ){ break; } if (curu>newu+var2) current.insert(current.end(),(curu-newu)/var2-1,T(0)); current.push_back(it->g); curu=newu; } if (curu>u) current.insert(current.end(),(curu-u)/var2,T(0)); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || !DivRem(current,pcontxn,pminptr,modulo,tmp,reste)) return false; typename vector::const_iterator jt=tmp.begin(),jtend=tmp.end(); for (int s=int(jtend-jt)-1;jt!=jtend;++jt,--s){ if (!is_zero(*jt)) res.push_back(T_unsigned(*jt,u+s*var2)); } u=(newu/var)*var; } p=res; return true; } template static bool operator < (const vector< T_unsigned > & v1,const vector< T_unsigned > & v2){ return v1.size() > & p,hashgcd_U var,hashgcd_U var2){ vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); hashgcd_U degxn=0; hashgcd_U u; for (;it!=itend;++it){ u=it->u%var; if (degxn static unsigned degree_xn(const vector< T_unsigned > & p,short int shift_var,short int shift_var2){ typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(),it1; unsigned degxn=0; unsigned u,uend; for (;it!=itend;++it){ uend = ((it->u >> shift_var) << shift_var); u = (it->u - uend) >> shift_var2 ; if (!u) continue; if (degxn=itend-it) continue; it1 = it+u; if (it1->u==uend) it = it1; } return degxn ; } /* template int degree(const vector< T_unsigned > & p, const std::vector & vars, index_t & res){ typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); std::vector::const_iterator jtbeg=vars.begin(),jtend=vars.end(),jt; hashgcd_U u,ur; res=index_t(jtend-jtbeg); index_t::iterator ktbeg=res.begin(),ktend=res.end(),kt; int totaldeg=0,current; for (;it!=itend;++it){ u=it->u; current=0; for (jt=jtbeg,kt=ktbeg;jt!=jtend;++jt,++kt){ ur=u%(*jt); u=u/(*jt); if (u>*kt) *kt=u; current += u; u=ur; } if (current>totaldeg) totaldeg=current; } return totaldeg; } */ template static int degree(const vector< T_unsigned > & p, const index_t & shift_vars, index_t & res){ typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); index_t::const_iterator jtbeg=shift_vars.begin(),jtend=shift_vars.end(),jt; int dim=int(jtend-jtbeg); if (dim==1){ int deg=it->u >> *jtbeg; res.clear(); res.push_back(deg); return deg; } hashgcd_U u,uq,uend,skip; short int lastvar=shift_vars[shift_vars.size()-2],shift_var2=shift_vars.back(); res=index_t(jtend-jtbeg); index_t::iterator ktbeg=res.begin(),ktend=res.end(),kt; int totaldeg=0,current; for (;it!=itend;){ u=it->u; // find degree for *it and skip monomials having x1..xn-1 in common uend = (u >> lastvar) << lastvar; skip = (u-uend) >> shift_var2; current=0; // find partial and total degree for (jt=jtbeg,kt=ktbeg;jt!=jtend;++jt,++kt){ uq = u >> *jt; u -= uq << *jt; if ((int)uq>*kt) *kt = uq; current += uq; } if (current>totaldeg) totaldeg=current; if ((int)skipu==uend){ it += skip; ++it; continue; } for (++it;it!=itend;++it){ if (it->u,hashgcd_U> > & p_orig,const vector< T_unsigned,hashgcd_U> > & q_orig,const std::vector & vars,const vector & pmin,int modulo, vector< T_unsigned,hashgcd_U> > & d, vector< T_unsigned,hashgcd_U> > & pcof,vector< T_unsigned,hashgcd_U> > & qcof,bool compute_pcof,bool compute_qcof, int nthreads); bool mod_gcd(const vector< T_unsigned,hashgcd_U> > & p_orig,const vector< T_unsigned,hashgcd_U> > & q_orig,const vector * pminptr,int modulo,const std::vector & vars, vector< T_unsigned,hashgcd_U> > & d, vector< T_unsigned,hashgcd_U> > & pcof,vector< T_unsigned,hashgcd_U> > & qcof,bool compute_cof, int nthreads){ return mod_gcd_ext(p_orig,q_orig,vars,*pminptr,modulo,d,pcof,qcof,compute_cof,compute_cof,nthreads)!=0; } bool mod_gcd(const vector< T_unsigned > & p_orig,const vector< T_unsigned > & q_orig,const vector * pminptr,int modulo,const std::vector & vars, vector< T_unsigned > & d, vector< T_unsigned > & pcofactor, vector< T_unsigned > & qcofactor,bool compute_cofactor,int nthreads){ return mod_gcd(p_orig,q_orig,modulo,d,pcofactor,qcofactor,vars,compute_cofactor,nthreads); } struct modred { int modulo; vector pmin; modred(int m,const vector & p):modulo(m),pmin(p) {} modred():modulo(0),pmin(0) {} }; static bool is_zero(const modred & r){ return r.modulo==0; } static int make_modulo(const vector * pminptr,int modulo,int){ return modulo; } static modred make_modulo(const vector * pminptr,int modulo,vector){ return modred(modulo,*pminptr); } struct Modred { int modulo; vecteur pmin; Modred(int m,const vecteur & p):modulo(m),pmin(p) {} Modred():modulo(0),pmin(0) {} }; static bool is_zero(const Modred & r){ return r.modulo==0; } // extract modular content of p wrt to all but the last variable template static bool pp_mod(vector< T_unsigned > & p,const vector * pminptr,int modulo,const std::vector & vars,vector< T_unsigned > & pcont,int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return false; #else typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); pcont.clear(); pcont.push_back(T_unsigned(make_unit(T(0)),0)); if (it==itend || vars.empty()) return true; std::vector varsn(vars); varsn.pop_back(); // hashgcd_U degxnu=0; hashgcd_U u,u0,var=varsn.back(),var2=vars.back(); short int shiftvar=find_shift(var),shiftvar2=find_shift(var2); unsigned degxn = degree_xn(p,shiftvar,shiftvar2); vector nterms(degxn+1); vector nonzero(degxn+1,false); // number of non constant terms int pos; for (it=p.begin();it!=itend;++it){ u = it->u; pos = (u%var) >> shiftvar2; nonzero[pos]=true; if (u>>shiftvar) ++nterms[pos]; } vector< vector< T_unsigned > > vp(degxn+1); for (int i=degxn;i>=0;--i){ if (!nterms[i] && nonzero[i]) // all terms would be constant in vp[i] return true; vp[i].reserve(nterms[i]+1); } for (it=p.begin();it!=itend;++it){ u = it->u; u0 = (u >> shiftvar) << shiftvar; vp[(u-u0)>>shiftvar2].push_back(T_unsigned(it->g, u0)); } // sort vp by size sort(vp.begin(),vp.end()); // compute gcd unsigned int i=0; for (;i<=degxn;++i){ if (!vp[i].empty()){ pcont=vp[i]; break; } } if (pcont.empty()) CERR << "empty" << endl; vector< T_unsigned > res,rem,pcof,qcof; for (++i;i<=degxn;++i){ if (pcont.size()==1 && pcont.front().u==0) return true; #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || !mod_gcd(pcont,vp[i],pminptr,modulo,varsn,pcont,pcof,qcof,false,nthreads)) return false; } if (pcont.size()==1 && pcont.front().u==0) return true; hashdivrem(p,pcont,res,rem,vars,make_modulo(pminptr,modulo,T(0)),0,false); swap(p,res); return true; #endif // NO_TEMPLATE_MULTGCD } // fast check if p is primitive with respect to the main var // p main var is y, inner var is x static bool is_front_primitive(vector< vector > & p,int modulo){ vector< vector >::iterator it=p.begin(),itend=p.end(); int degy=int(itend-it)-1; vector degrees; int degs=0,d; for (it=p.begin();it!=itend;++it,--degy){ d=int(it->size()); // there is a x^d*y^degy term -> set degree of x^0 to x^d to degy at least if (d>degs){ degrees.insert(degrees.end(),d-degs,degy); degs=d; } } vector::iterator jt=degrees.begin(),jtend=degrees.end(); for (;jt!=jtend;++jt){ if (!*jt) return true; } return false; } // eval p at x with respect to the last variable // keep only terms of degree <= maxdeg with respect to the main variable static bool horner(const vector< T_unsigned > & p,int x,const std::vector & vars,vector< T_unsigned > & px,int modulo,int maxdeg){ hashgcd_U var=vars[vars.size()-2]; hashgcd_U var2=vars.back(); vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(),it1,it2; vector::const_iterator jtend=vars.end()-1; hashgcd_U ucur,uend; if (maxdeg>=0){ uend=(maxdeg+1)*vars.front(); // dichotomy to find start position int pos1=0,pos2=int(itend-it),pos; for (;pos2-pos1>1;){ pos=(pos1+pos2)/2; if ((it+pos)->uu>=uend) ++it; } if (x==0){ px.clear(); for (;it!=itend;){ ucur=it->u; uend=(ucur/var)*var; if (ucur==uend){ register int g=smod(it->g,modulo); if (g!=0) px.push_back(T_unsigned(g,uend)); ++it; continue; } register int nterms = (ucur-uend)/var2; if (ntermsu==uend){ it += nterms; register int g=smod(it->g,modulo); if (g!=0) px.push_back(T_unsigned(g,uend)); ++it; continue; } for (++it;it!=itend;++it){ if (it->u<=uend){ if (it->u==uend){ register int g=smod(it->g,modulo); if (g!=0) px.push_back(T_unsigned(g,uend)); ++it; } break; } } } return true; } vector< T_unsigned >::iterator kt=px.begin(),ktend=px.end(); for (;it!=itend;){ ucur=it->u; uend=(ucur/var)*var; if (ucur==uend){ if (kt!=ktend){ *kt=*it; ++kt; } else px.push_back(*it); ++it; continue; } int g=0; int nterms=(ucur-uend)/var2+1; if (x==1 && nterms < RAND_MAX/modulo ){ for (;it!=itend;++it){ if (it->u(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } break; } else g += it->g; } if (it==itend){ g=smod(g,modulo); if (g!=0){ if (kt!=ktend){ *kt=T_unsigned(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } } continue; } else { // Check if the next group of monomials is dense wrt to xn it1=it+nterms; if (//false && ntermsu==uend ){ if (modulo<=46340){ if (x>=-14 && x<=14){ if (x>=-8 && x<=8){ it2=it+5*((it1-it)/5); for (;it!=it2;){ g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; g %= modulo; ++it; } } else { it2=it+4*((it1-it)/4); for (;it!=it2;){ g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; g %= modulo; ++it; } } } else { if (x>=-35 && x<=35){ it2=it+3*((it1-it)/3); for (;it!=it2;){ g *= x; g += it->g; ++it; g *= x; g += it->g; ++it; g *= x; g += it->g; g %= modulo; ++it; } } } for (;it!=it1;++it){ g = (g*x+it->g)%modulo; } } // end if (modulo<46430) else { // modulo>=46430, using longlong if (x>=-84 && x<=84){ longlong G=g; it2=it+5*((it1-it)/5); for (;it!=it2;){ G *= x; G += it->g; ++it; G *= x; G += it->g; ++it; G *= x; G += it->g; ++it; G *= x; G += it->g; ++it; G *= x; G += it->g; G %= modulo; ++it; } for (;it!=it1;++it){ G *=x; G += it->g; } g=G%modulo; } else { for (;it!=it1;++it){ g = (g*longlong(x)+it->g)%modulo; } } } g=smod(g,modulo); if (g!=0){ if (kt!=ktend){ *kt=T_unsigned(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } continue; } if (modulo>=46340){ for (;it!=itend;++it){ const hashgcd_U & u=it->u; if (u(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } } break; } if (ucur-u==var2) g = (g*longlong(x)+ it->g)%modulo; else g = (g*longlong(powmod(x,(ucur-u)/var2,modulo))+ it->g)%modulo; ucur=u; } // end for } // end if modulo>=46340 else { for (;it!=itend;++it){ const hashgcd_U & u=it->u; if (u(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } } break; } if (ucur-u==var2) g = (g*x+ it->g)%modulo; else g = (g*powmod(x,(ucur-u)/var2,modulo)+ it->g)%modulo; ucur=u; } // end for } } // end else x=1 if (it==itend){ if (g!=0){ g = smod(g*longlong(powmod(x,(ucur-uend)/var2,modulo)),modulo); if (g!=0){ if (kt!=ktend){ *kt=T_unsigned(g,uend); ++kt; } else px.push_back(T_unsigned(g,uend)); } } } } if (kt!=ktend) px.erase(kt,ktend); return true; } // v <- v*k % m void mulmod(vector & v,int k,int m){ if (k==1) return; vector::iterator it=v.begin(),itend=v.end(); for (;it!=itend;++it){ type_operator_times_reduce(*it,k,*it,m); // *it = ((*it)*k)%m; } } // v <- v*k % m static void mulmod(vector< vector > & v,int k,int m){ if (k==1) return; vector< vector >::iterator it=v.begin(),itend=v.end(); for (;it!=itend;++it){ mulmod(*it,k,m); } } inline void mulmod(int k,vector & v,int m){ mulmod(v,k,m); } // v <- v+w % m static void addmod(vector & v,const vector & w,int m){ vector::const_iterator jt=w.begin(),jtend=w.end(); vector::iterator it=v.begin(),itend=v.end(); int ws=int(jtend-jt),vs=int(v.size()); if (ws>vs){ if ((int)v.capacity() tmp(ws); copy(v.begin(),v.end(),tmp.begin()+ws-vs); swap(v,tmp); } else { v.insert(v.begin(),ws-vs,0); } it=v.begin(); itend=v.end(); } for (it=itend-ws;it!=itend;++jt,++it){ *it = (*it+*jt)%m; } // trim resulting polynomial for (it=v.begin();it!=itend;++it){ if (*it) break; } if (it!=v.begin()) v.erase(v.begin(),it); } // v <- v+w % m static void addmod(vector< vector > & v,const vector< vector > & w,int m){ vector< vector >::iterator it=v.begin(),itend=v.end(); vector< vector >::const_iterator jt=w.begin(),jtend=w.end(); int addv=int(jtend-jt)-int(itend-it); if (addv>0){ v.insert(v.begin(),addv,vector(0)); it=v.begin(); itend=v.end(); } for (it=itend-(jtend-jt);it!=itend;++jt,++it){ addmod(*it,*jt,m); } } // v <- v+w % m static void addmod(vecteur & v,const vecteur & w,int m){ vecteur::const_iterator jt=w.begin(),jtend=w.end(); vecteur::iterator it=v.begin(),itend=v.end(); int ws=int(jtend-jt),vs=int(v.size()); if (ws>vs){ if ((int)v.capacity() & v,const vector & w,int m){ vector::iterator it=v.begin(),itend=v.end(); vector::const_iterator jt=w.begin(),jtend=w.end(); int addv=int(jtend-jt)-int(itend-it); if (addv>0){ v.insert(v.begin(),addv,0); it=v.begin(); itend=v.end(); } for (it=itend-(jtend-jt);it!=itend;++jt,++it){ *it = (*it-*jt)%m; } for (it=v.begin();it!=itend;++it){ if (*it) break; } if (it!=v.begin()) v.erase(v.begin(),it); } // v <- k*(v-w) % m static void mulsubmod(int k,vector & v,const vector & w,int m){ vector::iterator it=v.begin(),itend=v.end(); vector::const_iterator jt=w.begin(),jtend=w.end(); int addv=int(jtend-jt)-int(itend-it); if (addv>0){ v.insert(v.begin(),addv,0); it=v.begin(); itend=v.end(); } itend -= (jtend-jt); if (m<=46340){ if (2*k>m) k -= m; for (;it!=itend;++it){ *it = ((*it)*k)%m; } for (itend=v.end();it!=itend;++jt,++it){ *it = ((*it-*jt)*k)%m; } } else { for (;it!=itend;++it){ type_operator_times_reduce(*it,k,*it,m); // *it = ((*it)*k)%m; } for (itend=v.end();it!=itend;++jt,++it){ *it -= *jt; type_operator_times_reduce(*it,k,*it,m); // *it = ((*it)*k)%m; } } for (it=v.begin();it!=itend;++it){ if (*it) break; } if (it!=v.begin()) v.erase(v.begin(),it); } // eval p at x with respect to all but the last variable template static bool horner(const std::vector< T_unsigned > & p,const std::vector & x,const std::vector & vars,std::vector & px,int modulo){ int s=int(x.size()); // int xback=x.back(); int vs=int(vars.size()); if (s+1!=vs || vs<2) return false; // setdimerr(); hashgcd_U var=vars[vs-2],var2=vars.back(); int shift_var=find_shift(var),shift_var2=find_shift(var2); typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); if (is_zero(x)){ int pos1=0,pos2=int(itend-it),pos; for (;pos2-pos1>1;){ pos=(pos1+pos2)/2; if ((it+pos)->uu>=var) ++it; if (it==itend) px.clear(); else { int pxdeg = it->u >> shift_var2; px=vector(pxdeg+1); for (;it!=itend;++it){ px[pxdeg - (it->u >> shift_var2)]=it->g; } } } else { int pxdeg=degree_xn(p,shift_var,shift_var2); px=vector(pxdeg+1); vector::const_iterator jtbeg=vars.begin(),jtend=vars.end(),jt; --jtend; vector::const_iterator ktbeg=x.begin(),kt; hashgcd_U u,oldu=0; int oldfact=0; int inverse=x.back()?invmod(x.back(),modulo):0; for (;it!=itend;){ // group next monomials with same powers in x1..xn-1 u=(it->u/var)*var; int tmpdeg=(it->u-u)/var2; vector tmp(tmpdeg+1); for (;;++it){ if (it==itend || it->uu-u)/var2]=it->g; } } } // end else is_zero(x) // trim px typename vector::iterator lt=px.begin(),ltend=px.end(); for (;lt!=ltend;++lt){ if (!is_zero(*lt)) break; } if (lt!=px.begin()) px.erase(px.begin(),lt); return true; } template static void convert_back(const vector & v,hashgcd_U var,vector< T_unsigned > & p){ p.clear(); typename vector::const_iterator it=v.begin(),itend=v.end(); unsigned s=unsigned(itend-it); p.reserve(s); hashgcd_U u=var*(s-1); for (;it!=itend;u-=var,++it){ if (!is_zero(*it)) p.push_back(T_unsigned(*it,u)); } } static void convert_back(const vector< vector > & v,hashgcd_U varxn,hashgcd_U var2,vector< T_unsigned > & p){ vector< vector >::const_iterator jt=v.begin(),jtend=v.end(); vector< T_unsigned >::iterator kt=p.begin(),ktend=p.end(); for (;jt!=jtend;++jt){ vector::const_iterator it=jt->begin(),itend=jt->end(); unsigned s=unsigned(itend-it); hashgcd_U u=var2*(s-1)+varxn*(unsigned(jtend-jt)-1); for (;it!=itend;u-=var2,++it){ if (*it!=0){ if (kt!=ktend){ *kt=T_unsigned(*it,u); ++kt; } else p.push_back(T_unsigned(*it,u)); } } } if (kt!=ktend) p.erase(kt,ktend); } static void convert(const vector< T_unsigned > & p,hashgcd_U var,vector & v,int modulo){ v.clear(); vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); if (it==itend) return; hashgcd_U u=it->u; unsigned s=u/var; v=vector(s+1); for (;it!=itend;++it){ v[s-it->u/var]=it->g<0?it->g+modulo:it->g; } } #if 0 static void convert(const vector< T_unsigned > & p,short int var,vector & v,int modulo){ vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); if (it==itend){ v.clear(); return; } hashgcd_U u=it->u; unsigned s=u >> var; if (v.size()==s+1) fill(v.begin(),v.end(),0); else v=vector(s+1); for (;it!=itend;++it){ v[s-(it->u >> var)]=it->g<0?it->g+modulo:it->g; } } #endif /* void convert(const vector< T_unsigned > & p,hashgcd_U var,vector & v){ v.clear(); vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); if (it==itend) return; hashgcd_U u=it->u,prevu=u; unsigned s=u/var; v.reserve(s+1); v.push_back(it->g); for (++it;it!=itend;++it){ u=it->u; prevu -= var; if (u==prevu) v.push_back(it->g); else { v.insert(v.end(),(prevu-u)/var,0); v.push_back(it->g); prevu=u; } } } */ // Find non zeros coeffs of p template static int find_nonzero(const vector & p,index_t & res){ res.clear(); typename vector::const_iterator it=p.begin(),itend=p.end(); if (it==itend) return 0; int nzeros=0; for (;it!=itend;++it){ bool test=is_zero(*it); res.push_back(test?0:1); if (test) ++nzeros; } return nzeros; } template static hashgcd_U lcoeff(const vector< T_unsigned > & p,hashgcd_U var,hashgcd_U var2,vector & lp){ lp.clear(); typename vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); if (it==itend) return 0; hashgcd_U u=it->u; int deg=(u%var)/var2; lp=vector(deg+1); u=(u/var)*var; for (;it!=itend;++it){ if (it->uu%var)/var2]=it->g; } return u; } static std::vector smod(const std::vector & v,int modulo){ std::vector res(v); std::vector::iterator it=res.begin(),itend=res.end(); for (;it!=itend;++it){ *it = smod(*it,modulo); } if (res.empty() || res.front()) return res; for (it=res.begin();it!=itend;++it){ if (*it) break; } return std::vector(it,itend); } static int hornermod(const vector & v,int alpha,int modulo,bool unsig=false){ vector::const_iterator it=v.begin(),itend=v.end(),it0; if (!alpha) return it==itend?0:v.back(); int res=0; if (alpha==1 && (itend-it)=-8 && alpha<=8){ it0=it+5*((itend-it)/5); for (;it!=it0;){ res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; res %= modulo; ++it; } for (;it!=itend;++it){ res *= alpha; res += *it; } return smod(res,modulo); } if (alpha>=-14 && alpha<=14){ it0=it+4*((itend-it)/4); for (;it!=it0;){ res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; res %= modulo; ++it; } for (;it!=itend;++it){ res *= alpha; res += *it; } return smod(res,modulo); } if (alpha>=-35 && alpha<=35){ it0=it+3*((itend-it)/3); for (;it!=it0;){ res *= alpha; res += *it; ++it; res *= alpha; res += *it; ++it; res *= alpha; res += *it; res %= modulo; ++it; } for (;it!=itend;++it){ res *= alpha; res += *it; } return smod(res,modulo); } if (alpha>=-214 && alpha<=214){ it0=it+2*((itend-it)/2); for (;it!=it0;){ res *= alpha; res += *it; ++it; res *= alpha; res += *it; res %= modulo; ++it; } for (;it!=itend;++it){ res *= alpha; res += *it; } return smod(res,modulo); } } // if (modulo<46430) else { longlong Res=res; if (alpha>=-8 && alpha<=8){ it0=it+5*((itend-it)/5); for (;it!=it0;){ Res *= alpha; Res += *it; ++it; Res *= alpha; Res += *it; ++it; Res *= alpha; Res += *it; ++it; Res *= alpha; Res += *it; ++it; Res *= alpha; Res += *it; Res %= modulo; ++it; } for (;it!=itend;++it){ Res *= alpha; Res += *it; } return smod(Res,modulo); } // end if alpha>=-84 and alpha<=84 } // end else (modulo>46430) #if defined _I386_ && !defined x86_64 if (unsig){ if (alpha<0) alpha += modulo; // it va dans ecx, itend dans ebx, modulo dans edi, alpha dans esi // eax::edx produit et division, eax contient res, asm volatile("movl $0x0,%%eax;\n\t" "cmpl %%ecx,%%ebx;\n\t" "je .Lend%=\n" ".Lloop%=:\t" "imul %%esi;\n\t" /* res<-res*alpha */ "addl (%%ecx),%%eax;\n\t" "adcl $0x0,%%edx; \n\t" /* res += *it */ "idivl %%edi; \n\t" /* res %= modulo */ "movl %%edx,%%eax; \n\t" "addl $4,%%ecx; \n\t" /* ++ *it */ "cmpl %%ecx,%%ebx;\n\t" /* it==itend */ "jne .Lloop%=\n" ".Lend%=:\t" :"=a"(res) :"c"(it),"b"(itend),"D"(modulo),"S"(alpha) :"%edx" ); if (res>modulo) return res-modulo; return res; } #endif if (modulo<46340){ if (unsig){ unsigned Alpha=alpha<0?alpha+modulo:alpha; unsigned res=0; for (;it!=itend;++it){ res = (res*Alpha+unsigned(*it))%modulo; } if (res>(unsigned)(modulo/2)) return int(res)-modulo; else return int(res); } for (;it!=itend;++it){ res = (res*alpha+*it)%modulo; } } else { for (;it!=itend;++it){ res = (res*longlong(alpha)+*it)%modulo; } } register int tmp=res+res; if (tmp>modulo) return res-modulo; if (tmp<=-modulo) return res+modulo; return res; } // distribute multiplication static void distmult(const vector< T_unsigned > & p,const vector & v,vector< T_unsigned > & pv,hashgcd_U var,int modulo){ if (&pv==&p){ vector< T_unsigned > tmp; distmult(p,v,tmp,var,modulo); swap(pv,tmp); return; } pv.clear(); vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); vector::const_iterator jtbeg=v.begin(),jtend=v.end(),jt; int vs=int(jtend-jtbeg),j; pv.reserve((itend-it)*vs); if (modulo>=46340){ for (;it!=itend;++it){ for (jt=jtbeg,j=1;jt!=jtend;++j,++jt){ if (*jt) pv.push_back(T_unsigned((it->g*longlong(*jt))%modulo,it->u+(vs-j)*var)); } } } else { for (;it!=itend;++it){ for (jt=jtbeg,j=1;jt!=jtend;++j,++jt){ if (*jt) pv.push_back(T_unsigned((it->g*(*jt))%modulo,it->u+(vs-j)*var)); } } } } // distribute multiplication static void distmult(const vector< T_unsigned,hashgcd_U> > & p,const vector & v,vector< T_unsigned,hashgcd_U> > & pv,hashgcd_U var,int modulo){ if (&pv==&p){ vector< T_unsigned,hashgcd_U> > tmp; distmult(p,v,tmp,var,modulo); swap(pv,tmp); return; } pv.clear(); vector< T_unsigned,hashgcd_U> >::const_iterator it=p.begin(),itend=p.end(); vector tmp; vector::const_iterator jtbeg=v.begin(),jtend=v.end(),jt; int vs=int(jtend-jtbeg),j; pv.reserve((itend-it)*vs); for (;it!=itend;++it){ for (jt=jtbeg,j=1;jt!=jtend;++j,++jt){ if (*jt){ tmp=it->g; mulmod(*jt,tmp,modulo); pv.push_back( T_unsigned,hashgcd_U>(tmp,it->u+(vs-j)*var)); } } } } #if 0 static void distmult(const vector & p,const vector & v,vector< vector > & pv,int modulo){ vector::const_iterator it=p.begin(),itend=p.end(),jt=v.begin(),jtend=v.end(); pv.clear(); pv.insert(pv.end(),(itend-it),vector(0)); if (modulo<=46340){ for (int i=0;it!=itend;++it,++i){ vector & w = pv[i]; const int & fact = *it; if (fact){ w.reserve(v.size()); for (jt=v.begin();jt!=jtend;++jt) w.push_back(fact*(*jt) % modulo); } } } else { for (int i=0;it!=itend;++it,++i){ vector & w = pv[i]; const int & fact = *it; if (fact){ w.reserve(v.size()); for (jt=v.begin();jt!=jtend;++jt) w.push_back( (longlong(fact)*(*jt)) % modulo); } } } } static void mulmod(const vector & v,int m,vector & w,int modulo){ if (&v==&w){ mulmod(w,m,modulo); return; } vector::const_iterator jt=v.begin(),jtend=v.end(); w.clear(); w.reserve(jtend-jt); if (modulo>=46340){ for (;jt!=jtend;++jt) w.push_back(smod(*jt * longlong(m),modulo)); } else { for (;jt!=jtend;++jt) w.push_back(smod(*jt * m,modulo)); } } #endif // if all indices are == return -2 // if all indices corr. to u1 are <= to u2 return 1, // if all are >= to u2 return 0 // otherwise return -1 int compare(hashgcd_U u1,hashgcd_U u2,const vector & vars){ if (u1==u2) return -2; std::vector::const_iterator it=vars.begin(),itend=vars.end(); hashgcd_U r1,r2; int res=-2; for (;it!=itend;++it){ r1=u1%(*it); r2=u2%(*it); if (r1==r2) continue; if (res==-2){ res=r1= to u2 return 0 // otherwise return -1 static int compare(const index_t & u1,const index_t & u2){ if (u1==u2) return -2; index_t::const_iterator it=u1.begin(),itend=u1.end(),jt=u2.begin(); int res=-2; for (;it!=itend;++it,++jt){ if (*it==*jt) continue; if (res==-2){ res=*it<*jt; continue; } if (*it<*jt){ if (res) continue; return -1; } else { if (!res) continue; return -1; } } return res; } inline bool is_one(const vector< T_unsigned > & p){ return p.size()==1 && p.front().g==1 && p.front().u==0; } // Note that smallmult may fail if the degree of a and b * modulo^2 // overflows in a longlong, so keep modulo not too large static void smallmult(const vector::const_iterator & ita0,const vector::const_iterator & ita_end,const vector::const_iterator & itb0,const vector::const_iterator & itb_end,vector & new_coord,int modulo){ longlong test=longlong(modulo)*std::min(ita_end-ita0,itb_end-itb0); bool large=test/RAND_MAX>RAND_MAX/modulo; new_coord.clear(); vector::const_iterator ita_begin=ita0,ita=ita0,itb=itb0; for ( ; ita!=ita_end; ++ita ){ vector::const_iterator ita_cur=ita,itb_cur=itb; if (large){ int res=0; for (;itb_cur!=itb_end;--ita_cur,++itb_cur) { res = (res + *ita_cur * longlong(*itb_cur))%modulo ; if (ita_cur==ita_begin) break; } new_coord.push_back(smod(res,modulo)); } else { longlong res=0; for (;itb_cur!=itb_end;--ita_cur,++itb_cur) { res += *ita_cur * longlong(*itb_cur) ; if (ita_cur==ita_begin) break; } new_coord.push_back(smod(res,modulo)); } } --ita; ++itb; for ( ; itb!=itb_end;++itb){ vector::const_iterator ita_cur=ita,itb_cur=itb; if (large){ int res=0; for (;;) { res = (res + *ita_cur * longlong(*itb_cur))%modulo ; if (ita_cur==ita_begin) break; --ita_cur; ++itb_cur; if (itb_cur==itb_end) break; } new_coord.push_back(smod(res,modulo)); } else { longlong res= 0; for (;;) { res += *ita_cur * longlong(*itb_cur) ; if (ita_cur==ita_begin) break; --ita_cur; ++itb_cur; if (itb_cur==itb_end) break; } new_coord.push_back(smod(res,modulo)); } } } #if 0 static void smallmult1(const vector< T_unsigned > & p,const vector< T_unsigned > & q,vector< T_unsigned > & res,hashgcd_U var,int modulo){ vector p1,q1,r1; convert(p,var,p1,modulo); convert(q,var,q1,modulo); r1.reserve(p1.size()+q1.size()-1); smallmult(p1.begin(),p1.end(),q1.begin(),q1.end(),r1,modulo); convert_back(r1,var,res); } #endif static void convert(const vector< T_unsigned > & p,hashgcd_U var,hashgcd_U var2,vector< vector > & v,int modulo){ if (p.empty()){ v.clear(); return; } int dim1=p.front().u/var,dim2; if (int(v.size())!=dim1+1){ vector< vector > tmp(dim1+1); swap(tmp,v); } vector< T_unsigned >::const_iterator it=p.begin(),itend=p.end(); hashgcd_U u,uend; int deg,prevdeg=0; for (;it!=itend;){ u=it->u; deg=u/var; if (prevdeg){ for (--prevdeg;prevdeg>deg;--prevdeg){ v[dim1-prevdeg].clear(); } } prevdeg=deg; uend=deg*var; dim2=(u-uend)/var2; vector & vi = v[dim1-deg]; if (int(vi.size())!=dim2+1){ vector tmp(dim2+1); swap(vi,tmp); } else fill(vi.begin(),vi.end(),0); if (dim2u==uend){ vector::iterator jt=vi.begin(),jtend=vi.end(); for (;jt!=jtend;++it,++jt){ *jt = it->g; if (*jt>=0) continue; *jt += modulo; } continue; } for (;it!=itend;++it){ u=it->u; if (ug<0?it->g+modulo:it->g; } } for (--prevdeg;prevdeg>=0;--prevdeg) v[dim1-prevdeg].clear(); } static void horner_back(const vector< vector > & v,int x,vector & vx,int modulo,int maxdeg=-1,bool unsig=false){ vector< vector >::const_iterator it=v.begin(),itend=v.end(); if (maxdeg>=0 && maxdeg::iterator jt=vx.begin(),jtend=vx.end(); if (jtend-jt>=itend-it){ for (;it!=itend;++it,++jt){ *jt=hornermod(*it,x,modulo,unsig); } if (jt!=jtend) vx.erase(jt,jtend); } else { vx.clear(); vx.reserve(itend-it); for (;it!=itend;++it){ vx.push_back(hornermod(*it,x,modulo,unsig)); } } } static void horner_front(const vector< vector > & v,int x,vector & vx,int modulo){ if (!x && !v.empty()){ vx=v.back(); return; } vx.clear(); vector< vector >::const_iterator it=v.begin(),itend=v.end(); for (;it!=itend;++it){ mulmod(vx,x,modulo); addmod(vx,*it,modulo); } // trim(vx); } static bool is_equal_modulo(const vector & v,const vector & w,int modulo){ vector::const_iterator it=v.begin(),itend=v.end(),jt=w.begin(),jtend=w.end(); if (itend-it!=jtend-jt) return false; for (;it!=itend;++jt,++it){ if ((*it-*jt)%modulo) return false; } return true; } static bool is_p_a_times_b(const vector & p,const vector & a,const vector & b,int modulo,int maxdeg){ int as=int(a.size()),bs=int(b.size()),ps=int(p.size()); if (ps!=as+bs-1 && ps!=maxdeg+1) return false; if (ps<=maxdeg){ vector r; smallmult(a.begin(),a.end(),b.begin(),b.end(),r,modulo); return is_equal_modulo(r,p,modulo); } longlong test=longlong(modulo)*std::min(maxdeg+1,std::min(as,bs)); bool large=test/RAND_MAX>RAND_MAX/modulo; int j; vector::const_iterator it,itbeg=a.begin(),jt,jtend=b.end(); for (int i=0;i<=maxdeg;++i){ // degree i==? p[ps-i-1]= sum_j a[as-1-j]*b[bs-1-i+j] // starting value for j: 0 or i+1-bs j=i+1-bs; if (j<0) j=0; it=itbeg+as-1-j; jt=b.begin()+bs-1-i+j; // end value: a.begin() or b.end() if (large){ int res = - p[ps-i-1]; for (;jt!=jtend;--it,++jt){ res = (res+(*it)* longlong(*jt))%modulo; if (it==itbeg) break; } if (res) return false; } else { longlong res = - p[ps-i-1]; for (;jt!=jtend;--it,++jt){ res += (*it)*longlong(*jt); if (it==itbeg) break; } if (res%modulo) return false; } } return true; } static bool is_p_a_times_b(const vector< T_unsigned > & p,const vector< T_unsigned > & a,const vector< T_unsigned > & b,const std::vector & vars,int modulo,int maxtotaldeg){ if (a.empty() || b.empty()) return p.empty(); if (p.empty()) return false; int dim=int(vars.size())-1; if (dim<=0) return false; // setdimerr(); // double as=a.size(),bs=b.size(); double ps=double(p.size()); hashgcd_U var2=vars[dim-1]; if (dim==1){ // ? dense vector multiplication or sparse mult ? unsigned adeg=a.front().u/var2,bdeg=b.front().u/var2,pdeg=p.front().u/var2; if (pdeg!=adeg+bdeg && int(pdeg)!=maxtotaldeg) return false; // double timesparse=as*bs*std::log(ps); // double timedense=(adeg+1)*(bdeg+1); if (true // timedense p1,a1,b1,r1; convert(a,var2,a1,modulo); convert(b,var2,b1,modulo); convert(p,var2,p1,modulo); return is_p_a_times_b(p1,a1,b1,modulo,maxtotaldeg); } // FIXME take maxdeg in account vector< T_unsigned > r; smallmult(a,b,r,modulo,size_t(ps+1)); smallsub(r,p,r,modulo); return r.empty(); } hashgcd_U var=vars[dim-2]; int shift_var=find_shift(var),shift_var2=find_shift(var2); unsigned adeg=degree_xn(a,shift_var,shift_var2),bdeg=degree_xn(b,shift_var,shift_var2),pdeg=degree_xn(p,shift_var,shift_var2); int ntests=std::min(int(pdeg),maxtotaldeg); if (pdeg!=adeg+bdeg && int(pdeg)!=maxtotaldeg) return false; // multiplication requires as*bs*ln(ps) operations // alternative is checking pdeg+1 values of x=xn // time for 1 check: ps+as+bs (horner)+ dim^2*as/adeg*bs/bdeg*ln(ps) // double timemult=3*as*bs*std::log(ps); // double timeeval=(ntests+1)*(ps+as+bs+dim*dim*as/adeg*bs/bdeg*std::log(ps)); if (false // timemult > r; smallmult(a,b,r,modulo,size_t(ps)); smallsub(r,p,r,modulo); // FIXME take maxdeg in account return r.empty(); } if (dim==2){ vector< vector > pv,av,bv; vector pi,ai,bi,ri; convert(p,var,var2,pv,modulo); convert(a,var,var2,av,modulo); convert(b,var,var2,bv,modulo); for (int i=0;i<=ntests;++i){ int alpha=i%2?-(i+1)/2:i/2; horner_back(pv,alpha,pi,modulo,maxtotaldeg,true); horner_back(av,alpha,ai,modulo,maxtotaldeg,true); horner_back(bv,alpha,bi,modulo,maxtotaldeg,true); if (!is_p_a_times_b(pi,ai,bi,modulo,maxtotaldeg)) return false; --maxtotaldeg; } return true; } vector< T_unsigned > pi,ai,bi; std::vector vars_truncated(vars); vars_truncated.pop_back(); for (int i=0;i<=ntests;++i){ if (!horner(p,i,vars_truncated,pi,modulo,maxtotaldeg) || !horner(a,i,vars_truncated,ai,modulo,maxtotaldeg) || !horner(b,i,vars_truncated,bi,modulo,maxtotaldeg)) return false; if (!is_p_a_times_b(pi,ai,bi,vars_truncated,modulo,maxtotaldeg)) return false; --maxtotaldeg; } return true; } /* template void divided_differences(const vector & x,vector< vector< T_unsigned > > & res,int modulo){ int s=x.size(); int fact; for (int k=1;k=k;--j){ smallsub(res[j],res[j-1],res[j],modulo); fact=invmod(x[j]-x[j-k],modulo); if (fact!=1) smallmult(fact,res[j],res[j],modulo); } } } */ template static void divided_differences(const vector & x,vector< vector< T_unsigned > > & res,int modulo){ int s=int(x.size()); int fact; vector< T_unsigned > tmp; for (int k=1;k=k;--j){ smallsub(res[j],res[j-1],tmp,modulo); swap(tmp,res[j]); fact=invmod(x[j]-x[j-k],modulo); if (fact!=1) smallmult(fact,res[j],res[j],modulo); } } } // Lagrange interpolation at x/y template static void interpolate(const vector & x,vector< vector< T_unsigned > > & diff,vector< T_unsigned > & res,hashgcd_U varx,int modulo){ divided_differences(x,diff,modulo); int s=int(diff.size()); vector interp(1,1); res=diff.front(); int alpha; vector< T_unsigned > tmp,tmp2; for (int j=1;j0;--i){ interp[i]=smod(-longlong(alpha)*interp[i-1]+interp[i],modulo); } distmult(diff[j],interp,tmp,varx,modulo); smalladd(res,tmp,tmp2,modulo); swap(tmp2,res); } /* res=diff[s-1]; vector< T_unsigned > res_shift,res_times; for (int i=s-2;i>=0;--i){ // res = res*(x-x[i])+diff[i]; smallshift(res,varx,res_shift); smallmult(-x[i],res,res_times,modulo); smalladd(res_times,diff[i],res_times,modulo); smalladd(res_shift,res_times,res,modulo); } */ } static void divided_differences(const vector & x,vector< vector > & res,int modulo){ int s=int(x.size()); int fact; for (int k=1;k=k;--j){ fact=invmod(x[j]-x[j-k],modulo); mulsubmod(fact,res[j],res[j-1],modulo); /* submod(res[j],res[j-1],modulo); if (fact!=1) mulmod(res[j],fact,modulo); */ } } } // Lagrange interpolation at x/y static void interpolate(const vector & x,vector< vector > & diff,vector< vector > & res,int modulo){ divided_differences(x,diff,modulo); // CERR << "end diff div " << CLOCK() << endl; res.clear(); int s=int(diff.size()),alpha; int ysize=0,cur; for (int i=0;iysize ) ysize=cur; } res.reserve(ysize); for (int i=ysize-1;i>=0;--i){ res.push_back(vector(0)); vector & curx = res.back(); vector & cury = diff[s-1]; if ( (cur=int(cury.size())) >i) curx.push_back(cury[cur-1-i]); for (int j=s-2;;--j){ // multiply curx by (x-x[j]) alpha=-x[j]; if (!curx.empty()){ curx.push_back( (longlong(alpha)*curx.back()) % modulo ); vector::iterator it=curx.end()-2,itbeg=curx.begin(); for (;it!=itbeg;){ int & curxk=*it; --it; type_operator_plus_times_reduce(alpha,*it,curxk,modulo); } } // add diff[j] vector & cury = diff[j]; if ( (cur=int(cury.size())) >i){ if (curx.empty()) curx.push_back(cury[cur-1-i]); else curx.back() = (curx.back() + cury[cur-1-i])%modulo; } if (!j) break; } } /* vector interp(1,1); distmult(diff[0],interp,res,modulo); vector< vector >::iterator it=res.begin(),itend=res.end(); for (;it!=itend;++it) it->reserve(s); vector< vector > tmp; for (int j=1;j0;--i){ interp[i]=smod(-longlong(alpha)*interp[i-1]+interp[i],modulo); } distmult(diff[j],interp,tmp,modulo); addmod(res,tmp,modulo); } */ // unsigne res s=int(res.size()); for (int j=0;j & cur=res[j]; vector::iterator it=cur.begin(),itend=cur.end(); for (;it!=itend;++it){ if (*it<0) *it += modulo; } } } static void interpolate(const vector & x,vector< vector > & y,vector< T_unsigned > & res,int varxn,int var2,int modulo){ vector< vector > tmp; interpolate(x,y,tmp,modulo); convert_back(tmp,varxn,var2,res); } template struct gcd_call_param { vector * Delta ; vector * lcoeffp ; vector * lcoeffq ; vector * alphav; vector< vector > * pv ; vector< vector > * qv ; vector< vector > * dim2gcdv ; vector< vector > * dim2pcofactorv ; vector< vector > * dim2qcofactorv ; vector< T_unsigned > * p ; vector< T_unsigned > * q ; vector< vector< T_unsigned > > * gcdv ; vector< vector< T_unsigned > > * pcofactorv ; vector< vector< T_unsigned > > * qcofactorv ; index_t * pdeg ; index_t * qdeg ; const vector * vars ; vector * vars_truncated ; index_t * shift_vars ; index_t * shift_vars_truncated ; bool compute_cof ; bool compute_qcofactor ; bool dim2; const vector * pminptr; int modulo ; int vpos ; int nthreads ; int ext_gcd_ok ; // used for gcd over algebraic extension of Q }; static bool mod_gcd(const vector< T_unsigned > & p_orig,const vector< T_unsigned > & q_orig,int modulo,vector< T_unsigned > & d, vector< T_unsigned > & pcofactor, vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_pcofactor,bool compute_qcofactor,bool & divtest,vector< vector > & pv,vector< vector > & qv,int nthreads); #ifndef NO_TEMPLATE_MULTGCD static void * do_recursive_gcd_call(void * ptr_){ #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted) return 0; gcd_call_param * ptr = (gcd_call_param *) ptr_; vector & Delta = *ptr->Delta; vector & lcoeffp = *ptr->lcoeffp; vector & lcoeffq = *ptr->lcoeffq; vector & alphav = * ptr->alphav; vector< vector > & pv = *ptr->pv; vector< vector > & qv = *ptr->qv; vector dim2palpha; vector dim2qalpha; vector< vector > & dim2gcdv = *ptr->dim2gcdv; vector< vector > & dim2pcofactorv = *ptr->dim2pcofactorv; vector< vector > & dim2qcofactorv = *ptr->dim2qcofactorv; vector< T_unsigned > & p = * ptr->p; vector< T_unsigned > & q = * ptr->q; vector< vector< T_unsigned > > & gcdv = * ptr->gcdv; vector< vector< T_unsigned > > & pcofactorv = * ptr->pcofactorv; vector< vector< T_unsigned > > & qcofactorv = * ptr->qcofactorv; index_t & pdeg = * ptr->pdeg; index_t & qdeg = * ptr->qdeg; vector< T_unsigned > palpha,qalpha; index_t pdegalpha,qdegalpha; int vpos=ptr->vpos; int alpha1 = alphav[vpos]; int modulo = ptr->modulo; const vector & vars = * ptr->vars; vector & vars_truncated = * ptr->vars_truncated; // index_t & shift_vars = *ptr->shift_vars; index_t & shift_vars_truncated = *ptr->shift_vars_truncated; bool compute_cof = ptr->compute_cof; bool compute_qcofactor = ptr->compute_qcofactor; bool dim2 = ptr->dim2; int nthreads = ptr->nthreads; // Eval p and q at xn=alpha if (dim2){ horner_back(pv,alpha1,dim2palpha,modulo,-1,true); if ( int(dim2palpha.size())-1 != pdeg.front()) return 0; // convert(dim2palpha,varxn,palpha); horner_back(qv,alpha1,dim2qalpha,modulo,-1,true); if ( int(dim2qalpha.size())-1 != qdeg.front()) return 0; // convert(dim2qalpha,varxn,qalpha); } else { if (!horner(p,alpha1,vars,palpha,modulo,-1)) return 0; degree(palpha,shift_vars_truncated,pdegalpha); pdegalpha.push_back(pdeg.back()); if (pdegalpha!=pdeg) return 0; if (!horner(q,alpha1,vars,qalpha,modulo,-1)) return 0; degree(qalpha,shift_vars_truncated,qdegalpha); qdegalpha.push_back(qdeg.back()); if (qdegalpha!=qdeg) return 0; } if (dim2){ gcdsmallmodpoly(dim2palpha,dim2qalpha,modulo,dim2gcdv[vpos],compute_cof?&dim2pcofactorv[vpos]:0,(compute_cof && compute_qcofactor)?&dim2qcofactorv[vpos]:0); mulmod(dim2gcdv[vpos],smod(hornermod(Delta,alpha1,modulo)*longlong(invmod(dim2gcdv[vpos].front(),modulo)),modulo),modulo); if (compute_cof){ mulmod(dim2pcofactorv[vpos],smod(hornermod(lcoeffp,alpha1,modulo)*longlong(invmod(dim2pcofactorv[vpos].front(),modulo)),modulo),modulo); if (compute_qcofactor){ mulmod(dim2qcofactorv[vpos],smod(hornermod(lcoeffq,alpha1,modulo)*longlong(invmod(dim2qcofactorv[vpos].front(),modulo)),modulo),modulo); } } } else { vector< T_unsigned > & g=gcdv[vpos]; vector< T_unsigned > & gp=pcofactorv[vpos]; vector< T_unsigned > & gq=qcofactorv[vpos]; bool tmptestdiv; if (&pv && &qv){ if (!mod_gcd(palpha,qalpha,modulo,g,gp,gq,vars_truncated,compute_cof,compute_qcofactor,tmptestdiv,pv,qv,nthreads)){ g.clear(); return 0; } } else { vector< vector > pv1,qv1; if (!mod_gcd(palpha,qalpha,modulo,g,gp,gq,vars_truncated,compute_cof,compute_qcofactor,tmptestdiv,pv1,qv1,nthreads)){ g.clear(); return 0; } } smallmult(smod(hornermod(Delta,alpha1,modulo)*longlong(invmod(g.front().g,modulo)),modulo),g,g,modulo); if (compute_cof){ // adjust gp lcoeff smallmult(smod(hornermod(lcoeffp,alpha1,modulo)*longlong(invmod(gp.front().g,modulo)),modulo),gp,gp,modulo); if (compute_qcofactor){ // adjust gq cofactor smallmult(smod(hornermod(lcoeffq,alpha1,modulo)*longlong(invmod(gq.front().g,modulo)),modulo),gq,gq,modulo); } } } return ptr; } #endif //NO_TEMPLATE_MULTGCD // Modular gcd in "internal form" static bool mod_gcd(const vector< T_unsigned > & p_orig,const vector< T_unsigned > & q_orig,int modulo,vector< T_unsigned > & d, vector< T_unsigned > & pcofactor, vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_pcofactor,bool compute_qcofactor,bool & divtest,vector< vector > & pv,vector< vector > & qv,int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return false; #else divtest=true; if (p_orig.empty() || is_one(q_orig) ){ d=q_orig; if (compute_pcofactor) pcofactor=p_orig; if (compute_qcofactor){ if (p_orig.empty()){ qcofactor.clear(); qcofactor.push_back(T_unsigned(1,0)); } else qcofactor=q_orig; } return true; } if (q_orig.empty() || is_one(p_orig) ){ d=p_orig; if (compute_qcofactor) qcofactor=q_orig; if (compute_pcofactor){ if (q_orig.empty()){ pcofactor.clear(); pcofactor.push_back(T_unsigned(1,0)); } else pcofactor=p_orig; } return true; } if (&p_orig==&d || &q_orig==&d){ vector< T_unsigned > res; bool b=mod_gcd(p_orig,q_orig,modulo,res,pcofactor,qcofactor,vars,compute_pcofactor,compute_qcofactor,divtest,qv,pv,nthreads); swap(res,d); return b; } d.clear(); // Check dim, if 1 -> int dim=int(vars.size()); if (dim==1){ hashgcd_U var=vars.front(); vector pv,qv,dv,pvcof,qvcof; convert(p_orig,var,pv,modulo); convert(q_orig,var,qv,modulo); gcdsmallmodpoly(pv,qv,modulo,dv,compute_pcofactor?&pvcof:0,compute_qcofactor?&qvcof:0); convert_back(dv,var,d); if (compute_pcofactor){ convert_back(pvcof,var,pcofactor); } if (compute_qcofactor){ convert_back(qvcof,var,qcofactor); } return true; } std::vector vars_truncated(vars); vars_truncated.pop_back(); hashgcd_U varxn=vars_truncated.back(),var2=vars.back(); index_t shift_vars,shift_vars_truncated; if (!find_shift(vars,shift_vars)) return false; // setsizeerr(); shift_vars_truncated=shift_vars; shift_vars_truncated.pop_back(); short int shiftxn=shift_vars_truncated.back(),shift2=shift_vars.back(); // Make p and q primitive as polynomials in x1,...,xn-1 // with coeff polynomial in xn vector< T_unsigned > p(p_orig),q(q_orig),pcont,qcont,dcont,tmp; vector pcontxn,qcontxn,dcontxn,pcofcontxn,qcofcontxn; if (debug_infolevel>20-dim) CERR << "gcdmod threads " << nthreads << " content begin " << "dim " << dim << " " << CLOCK() << endl; pp_mod(p,0,modulo,varxn,var2,pcontxn); pp_mod(q,0,modulo,varxn,var2,qcontxn); gcdsmallmodpoly(pcontxn,qcontxn,modulo,dcontxn,compute_pcofactor?&pcofcontxn:0,compute_qcofactor?&qcofcontxn:0); if (debug_infolevel>20-dim) CERR << "gcdmod content in " << "dim " << dim << " " << CLOCK() << endl; // Make p and q primitive as polynomial in xn with coeff in x1...xn-1 vector< vector > dv,dpv,dim2gcdv,dim2pcofactorv,dim2qcofactorv; if (dim==2){ convert(p,varxn,var2,pv,modulo); if (is_front_primitive(pv,modulo)) pcont.push_back(T_unsigned(1,0)); else { pp_mod(p,0,modulo,vars,pcont,nthreads); if (!is_one(pcont)) convert(p,varxn,var2,pv,modulo); } convert(q,varxn,var2,qv,modulo); if (is_front_primitive(qv,modulo)) qcont.push_back(T_unsigned(1,0)); else { pp_mod(q,0,modulo,vars,qcont,nthreads); if (!is_one(qcont)) convert(q,varxn,var2,qv,modulo); } } else { pp_mod(p,0,modulo,vars,pcont,nthreads); pp_mod(q,0,modulo,vars,qcont,nthreads); } mod_gcd(pcont,qcont,modulo,dcont,pcofactor,qcofactor,vars_truncated,compute_pcofactor,compute_qcofactor,nthreads); // don't use pv and qv here! // multiply pcofactor and qcofactor by the initial contents dep. on xn if (debug_infolevel>20-dim) CERR << "gcdmod content end " << "dim " << dim << " " << CLOCK() << endl; if (compute_pcofactor){ convert_back(pcofcontxn,vars.back(),tmp); smallmult(pcofactor,tmp,pcofactor,modulo,0); } if (compute_qcofactor){ convert_back(qcofcontxn,vars.back(),tmp); smallmult(qcofactor,tmp,qcofactor,modulo,0); } distmult(dcont,dcontxn,dcont,var2,modulo); // ready for gcd computation by interpolation with respect to xn // first find degree of gcd with respect to xn int pxndeg=degree_xn(p,shiftxn,shift2),qxndeg=degree_xn(q,shiftxn,shift2),gcddeg=0; vector pb(pxndeg+1),qb(qxndeg+1),db,b(dim-1),bnext(dim-1); index_t vzero; // coeff of vzero correspond to zero or non zero int nzero=1; // Number of zero coeffs for (int essai=0;essai<2;){ if (debug_infolevel>20-dim) CERR << "gcdmod degree? " << essai << " dim " << dim << " " << CLOCK() << endl; if (dim==2){ horner_front(pv,b.front(),pb,modulo); horner_front(qv,b.front(),qb,modulo); } else { if (!horner(p,b,vars,pb,modulo) || !horner(q,b,vars,qb,modulo)) return false; } for (;;){ for (int i=0;i lcoeffp,lcoeffq,lcoeffg,Delta; hashgcd_U lcoeffpu,lcoeffqu; if (debug_infolevel>20-dim) CERR << "gcdmod lcoeff begin " << "dim " << dim << " " << CLOCK() << endl; lcoeffpu=lcoeff(p,varxn,var2,lcoeffp); lcoeffqu=lcoeff(q,varxn,var2,lcoeffq); gcdsmallmodpoly(lcoeffp,lcoeffq,modulo,Delta); if (debug_infolevel>20-dim){ CERR << "lcoeff p, q, gcd" << lcoeffp << "," << lcoeffq << "," << Delta << endl; CERR << "gcdmod lcoeff end " << "dim " << dim << " " << CLOCK() << endl; } // estimate time for full lift or division try // size=p.size()+q.size() // sumdeg=pxndeg+qxndeg // %age=gcddeg/min(pxndeg,qxndeg) // %age^dim*(1-%age)^dim*size^2 estimates the time for division try // gcddeg*size estimates the time for lifting to gcddeg // sumdeg*size estimates the time for full lifting // if sumdeg<(gcddeg+%age^dim*(1-%age)^dim*size) do full lifting int Deltadeg = int(Delta.size())-1,liftdeg=(compute_qcofactor?qxndeg:pxndeg)+Deltadeg; int gcddeg_plus_delta=gcddeg+Deltadeg; int liftdeg0=giacmax(liftdeg-gcddeg,gcddeg_plus_delta); // once liftdeg0 is reached we can replace g/gp/gq computation // by a check that d*dp=dxn*lcoeff(d*dp)/Delta at alpha // and d*dq=dxn*lcoeff(d*dq)/lcoeff(qxn) at alpha int sumdeg = pxndeg+qxndeg; double percentage = double(gcddeg)/giacmin(pxndeg,qxndeg); int sumsize = int(p.size()+q.size()); // ? add a malus factor for division double gcdlift=gcddeg+std::pow(percentage,dim)*std::pow(1-percentage,dim)*sumsize; bool compute_cof = dim==2 || sumdeg20-dim) CERR << "gcdmod degree begin " << "dim " << dim << " " << CLOCK() << " compute_cof " << compute_cof << "(" << sumdeg/gcdlift << ")" << endl; int ptotaldeg=degree(p,shift_vars,pdeg); int qtotaldeg=degree(q,shift_vars,qdeg); if (debug_infolevel>20-dim){ CERR << "pdeg " << pdeg << " " << ptotaldeg << endl; CERR << "qdeg " << qdeg << " " << qtotaldeg << endl; } if (debug_infolevel>20-dim) CERR << "gcdmod degree end " << "dim " << dim << " " << CLOCK() << endl; int spdeg=0,sqdeg=0; for (int i=0;i20-dim) CERR << "gcdmod find alpha dim " << dim << " " << CLOCK() << endl; if (debug_infolevel>25-dim) CERR << " p " << p << " q " << q << endl; vector< T_unsigned > palpha,qalpha,dp,dq; // d, dp and dq are the current interpolated values of gcd and cofactors vector alphav,dim2palpha,dim2qalpha,tmpcont; vector< vector< T_unsigned > > gcdv,pcofactorv,qcofactorv; // for dim 2 bool dim2 = dim==2 && compute_cof; if (dim2){ dim2gcdv.reserve(liftdeg0+1); dim2pcofactorv.reserve(liftdeg0+1); dim2qcofactorv.reserve(liftdeg0+1); } else { gcdv.reserve(liftdeg0+1); pcofactorv.reserve(liftdeg0+1); qcofactorv.reserve(liftdeg0+1); } gcd_call_param gcd_par; gcd_par.Delta=Δ gcd_par.lcoeffp=&lcoeffp; gcd_par.lcoeffq=&lcoeffq; gcd_par.alphav=&alphav; gcd_par.pv=&pv; gcd_par.qv=&qv; gcd_par.dim2gcdv=&dim2gcdv; gcd_par.dim2pcofactorv=&dim2pcofactorv; gcd_par.dim2qcofactorv=&dim2qcofactorv; gcd_par.p=&p; gcd_par.q=&q; gcd_par.gcdv=&gcdv; gcd_par.pcofactorv=&pcofactorv; gcd_par.qcofactorv=&qcofactorv; gcd_par.pdeg=&pdeg; gcd_par.qdeg=&qdeg; gcd_par.vars=&vars; gcd_par.vars_truncated=&vars_truncated; gcd_par.shift_vars=&shift_vars; gcd_par.shift_vars_truncated=&shift_vars_truncated; gcd_par.compute_cof=compute_cof; gcd_par.compute_qcofactor=compute_qcofactor; gcd_par.dim2=dim2; gcd_par.modulo=modulo; if (dim>3 || (dim==3 && sumsize*ptotaldeg*4 > modgcd_cachesize )){ gcd_par.nthreads=nthreads; nthreads=1; } else gcd_par.nthreads=1; if (debug_infolevel>20 && nthreads>1) CERR << "nthreads " << nthreads << " dim " << dim << " " << sumsize << " " << sumsize*ptotaldeg << endl; if (nthreads>gcddeg_plus_delta) nthreads=gcddeg_plus_delta+1; if (nthreads>1){ int todo=compute_cof?(liftdeg0+1):(gcddeg_plus_delta+1); double nth=todo/double(nthreads); // if (nthreads>=4 && (todo%nthreads==0)) nth = (todo+1)/double(nthreads); // keep one proc for a bad prime nth=std::ceil(nth); nthreads=int(std::ceil(todo/nth)); if (debug_infolevel>20-dim) CERR << "Using " << nthreads << " threads " << nth << " " << todo/nth << endl; } for (alpha=-1;;){ // Possible improvement: if gcddeg is high, cofactors will stabilize // soon, then gcd could be obtained by division instead of interp // First check if we are ready to interpolate if (!compute_cof && e>gcddeg_plus_delta){ if (dim2) interpolate(alphav,dim2gcdv,d,varxn,var2,modulo); else interpolate(alphav,gcdv,d,var2,modulo); if (debug_infolevel>20-dim) CERR << "gcdmod pp1mod dim " << dim << " " << CLOCK() << " d " << d << endl; vector< T_unsigned > pquo,qquo,tmprem,pD(d); pp_mod(pD,0,modulo,varxn,var2,tmpcont); // This removes the polynomial in xn that we multiplied by // (it was necessary to know the lcoeff of the interpolated poly) if (debug_infolevel>20-dim) CERR << "gcdmod check dim " << dim << " " << CLOCK() << endl; // Now, gcd divides pD for gcddeg+1 values of x1 // degree(pD)<=degree(gcd) if (hashdivrem(p,pD,pquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ // If pD divides both P and Q, then the degree wrt variables // x1,...,xn-1 is the right one (because it is <= since pD // divides the gcd and >= since pD(xn=one of the try) was a gcd // The degree in xn is the right one because of the condition // on the lcoeff // Note that the division test might be much longer than the // interpolation itself (e.g. if the degree of the gcd is small) // but it seems unavoidable, for example if // P=Y-X+X(X-1)(X-2)(X-3) // Q=Y-X+X(X-1)(X-2)(X-4) // then gcd(P,Q)=1, but if we take Y=0, Y=1 or Y=2 // we get gcddeg=1 (probably degree 1 for the gcd) // interpolation at X=0 and X=1 will lead to Y-X as candidate gcd // and even adding X=2 will not change it // We might remove division if we compute the cofactors of P and Q // if P=pD*cofactor is true for degree(P) values of x1 // and same for Q, and the degrees wrt xn of pD and cofactors // have sum equal to degree of P or Q + lcoeff then pD is the gcd if (hashdivrem(q,pD,qquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ smallmult(pD,dcont,d,modulo,0); smallmult(invmod(d.front().g,modulo),d,d,modulo); if (compute_pcofactor){ smallmult(pcofactor,pquo,pcofactor,modulo,0); smallmult(smod(longlong(p_orig.front().g)*invmod(pcofactor.front().g,modulo),modulo),pcofactor,pcofactor,modulo); } if (compute_qcofactor){ smallmult(qcofactor,qquo,qcofactor,modulo,0); smallmult(smod(longlong(q_orig.front().g)*invmod(qcofactor.front().g,modulo),modulo),qcofactor,qcofactor,modulo); } if (debug_infolevel>20-dim) CERR << "gcdmod found dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); return true; } // end if hashdivrem(q,...) } // end if hashdivrem(p,...) if (debug_infolevel>20-dim) CERR << "Gcdmod bad guess " << endl; // restart gcdv.clear(); alphav.clear(); dim2gcdv.clear(); pcofactorv.clear(); qcofactorv.clear(); dim2pcofactorv.clear(); dim2qcofactorv.clear(); e=0; } // end if (e>gcddeg+delta) if (compute_cof && e>liftdeg0 ){ // interpolate d and dp if (dim2){ interpolate(alphav,dim2gcdv,dv,modulo); if (debug_infolevel>20-dim) CERR << "end interpolate gcd " << CLOCK() << endl; convert_back(dv,varxn,var2,d); interpolate(alphav,dim2pcofactorv,dpv,modulo); convert_back(dpv,varxn,var2,dp); if (debug_infolevel>20-dim) CERR << "end interpolate p cof " << CLOCK() << endl; } else { interpolate(alphav,gcdv,d,var2,modulo); if (debug_infolevel>20-dim) CERR << "end interpolate gcd " << CLOCK() << endl; interpolate(alphav,pcofactorv,dp,var2,modulo); if (debug_infolevel>20-dim) CERR << "end interpolate p cof " << CLOCK() << endl; } // check that d(alpha)*dp(alpha)=palpha with lcoeff adjusted // for e<=liftdeg if (dim2){ vector dv1,dpv1; for (++alpha;e<=liftdeg;++e,++alpha){ alpha1=alpha%2?-(alpha+1)/2:alpha/2; while (hornermod(lcoeffp,alpha1,modulo)==0){ ++alpha; alpha1=alpha%2?-(alpha+1)/2:alpha/2; } #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || alpha>=modulo) return false; int maxtotaldeg=ptotaldeg+1-e; horner_back(pv,alpha1,dim2palpha,modulo,maxtotaldeg,true); horner_back(dv,alpha1,dv1,modulo,maxtotaldeg,true); horner_back(dpv,alpha1,dpv1,modulo,maxtotaldeg,true); if (debug_infolevel>20){ CERR << "palpha " << dim2palpha << endl; CERR << "gcd alpha " << dv1 << endl; CERR << "p-cof alpha " << dpv1 << endl; } mulmod(dpv1,smod(longlong(hornermod(pv.front(),alpha1,modulo))*invmod((hornermod(dv.front(),alpha1,modulo)*longlong(hornermod(dpv.front(),alpha1,modulo)))%modulo,modulo),modulo),modulo); if (!is_p_a_times_b(dim2palpha,dpv1,dv1,modulo,maxtotaldeg)){ e=liftdeg0+1; break; } } } else { vector< T_unsigned > g,gp; for (++alpha;e<=liftdeg;++e,++alpha){ alpha1=alpha%2?-(alpha+1)/2:alpha/2; while (hornermod(lcoeffp,alpha1,modulo)==0){ ++alpha; alpha1=alpha%2?-(alpha+1)/2:alpha/2; } #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || alpha>=modulo) return false; int maxtotaldeg=ptotaldeg+1-e; if (!horner(p,alpha1,vars,palpha,modulo,maxtotaldeg)) return false; if (debug_infolevel>20-dim) CERR << "gcdmod horner d " << alpha << " dim " << dim << " " << CLOCK() << endl; if (!horner(d,alpha1,vars,g,modulo,maxtotaldeg)) return false; if (debug_infolevel>20-dim) CERR << "gcdmod horner dp " << alpha << " dim " << dim << " " << CLOCK() << endl; if (!horner(dp,alpha1,vars,gp,modulo,maxtotaldeg)) return false; smallmult(smod(longlong(palpha.back().g)*invmod((gp.back().g*longlong(g.back().g))%modulo,modulo),modulo),gp,gp,modulo); if (!is_p_a_times_b(palpha,gp,g,vars,modulo,maxtotaldeg)){ // Bad guess, go find some new gcd and interpolate e=liftdeg0+1; break; } } // end for ( e loop ) } if (e>liftdeg){ // enough evaluation point // divide d,dp,dq by their content in xn pp_mod(d,0,modulo,varxn,var2,tmpcont); pp_mod(dp,0,modulo,varxn,var2,tmpcont); // check xn degrees of d+dp=degree(pxn), d+dq=degree(qxn) int dxndeg=degree_xn(d,shiftxn,shift2),dpxndeg=degree_xn(dp,shiftxn,shift2); // int dqxndeg=degree_xn(dq,shiftxn,shift2); if ( dxndeg+dpxndeg==pdeg.back() ){ smallmult(d,dcont,d,modulo,0); if (compute_pcofactor){ smallmult(dp,pcofactor,pcofactor,modulo,0); smallmult(smod(longlong(p_orig.front().g)*invmod(pcofactor.front().g,modulo),modulo),pcofactor,pcofactor,modulo); } if (compute_qcofactor){ if (dim2) interpolate(alphav,dim2qcofactorv,dq,varxn,var2,modulo); else interpolate(alphav,qcofactorv,dq,var2,modulo); pp_mod(dq,0,modulo,varxn,var2,tmpcont); smallmult(dq,qcofactor,qcofactor,modulo,0); smallmult(smod(longlong(q_orig.front().g)*invmod(qcofactor.front().g,modulo),modulo),qcofactor,qcofactor,modulo); } if (debug_infolevel>20-dim) CERR << "gcdmod end dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); divtest=false; return true; } // failure, restart gcdv.clear(); alphav.clear(); dim2gcdv.clear(); pcofactorv.clear(); qcofactorv.clear(); dim2pcofactorv.clear(); dim2qcofactorv.clear(); e=0; } // end if (e>lifdeg) } // end if (compute_cof && e in liftdeg0..liftdeg) // *************************************************** // // Not ready to interpolate, try to eval new points // *************************************************** // for (int thread=0;thread20-dim) CERR << "gcdmod eval alpha1=" << alpha1 << " dim " << dim << " " << CLOCK() << endl; break; } // end for (;;) // alpha is probably admissible // (we will test later if degree of palpha/qalpha wrt x1...xn-1 is max) // prepare room for gcd and cofactors if (debug_infolevel>25-dim) CERR << "dim " << dim << " palpha " << palpha << " qalpha " << qalpha << endl ; alphav.push_back(alpha1); if (dim2){ dim2gcdv.push_back(vector(0)); dim2pcofactorv.push_back(vector(0)); dim2qcofactorv.push_back(vector(0)); } else { gcdv.push_back(vector< T_unsigned >(0)); pcofactorv.push_back(vector< T_unsigned >(0)); qcofactorv.push_back(vector< T_unsigned >(0)); } } // end for (int thread=0;thread > gcd_call_param_v(nthreads,gcd_par); #ifdef HAVE_PTHREAD_H pthread_t tab[nthreads-1]; #endif for (int thread=0;thread2 && gcddeg-nzero==e){ // We have enough evaluations, let's try SPMOD // Build the matrix, each line has coeffs / vzero vector< vector > m,minverse; for (int j=0;j<=e;++j){ index_t::reverse_iterator it=vzero.rbegin(),itend=vzero.rend(); vector line; for (int p=alphav[j],pp=1;it!=itend;++it,pp=smod(p*pp,modulo)){ if (*it) line.push_back(pp); } reverse(line.begin(),line.end()); m.push_back(line); } // assume gcd is the vector of non zero coeffs of the gcd in x^n // we have the relation // m*gcd=gcdv // invert m (if invertible) longlong det_mod_p; if (smallmodinv(m,minverse,modulo,det_mod_p) && det_mod_p){ // hence gcd=minverse*gcdv, where the i-th component of gcd // must be "multiplied" by xn^degree_corresponding_vzero[i] vector< vector< T_unsigned > > minversegcd(e+1); for (int j=0;j<=e;++j){ for (int k=0;k<=e;++k){ vector< T_unsigned > tmp(gcdv[k]); smallmult(minverse[j][k],tmp,tmp,modulo); smalladd(minversegcd[j],tmp,minversegcd[j],modulo); // CERR << minversegcd[j] << endl; } } vector< T_unsigned > trygcd,pquo,qquo,tmprem; index_t::const_iterator it=vzero.begin(),itend=vzero.end(); int deg=int(itend-it)-1; for (int j=0;it!=itend;++it,--deg){ if (!*it) continue; smallshift(minversegcd[j],deg*var2,minversegcd[j]); smalladd(trygcd,minversegcd[j],trygcd,modulo); ++j; } // Check if trygcd is the gcd! pp_mod(trygcd,0,modulo,varxn,var2,tmpcont); if (hashdivrem(p,trygcd,pquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ if (hashdivrem(q,trygcd,qquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ smallmult(trygcd,dcont,d,modulo,0); smallmult(invmod(d.front().g,modulo),d,d,modulo); if (compute_pcofactor){ smallmult(pcofactor,pquo,pcofactor,modulo,0); smallmult(smod(longlong(p_orig.front().g)*invmod(pcofactor.front().g,modulo),modulo),pcofactor,pcofactor); } if (compute_qcofactor){ smallmult(qcofactor,qquo,qcofactor,modulo,0); smallmult(smod(longlong(q_orig.front().g)*invmod(qcofactor.front().g,modulo),modulo),qcofactor,qcofactor); } if (debug_infolevel>20-dim) CERR << "gcdmod found dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); return true; } // end q divisible by trygcd } // end p divisible by trygcd } // end m invertible } // end if (dim>2 && ) } if (debug_infolevel>20-dim) CERR << "gcdmod interp dim " << dim << " " << CLOCK() << endl; ++e; continue; } // end gdeg==delta if (comp==0 || comp==-1){ if (debug_infolevel>20-dim) CERR << "Bad reduction " << alphav[vpos] << endl; // bad reduction: all indices of gdeg are >= to delta and gdeg!=delta alphav.erase(alphav.begin()+vpos); if (dim2){ dim2gcdv.erase(dim2gcdv.begin()+vpos); dim2pcofactorv.erase(dim2pcofactorv.begin()+vpos); dim2qcofactorv.erase(dim2qcofactorv.begin()+vpos); } else { gcdv.erase(gcdv.begin()+vpos); pcofactorv.erase(pcofactorv.begin()+vpos); qcofactorv.erase(qcofactorv.begin()+vpos); } if (comp==0) continue; } // previous alpha where bad reduction if (debug_infolevel>20-dim && vpos) CERR << "Bads reductions " << alphav[vpos-1] << endl; alphav.erase(alphav.begin(),alphav.begin()+vpos); if (dim2){ dim2gcdv.erase(dim2gcdv.begin(),dim2gcdv.begin()+vpos); dim2pcofactorv.erase(dim2pcofactorv.begin(),dim2pcofactorv.begin()+vpos); dim2qcofactorv.erase(dim2qcofactorv.begin(),dim2qcofactorv.begin()+vpos); } else { gcdv.erase(gcdv.begin(),gcdv.begin()+vpos); pcofactorv.erase(pcofactorv.begin(),pcofactorv.begin()+vpos); qcofactorv.erase(qcofactorv.begin(),qcofactorv.begin()+vpos); } if (comp==-1){ e=0; continue; } // restart everything with this value of alpha // this will happen (almost all the time) at first iteration delta=gdeg; e=1; continue; } // end for (int thread=0;thread > & p_orig,const std::vector< T_unsigned > & q_orig,int modulo,std::vector< T_unsigned > & d, std::vector< T_unsigned > & pcofactor, std::vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_pcofactor,bool compute_qcofactor,int nthreads){ bool divtest; vector< vector > pv,qv; return mod_gcd(p_orig,q_orig,modulo,d,pcofactor,qcofactor,vars,compute_pcofactor,compute_qcofactor,divtest,pv,qv,nthreads); } bool mod_gcd(const std::vector< T_unsigned > & p_orig,const std::vector< T_unsigned > & q_orig,int modulo,std::vector< T_unsigned > & d, std::vector< T_unsigned > & pcofactor, std::vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_cofactors,int nthreads){ return mod_gcd(p_orig,q_orig,modulo,d,pcofactor,qcofactor,vars,compute_cofactors,compute_cofactors,nthreads); } #if 0 static void smod(vector< T_unsigned > & p_orig,int modulo){ vector< T_unsigned >::iterator it=p_orig.begin(),itend=p_orig.end(); for (;it!=itend;++it){ it->g=smod(it->g,modulo); } } #endif static void smod(const vector< T_unsigned > & p_orig,int modulo,vector< T_unsigned > & p){ p.clear(); vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); p.reserve(itend-it); int x; for (;it!=itend;++it){ x=smod(it->g,modulo).val; if (x) p.push_back(T_unsigned(x,it->u)); } } static void unmod(const vector< T_unsigned > & p_orig,vector< T_unsigned > & p,int modulo){ p.clear(); vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); for (;it!=itend;++it){ p.push_back(T_unsigned(smod(it->g,modulo),it->u)); } } static void complex_unmod(const vector< T_unsigned > & p1,const vector< T_unsigned > & p2,vector< T_unsigned > & p,int i,int modulo){ p.clear(); vector< T_unsigned >::const_iterator it=p1.begin(),itend=p1.end(); vector< T_unsigned >::const_iterator jt=p2.begin(),jtend=p2.end(); // real part is halfsum of p1,p2, imaginary part halfdiff of p1-p2 divided by i (mod modulo) longlong inv2=invmod(2,modulo), inv2i=invmod(2*i,modulo); for (;it!=itend && jt!=jtend;){ if (it->u==jt->u){ p.push_back(T_unsigned(gen(smod(inv2*(it->g+jt->g),modulo),smod(inv2i*(it->g-jt->g),modulo)),it->u)); ++it; ++jt; continue; } if (it->u>jt->u){ p.push_back(T_unsigned(gen(smod(inv2*(it->g),modulo),smod(inv2i*(it->g),modulo)),it->u)); ++it; } else { p.push_back(T_unsigned(gen(smod(inv2*(jt->g),modulo),smod(inv2i*(-jt->g),modulo)),it->u)); ++jt; } } for (;it!=itend;++it){ p.push_back(T_unsigned(gen(smod(inv2*(it->g),modulo),smod(inv2i*(it->g),modulo)),it->u)); } for (;jt!=jtend;++jt){ p.push_back(T_unsigned(gen(smod(inv2*(jt->g),modulo),smod(inv2i*(-jt->g),modulo)),it->u)); } } static void complex_unmod_ext(const vector< T_unsigned > & p1,const vector< T_unsigned > & p2,vector< T_unsigned > & p,int i,int modulo){ p.clear(); vector< T_unsigned >::const_iterator it=p1.begin(),itend=p1.end(); vector< T_unsigned >::const_iterator jt=p2.begin(),jtend=p2.end(); // real part is halfsum of p1,p2, imaginary part halfdiff of p1-p2 divided by i (mod modulo) gen inv2=invmod(2,modulo), inv2i=invmod(2*i,modulo); for (;it!=itend && jt!=jtend;){ if (it->u==jt->u){ p.push_back(T_unsigned( smod(inv2*(it->g+jt->g),modulo)+cst_i*smod(inv2i*(it->g-jt->g),modulo) ,it->u)); ++it; ++jt; continue; } if (it->u>jt->u){ p.push_back(T_unsigned( smod(inv2*(it->g),modulo)+cst_i*smod(inv2i*(it->g),modulo) ,it->u)); ++it; } else { p.push_back(T_unsigned( smod(inv2*(jt->g),modulo)+cst_i*smod(inv2i*(-jt->g),modulo),it->u)); ++jt; } } for (;it!=itend;++it){ p.push_back(T_unsigned(smod(inv2*(it->g),modulo)+cst_i*smod(inv2i*(it->g),modulo),it->u)); } for (;jt!=jtend;++jt){ p.push_back(T_unsigned(smod(inv2*(jt->g),modulo)+cst_i*smod(inv2i*(-jt->g),modulo),it->u)); } } static void ichinrem(const vector< T_unsigned > & p_orig,const gen & modulo,vector< T_unsigned > & p,const gen & pimod){ vector< T_unsigned > q; vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); vector< T_unsigned >::const_iterator jt=p.begin(),jtend=p.end(); q.reserve(jtend-jt); for (;it!=itend && jt!=jtend;){ if (it->u==jt->u){ q.push_back(T_unsigned(ichinrem(it->g,jt->g,modulo,pimod),it->u)); ++it; ++jt; continue; } if (it->uu){ q.push_back(T_unsigned(ichinrem(0,jt->g,modulo,pimod),jt->u)); ++jt; } else { q.push_back(T_unsigned(ichinrem(it->g,zero,modulo,pimod),it->u)); ++it; } } for (;it!=itend;++it) q.push_back(T_unsigned(ichinrem(it->g,zero,modulo,pimod),it->u)); for (;jt!=jtend;++jt) q.push_back(T_unsigned(ichinrem(0,jt->g,modulo,pimod),jt->u)); swap(p,q); } static void ichinrem(const vector< T_unsigned > & p_orig,const gen & modulo,vector< T_unsigned > & p,const gen & pimod){ vector< T_unsigned > q; vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); vector< T_unsigned >::const_iterator jt=p.begin(),jtend=p.end(); q.reserve(jtend-jt); for (;it!=itend && jt!=jtend;){ if (it->u==jt->u){ q.push_back(T_unsigned(ichinrem(it->g,jt->g,modulo,pimod),it->u)); ++it; ++jt; continue; } if (it->uu){ q.push_back(T_unsigned(ichinrem(0,jt->g,modulo,pimod),jt->u)); ++jt; } else { q.push_back(T_unsigned(ichinrem(it->g,zero,modulo,pimod),it->u)); ++it; } } for (;it!=itend;++it) q.push_back(T_unsigned(ichinrem(it->g,zero,modulo,pimod),it->u)); for (;jt!=jtend;++jt) q.push_back(T_unsigned(ichinrem(0,jt->g,modulo,pimod),jt->u)); swap(p,q); } static void complex_ichinrem(const vector< T_unsigned > & p1,const vector< T_unsigned > & p2,int i,int modulo,vector< T_unsigned > & p,const gen & pimod){ vector< T_unsigned > p_orig,q; complex_unmod(p1,p2,p_orig,i,modulo); ichinrem(p_orig,modulo,p,pimod); } static void complex_ichinrem_ext(const vector< T_unsigned > & p1,const vector< T_unsigned > & p2,int i,int modulo,vector< T_unsigned > & p,const gen & pimod){ vector< T_unsigned > p_orig; complex_unmod_ext(p1,p2,p_orig,i,modulo); ichinrem(p_orig,modulo,p,pimod); } static gen max(const vector< T_unsigned > & p,GIAC_CONTEXT){ vector< T_unsigned >::const_iterator jt=p.begin(),jtend=p.end(); gen g,res; for (;jt!=jtend;++jt){ g=abs(jt->g,contextptr); if (is_strictly_greater(g,res,contextptr)) res=g; } return res; } // if res==0 set res to gcd of coeffs of p // if divide, divide p by res static gen ppz(vector< T_unsigned > & p,const gen & pgcd,bool divide){ vector< T_unsigned >::iterator jt=p.begin(),jtend=p.end(); gen res=pgcd; if (res==0){ for (;jt!=jtend;++jt){ if (jt->g.type==_VECT){ const_iterateur it=jt->g._VECTptr->begin(),itend=jt->g._VECTptr->end(); for (;it!=itend;++it){ res=gcd(res,*it,context0); } } else res=gcd(res,jt->g,context0); if (is_one(res)) return res; } } if (!divide) return res; for (jt=p.begin();jt!=jtend;++jt){ jt->g = jt->g/res; } return res; } static gen lcmdeno(vector< T_unsigned > & p){ vector< T_unsigned >::iterator jt=p.begin(),jtend=p.end(); gen res=1; for (;jt!=jtend;++jt){ if (jt->g.type==_FRAC){ res=lcm(res,jt->g._FRACptr->den); continue; } if (jt->g.type==_VECT){ const_iterateur it=jt->g._VECTptr->begin(),itend=jt->g._VECTptr->end(); for (;it!=itend;++it){ if (it->type==_FRAC) res=lcm(res,it->_FRACptr->den); if (is_undef(res)) return res; if (it->type==_POLY){ vector< monomial >::const_iterator kt=it->_POLYptr->coord.begin(),ktend=it->_POLYptr->coord.end(); for (;kt!=ktend;++kt){ if (kt->value.type==_FRAC) res=lcm(res,kt->value._FRACptr->den); } } } continue; } if (jt->g.type==_POLY){ vector< monomial >::const_iterator it=jt->g._POLYptr->coord.begin(),itend=jt->g._POLYptr->coord.end(); for (;it!=itend;++it){ if (it->value.type==_FRAC) res=lcm(res,it->value._FRACptr->den); } continue; } } if (is_one(res)) return res; for (jt=p.begin();jt!=jtend;++jt){ jt->g = jt->g * res; } return res; } // 1: integer, 2: gaussian integer, 3: ext, 4: ext with gaussint coeff, 0: other static int is_integer(const vector< T_unsigned > & p,gen & coefft){ vector< T_unsigned >::const_iterator jt=p.begin(),jtend=p.end(); int t=1; for (;jt!=jtend;++jt){ if (jt->g.is_integer()) continue; if (jt->g.type==_EXT){ if (t<3) t=t==1?3:4; if (coefft.type==_EXT){ if (*(coefft._EXTptr+1)!=*(jt->g._EXTptr+1)) return 0; } else { coefft=jt->g; } if (t==3 && !is_zero(im(*jt->g._EXTptr,context0))) t=4; continue; } if (jt->g.type==_POLY){ vector< monomial >::const_iterator it=jt->g._POLYptr->coord.begin(),itend=jt->g._POLYptr->coord.end(); for (;it!=itend;++it){ if (it->value.is_integer()) continue; if (!it->value.is_cinteger()) return 0; if (t!=3) t=2; else t=4; } continue; } if (!jt->g.is_cinteger()) return 0; if (t!=3 && t!=4) t=2; else t=4; } return t; } // reduce g mod modulo using i as square root of -1 static int complex_smod(const gen & g,int i,int modulo){ gen tmp=smod(re(g,context0)+i*im(g,context0),modulo); #ifndef NO_STDEXCEPT if (tmp.type!=_INT_) setsizeerr(gettext("complex_smod")); #endif return tmp.val; } static gen complex_smod_ext(const gen & g,int i,int modulo){ return smod(re(g,context0)+i*im(g,context0),modulo); } static void complex_smod(const vector< T_unsigned > & p_orig,int i,int modulo,vector< T_unsigned > & p){ vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); p.clear(); p.reserve(itend-it); for (;it!=itend;++it){ p.push_back(T_unsigned(complex_smod(it->g,i,modulo),it->u)); } } static bool complex_smod_ext(const vector< T_unsigned > & p_orig,int i,int modulo,vector< T_unsigned > & p){ vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); p.clear(); p.reserve(itend-it); vecteur x; gen tmp,itg; const_iterateur vt,vtend; for (;it!=itend;++it){ x.clear(); if (it->g.type!=_EXT){ tmp=complex_smod_ext(it->g,i,modulo); if (!is_zero(tmp)) x.push_back(tmp); } else { tmp=*it->g._EXTptr; if (tmp.type!=_VECT || tmp._VECTptr->empty()) return false; vt=tmp._VECTptr->begin(),vtend=tmp._VECTptr->end(); x.reserve(vtend-vt); for (;vt!=vtend;++vt){ itg=complex_smod_ext(*vt,i,modulo); if (is_zero(itg) && x.empty()) continue; x.push_back(itg.val); } } if (!x.empty()) p.push_back( T_unsigned (x,it->u) ); } return true; } static int modsqrtminus1(int modulo){ int i; for (int j=2;j > & p_orig,const vector< T_unsigned > & q_orig,vector< T_unsigned > & d, vector< T_unsigned > & pcofactor, vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_cofactors,int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return false; #else index_t shift_vars; gen coefft; int tp=is_integer(p_orig,coefft),tq=is_integer(q_orig,coefft); if ( // tp!=1 || tq!=1 tp==0 || tq==0 || !find_shift(vars,shift_vars) ) return false; bool is_complex=tp==2 || tq==2; compute_cofactors=true; // FIXME index_t pdeg,qdeg,pdegmod,qdegmod,gdegmod,gdegmod2,gdeg; degree(p_orig,shift_vars,pdeg); degree(q_orig,shift_vars,qdeg); gdeg=min(pdeg,qdeg); // gen m=30000,pimod=1; gen m=536871000,pimod=1; gen lcoeffp=p_orig.front().g,lcoeffq=q_orig.front().g,gcdlcoeff=gcd(lcoeffp,lcoeffq,context0); d.clear(); pcofactor.clear(); qcofactor.clear(); bool divtest; vector< vector > pv,qv; if (is_complex){ vector< T_unsigned > p1,q1,g1,pcof1,qcof1,p2,q2,g2,pcof2,qcof2; for (;;){ m=nextprime(m+1); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || m.type!=_INT_) return false; int modulo=m.val; // computing gcd in Z[i]: use a modulo =1[4], find gcd for both roots of -1 mod modulo if (modulo%4==3) continue; // find square root of -1 int i=modsqrtminus1(modulo); int lg1=complex_smod(gcdlcoeff,i,modulo),lg2=complex_smod(gcdlcoeff,-i,modulo); int lp1=complex_smod(lcoeffp,i,modulo),lp2=complex_smod(lcoeffp,-i,modulo); int lq1=complex_smod(lcoeffq,i,modulo),lq2=complex_smod(lcoeffq,-i,modulo); if (!lg1 || !lp1 || !lq1 || !lg2 || !lp2 || !lq2) continue; complex_smod(p_orig,i,modulo,p1); degree(p1,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; complex_smod(q_orig,i,modulo,q1); degree(q1,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; complex_smod(p_orig,-i,modulo,p2); degree(p2,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; complex_smod(q_orig,-i,modulo,q2); degree(q2,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; if (!mod_gcd(p1,q1,modulo,g1,pcof1,qcof1,vars,compute_cofactors,compute_cofactors,divtest,pv,qv,nthreads)) continue; // normalize g, pcof, qcof smallmult(smod(longlong(lg1)*invmod(g1.front().g,modulo),modulo),g1,g1,modulo); smallmult(smod(longlong(lp1)*invmod(pcof1.front().g,modulo),modulo),pcof1,pcof1,modulo); smallmult(smod(longlong(lq1)*invmod(qcof1.front().g,modulo),modulo),qcof1,qcof1,modulo); degree(g1,shift_vars,gdegmod); if (!mod_gcd(p2,q2,modulo,g2,pcof2,qcof2,vars,compute_cofactors,compute_cofactors,divtest,pv,qv,nthreads)) continue; smallmult(smod(longlong(lg2)*invmod(g2.front().g,modulo),modulo),g2,g2,modulo); smallmult(smod(longlong(lp2)*invmod(pcof2.front().g,modulo),modulo),pcof2,pcof2,modulo); smallmult(smod(longlong(lq2)*invmod(qcof2.front().g,modulo),modulo),qcof2,qcof2,modulo); degree(g2,shift_vars,gdegmod2); if (gdegmod2!=gdegmod) continue; int cmp=compare(gdegmod,gdeg); if (cmp==0){ // bad reduction continue; } if (cmp==-1){ // restart d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; continue; } if (cmp==-2){ // same deg, chinese remainder // same degrees, chinese remainder complex_ichinrem(g1,g2,i,modulo,d,pimod); if (compute_cofactors){ complex_ichinrem(pcof1,pcof2,i,modulo,pcofactor,pimod); complex_ichinrem(qcof1,qcof2,i,modulo,qcofactor,pimod); } pimod = modulo * pimod; } else { // restart with this gcd complex_unmod(g1,g2,d,i,modulo); if (compute_cofactors){ complex_unmod(pcof1,pcof2,pcofactor,i,modulo); complex_unmod(qcof1,qcof2,qcofactor,i,modulo); } pimod=modulo; gdeg=gdegmod; } // finished?? if (compute_cofactors){ // d*pcofactor=p_orig*gcdlcoeff mod pimod // if |gcdlcoeff*p_orig|+|g|*|pcof|*min(size(g),size(pcof))< pimod we are done int mgp=int(std::min(d.size(),pcofactor.size())); int mgq=int(std::min(d.size(),qcofactor.size())); gen maxg=max(d,context0),maxp=max(p_orig,context0),maxq=max(q_orig,context0),maxpcof=max(pcofactor,context0),maxqcof=max(qcofactor,context0); gen dz=ppz(d,0,false),pz=ppz(pcofactor,0,false),qz=ppz(qcofactor,0,false); maxg = maxg/abs(dz,context0); maxpcof = maxpcof/abs(pz,context0); maxqcof = maxqcof/abs(qz,context0); if (is_strictly_greater(pimod,abs(gcdlcoeff,context0)*maxp+mgp*maxg*maxpcof,context0) && is_strictly_greater(pimod,abs(gcdlcoeff,context0)*maxq+mgq*maxg*maxqcof,context0) ){ // divide g,pcofactor,qcofactor by their integer content ppz(d,dz,true); if (compute_cofactors){ ppz(pcofactor,pz,true); ppz(qcofactor,qz,true); } return true; } } if (divtest) { // division test vector< T_unsigned > dtest(d),pquo,qquo,rem; ppz(d,0,true); if (hashdivrem(p_orig,dtest,pquo,rem,vars,0 /* reduce */,0/*qmax*/,false)==1 && rem.empty()){ if (hashdivrem(q_orig,dtest,qquo,rem,vars,0 /* reduce */,0/*qmax*/,false)==1 && rem.empty()){ pcofactor=pquo; qcofactor=qquo; return true; } } d=dtest; } } // end for (;;) return false; } // end if (is_complex) // gcd in Z[x1,..,xn] vector< T_unsigned > p,q,g,pcof,qcof; for (;;){ m=nextprime(m+1); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || m.type!=_INT_) return false; int modulo=m.val; int lg=smod(gcdlcoeff,modulo).val,lp=smod(lcoeffp,modulo).val,lq=smod(lcoeffq,modulo).val; if (!lg || !lp || !lq) continue; smod(p_orig,modulo,p); degree(p,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; smod(q_orig,modulo,q); degree(q,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; if (!mod_gcd(p,q,modulo,g,pcof,qcof,vars,compute_cofactors,compute_cofactors,divtest,pv,qv,nthreads)) continue; degree(g,shift_vars,gdegmod); // normalize g, pcof, qcof smallmult(smod(longlong(lg)*invmod(g.front().g,modulo),modulo),g,g,modulo); smallmult(smod(longlong(lp)*invmod(pcof.front().g,modulo),modulo),pcof,pcof,modulo); smallmult(smod(longlong(lq)*invmod(qcof.front().g,modulo),modulo),qcof,qcof,modulo); // CERR << " g " << g << " pcof " << pcof << "qcof " << qcof << endl; int cmp=compare(gdegmod,gdeg); if (cmp==0){ // bad reduction continue; } if (cmp==-1){ // restart d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; continue; } if (cmp==-2){ // same degrees, chinese remainder ichinrem(g,modulo,d,pimod); if (compute_cofactors){ ichinrem(pcof,modulo,pcofactor,pimod); ichinrem(qcof,modulo,qcofactor,pimod); } pimod = modulo * pimod; } else { // restart with this modular gcd unmod(g,d,modulo); if (compute_cofactors){ unmod(pcof,pcofactor,modulo); unmod(qcof,qcofactor,modulo); } pimod=modulo; gdeg=gdegmod; } // are we finished? if (compute_cofactors){ // d*pcofactor=p_orig*gcdlcoeff mod pimod // if |gcdlcoeff*p_orig|+|g|*|pcof|*min(size(g),size(pcof))< pimod we are done int mgp=int(std::min(d.size(),pcofactor.size())); int mgq = int(std::min(d.size(), qcofactor.size())); gen maxg=max(d,context0),maxp=max(p_orig,context0),maxq=max(q_orig,context0),maxpcof=max(pcofactor,context0),maxqcof=max(qcofactor,context0); gen dz=ppz(d,0,false),pz=ppz(pcofactor,0,false),qz=ppz(qcofactor,0,false); maxg = maxg/ dz; maxpcof = maxpcof/pz; maxqcof = maxqcof/qz; if (is_strictly_greater(pimod,abs(gcdlcoeff,context0)*maxp+mgp*maxg*maxpcof,context0) && is_strictly_greater(pimod,abs(gcdlcoeff,context0)*maxq+mgq*maxg*maxqcof,context0) ){ // divide g,pcofactor,qcofactor by their integer content ppz(d,dz,true); if (compute_cofactors){ ppz(pcofactor,pz,true); ppz(qcofactor,qz,true); } return true; } } if (divtest) { // division test vector< T_unsigned > dtest(d),pquo,qquo,rem; ppz(d,0,true); if (hashdivrem(p_orig,dtest,pquo,rem,vars,0 /* reduce */,0/*qmax*/,false)==1 && rem.empty()){ if (hashdivrem(q_orig,dtest,qquo,rem,vars,0 /* reduce */,0/*qmax*/,false)==1 && rem.empty()){ pcofactor=pquo; qcofactor=qquo; return true; } } d=dtest; } } #endif // NO_TEMPLATE_MULTGCD } /* ************************************************* MODULAR ALGEBRAIC EXTENSIONS GCD ************************************************* */ void mulsmall(const vector::const_iterator & ita0,const vector::const_iterator & ita_end,const vector::const_iterator & itb0,const vector::const_iterator & itb_end,int modulo,vector & new_coord){ new_coord.clear(); if (ita0==ita_end || itb0==itb_end) return; new_coord.reserve((ita_end-ita0)+(itb_end-itb0)-1); vector::const_iterator ita_begin=ita0,ita=ita0,itb=itb0; for ( ; ita!=ita_end; ++ita ){ vector::const_iterator ita_cur=ita,itb_cur=itb; longlong res=0; for (;itb_cur!=itb_end;--ita_cur,++itb_cur) { res += longlong(*ita_cur) * *itb_cur ; if (ita_cur==ita_begin) break; } new_coord.push_back(smod(res,modulo)); } --ita; ++itb; for ( ; itb!=itb_end;++itb){ longlong res= 0; vector::const_iterator ita_cur=ita,itb_cur=itb; for (;;) { res += longlong(*ita_cur) * *itb_cur ; if (ita_cur==ita_begin) break; --ita_cur; ++itb_cur; if (itb_cur==itb_end) break; } new_coord.push_back(smod(res,modulo)); } } // res=a*b mod pmin,modulo void mulext(const vector & a,const vector & b,const vector & pmin,int modulo,vector & res){ if (b.empty()){ res.clear(); return; } if (b.size()==1 && b.front()==1){ res=a; return; } vector q,tmp; mulsmall(a.begin(),a.end(),b.begin(),b.end(),modulo,tmp); DivRem(tmp,pmin,modulo,q,res); } // res=a*b mod pmin,modulo with temporary passed as arg static void mulext(const vector & a,const vector & b,const vector & pmin,int modulo,vector & res,vector & tmp,vector & q){ if (b.empty()){ res.clear(); return; } if (b.size()==1 && b.front()==1){ res=a; return; } mulsmall(a.begin(),a.end(),b.begin(),b.end(),modulo,tmp); DivRem(tmp,pmin,modulo,q,res); } static gen mulextaux2(const gen & a,const gen & b,int modulo){ gen res= a*b; if (res.type==_FRAC) return smod(res,modulo); else return res; } static void mulextaux(const modpoly & a,const modpoly &b,int modulo,modpoly & new_coord){ #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted) { interrupted = true; ctrl_c=false; new_coord=vecteur(1, gensizeerr(gettext("Stopped by user interruption."))); return; } modpoly::const_iterator ita=a.begin(),ita_end=a.end(),itb=b.begin(),itb_end=b.end(); new_coord.clear(); ita=a.begin(); itb=b.begin(); if (ita==ita_end || itb==itb_end) return; modpoly::const_iterator ita_begin=ita; for ( ; ita!=ita_end; ++ita ){ modpoly::const_iterator ita_cur=ita,itb_cur=itb; gen res; for (;;) { res += mulextaux2(*ita_cur,*itb_cur,modulo); // res = res + (*ita_cur) * (*itb_cur); if (ita_cur==ita_begin) break; --ita_cur; ++itb_cur; if (itb_cur==itb_end) break; } new_coord.push_back(res); } --ita; ++itb; for ( ; itb!=itb_end;++itb){ modpoly::const_iterator ita_cur=ita,itb_cur=itb; gen res; for (;;) { res += mulextaux2(*ita_cur,*itb_cur,modulo); // res = res + (*ita_cur) * (*itb_cur); if (ita_cur==ita_begin) break; --ita_cur; ++itb_cur; if (itb_cur==itb_end) break; } new_coord.push_back(res); } } static void mulext(const vecteur & a,const vecteur & b,const vecteur & pmin,int modulo,vecteur & res){ if (b.empty()){ res.clear(); return; } if (b.size()==1 && b.front()==1){ res=a; return; } vecteur q,tmp; environment env; env.modulo=modulo; env.moduloon=true; mulextaux(a,b,modulo,tmp); // operator_times(a,b,0,tmp); // bug if a and b contains fraction DivRem(tmp,pmin,&env,q,res); } // a *=b mod pmin, modulo static void mulext(vector & a,const vector & b,const vector & pmin,int modulo){ if (b.empty()){ a.clear(); return; } if (b.size()==1 && b.front()==1) return; vector q,tmp; mulsmall(a.begin(),a.end(),b.begin(),b.end(),modulo,tmp); DivRem(tmp,pmin,modulo,q,a); } static void mulext(vector< vector > & a,const vector & b,const vector & pmin,int modulo,vector & tmp,vector & q){ if (b.empty()){ a.clear(); return; } if (b.size()==1 && b.front()==1) return; vector< vector >::iterator it=a.begin(),itend=a.end(); for (;it!=itend;++it){ mulsmall(it->begin(),it->end(),b.begin(),b.end(),modulo,tmp); DivRem(tmp,pmin,modulo,q,*it); } } // inverse of rootof(a:pmin) modulo // pmin=1*pmin+0*a // a_=0*a+1*b // ... // 1=?*a+ainv*b static bool invmodext_(const vector & a_,const vector & pmin,int modulo,vector & u0,vector &a,vector & b,vector &q,vector &r,vector &u1,vector & u2){ if (a_.empty()) return false; if (a_.size()==1){ if (gcd(modulo,a_[0])!=1) return false; u0.resize(1); u0[0]=invmod(a_[0],modulo); return true; } if (debug_infolevel>10) CERR << a_ << " inv " << pmin << " mod " << modulo << endl; a=pmin; b=a_; vector::iterator it,itend; u0.clear(); u1.push_back(1); for (;!b.empty();){ if (gcd(modulo,b.front())!=1) return false; DivRem(a,b,modulo,q,r); swap(a,b); // a,b,r -> b,a,r swap(b,r); // -> b,r,a // u2=u0-q*u1 modulo mulext(q,u1,pmin,modulo,u2); // u2=q*u1 submod(u2,u0,modulo); // u2=q*u1-u0 for (it=u2.begin(),itend=u2.end();it!=itend;++it) *it = -*it; // done swap(u0,u1); swap(u1,u2); if (debug_infolevel>10) CERR << a << " " << b << " " << u0 << " " << u1 << endl; } if (gcd(modulo,a.front())!=1 || u0.empty()) return false; mulmod(u0,invmod(a.front(),modulo),modulo); return true; } static bool invmodext(const vector & a_,const vector & pmin,int modulo,vector & u0){ vector a,b,q,r,u1,u2; // u0=ainv return invmodext_(a_,pmin,modulo,u0,a,b,q,r,u1,u2); } static vector invmod(const vector & a,const modred & R){ vector ainv; if (!invmodext(a,R.pmin,R.modulo,ainv)){ #ifndef NO_STDEXCEPT setsizeerr(gettext("invmodext")); #endif ainv.clear(); } return ainv; } // c = a*b % reduce static void type_operator_reduce(const vector & a,const vector & b,vector & c,const modred & reduce){ mulext(a,b,reduce.pmin,reduce.modulo,c); } vector operator %(const vector & a,const modred & reduce){ vector q,res; DivRem(a,reduce.pmin,reduce.modulo,q,res); return res; } // c += a*b % reduce static void type_operator_plus_times_reduce(const vector & a,const vector & b,vector & c,const modred & reduce){ vector tmp; mulext(a,b,reduce.pmin,reduce.modulo,tmp); addmod(c,tmp,reduce.modulo); } static vecteur invmod(const vecteur & a,const Modred & R){ if (a.size()==1) return vecteur(1,invmod(a.front(),R.modulo)); // should not occur if division test is done with lcoeff=1 vecteur ainv,v,d; environment env; env.modulo=R.modulo; env.moduloon=true; egcd(a,R.pmin,&env,ainv,v,d); if (d.size()!=1){ #ifndef NO_STDEXCEPT setsizeerr(gettext("invmodext")); #endif return vecteur(1,gensizeerr(gettext("invmodext"))); } return ainv; } // c = a*b % reduce static void type_operator_reduce(const vecteur & a,const vecteur & b,vecteur & c,const Modred & reduce){ mulext(a,b,reduce.pmin,reduce.modulo,c); } static void type_operator_times(const vecteur & a,const vecteur &b,vecteur & c){ mulmodpoly(a,b,c); } static void type_operator_plus_times(const vecteur & a,const vecteur & b,vecteur & c){ vecteur tmp; mulmodpoly(a,b,tmp); addmodpoly(c,tmp,c); } static void type_operator_plus_times_reduce(const vecteur & a,const vecteur & b,vecteur & c,int reduce){ assert(reduce==0); vecteur tmp; mulmodpoly(a,b,tmp); addmodpoly(c,tmp,c); } // c += a*b % reduce static void type_operator_plus_times_reduce(const vecteur & a,const vecteur & b,vecteur & c,const Modred & reduce){ vecteur tmp; mulext(a,b,reduce.pmin,reduce.modulo,tmp); addmod(c,tmp,reduce.modulo); } vecteur operator %(const vecteur & a,const Modred & reduce){ vecteur q,res; environment env; env.modulo=reduce.modulo; env.moduloon=true; DivRem(a,reduce.pmin,&env,q,res); return res; } static void smallmult(int g,const std::vector< T_unsigned,hashgcd_U> > & v1,std::vector< T_unsigned,hashgcd_U> > & v,int reduce){ if (!g){ v.clear(); return; } std::vector< T_unsigned,hashgcd_U> >::const_iterator it1=v1.begin(),it1end=v1.end(); if (&v1==&v){ std::vector< T_unsigned,hashgcd_U> >::iterator it1=v.begin(),it1end=v.end(); for (;it1!=it1end;++it1){ mulmod(g,it1->g,reduce); // it1->g = (g*it1->g) % reduce; } } else { v.clear(); v.reserve(it1end-it1); // worst case vector res; for (;it1!=it1end;++it1){ res=it1->g; mulmod(g,res,reduce); v.push_back(T_unsigned,hashgcd_U>(res,it1->u)); } } } // a=b*q+r, modifies b (multiplication by inverse of lcoeff) // if you want q corresponding to original b, set q_orig_b to true static bool divrem_(vector< vector > & a,vector< vector > & b,const vector & pmin, int modulo, vector< vector > * qptr,bool set_q_orig_b,vector & b0,vector & b0inv,vector & tmp,vector & tmp1,vector & tmp2,vector & tmp3,vector & tmp4,vector & tmp5){ int as=int(a.size()),bs=int(b.size()); if (!bs) return false; if (set_q_orig_b) b0=b.front(); if (!invmodext_(b.front(),pmin,modulo,b0inv,tmp,tmp1,tmp2,tmp3,tmp4,tmp5)) return false; if (qptr) qptr->clear(); if (as >::iterator it,jt,jtend=b.end(); int nstep=as-bs+1,pos=0; #if 0 if (nstep==2){ // normal remainder sequence // compute quotient from a[0], a[1] and b[0]==1 b[1] // x coeff of q is q1=a[0], cst coeff of q is q0=a[1]-a[0]*b[1] // a -= q*b // x^k coeff ak -> ak- q0*bk - q1*bk-1, but beware a[0]==lcoeff(a) // a[k] = a[k]- (q1*b[k] + q0*b[k-1]) mulext(a[0],b[1],pmin,modulo,tmp,tmp1,tmp2); submod(a[1],tmp,modulo); // store q0 in a[1] (a[1] not used anymore) for (int k=2;k & q = a[pos]; if (qptr){ if (set_q_orig_b){ mulext(q,b0inv,pmin,modulo,tmp,tmp1,tmp2); qptr->push_back(tmp); } else qptr->push_back(q); } if (!q.empty()){ for (it=a.begin()+pos+1,jt=b.begin()+1;jt!=jtend;++it,++jt){ mulext(*jt,q,pmin,modulo,tmp,tmp1,tmp2); submod(*it,tmp,modulo); } } } while (pos > & p,const vector< vector > & q,const vector & pmin,int modulo,vector< vector > & d){ // lcoeffs of p/q should be invertible mod modulo // for q this is checked at the first divrem if (p.empty()){ d=q; return true; } if (q.empty()){ d=p; return true; } vector p0inv,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8; if (!invmodext(p.front(),pmin,modulo,p0inv)) return false; d=p; vector< vector > b(q); for (;!b.empty();){ if (!divrem_(d,b,pmin,modulo,0,false,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8)) return false; swap(d,b); } if (!invmodext(d.front(),pmin,modulo,p0inv)) return false; mulext(d,p0inv,pmin,modulo,tmp1,tmp2); return true; } static bool smod_ext(const vector< T_unsigned > & p_orig,int modulo,vector< T_unsigned > & p){ p.clear(); vector< T_unsigned >::const_iterator it=p_orig.begin(),itend=p_orig.end(); p.reserve(itend-it); vecteur x; gen tmp,itg; const_iterateur vt,vtend; for (;it!=itend;++it){ x.clear(); if (it->g.type!=_EXT){ tmp=smod(it->g,modulo); if (!is_zero(tmp)) x.push_back(tmp); } else { tmp=*it->g._EXTptr; if (tmp.type!=_VECT || tmp._VECTptr->empty()) return false; vt=tmp._VECTptr->begin(),vtend=tmp._VECTptr->end(); x.reserve(vtend-vt); for (;vt!=vtend;++vt){ itg=smod(*vt,modulo); if (is_zero(itg) && x.empty()) continue; x.push_back(itg); } } if (!x.empty()) p.push_back( T_unsigned (x,it->u) ); } return true; } static vector hornermod_ext(const vector< vector > & v,int alpha,int modulo){ vector< vector >::const_iterator it=v.begin(),itend=v.end(); if (!alpha) return it==itend?vector(0):v.back(); vector res; for (;it!=itend;++it){ mulmod(res,alpha,modulo); addmod(res,*it,modulo); } return res; } // convert one dimensional polynomial with ext coeff to full dense representation // back conversion is convert_back (template coefficients) static void convert(const vector< T_unsigned,hashgcd_U> > & p,hashgcd_U var,vector< vector > & v){ v.clear(); vector< T_unsigned,hashgcd_U> >::const_iterator it=p.begin(),itend=p.end(); if (it==itend) return; hashgcd_U u=it->u; unsigned s=u/var; v=vector< vector >(s+1); for (;it!=itend;++it){ v[s-it->u/var]=it->g; } } // find n and d such that n/d=a mod m, modulo // and deg(n) <= deg(m)/2, deg(d) <= deg(m)-deg(a) // n=m*u+a*d, // Bezout for m and a, m*u_k + a*v_k = r_k, stop when deg(r_k)<=deg(m)/2 // n=r_k and d=v_k, then check that d is prime with m and n static bool polyfracmod(const vector & a,const vector &m,int modulo,vector & n, vector & d){ vector r0(m),r1(a),r2,v0,v1(1,1),v2,q,tmp; // m*u0+a*v0=m=r0, m*u1+a*v1=a=r1, u0, u1 and u2 are not computed (not used) int N=int(m.size()-1)/2; for (;(int)r1.size()>N+1;){ DivRem(r0,r1,modulo,q,r2); smallmult(q.begin(),q.end(),v1.begin(),v1.end(),tmp,modulo); swap(v0,v2); // v2=v0; but v0 is not used anymore submod(v2,tmp,modulo); // move r1 to r0 and r2 to r1, same with v // A B C swap(r0,r1); // B A C swap(r1,r2); // B C A // unused B C swap(v0,v1); // B unused C swap(v1,v2); // B C unused } // deg(r1)<=N, store answer swap(r1,n); swap(v1,d); // checks gcdsmallmodpoly(d,m,modulo,r2); if (r2.size()>1) return false; gcdsmallmodpoly(n,d,modulo,r2); if (r2.size()>1) return false; return true; } static bool polyfracmod(const gen & g,const vector & m,int modulo,gen & res){ if (g.type==_POLY){ polynome p(g._POLYptr->dim); vector >::const_iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); for (;it!=itend;++it){ const gen tmp=it->value; gen x; if (tmp.type!=_VECT) p.coord.push_back(monomial(tmp,it->index)); else { polyfracmod(*tmp._VECTptr,m,modulo,x); p.coord.push_back(monomial(x,it->index)); } } res=p; return true; } if (g.type!=_VECT){ res=g; return true; } if (g._VECTptr->size()-1<=(m.size()-1)/2){ res=g; return true; } // FIXME coeffs might be polynomials instead of integers vector a,n,d; vecteur N,D; vecteur2vector_int(*g._VECTptr,modulo,a); if (!polyfracmod(a,m,modulo,n,d)) return false; vector_int2vecteur(n,N); vector_int2vecteur(d,D); res=fraction(N,D); return true; } static bool polyfracmod(vecteur & res,const vector & m,int modulo){ iterateur it=res.begin(),itend=res.end(); for (;it!=itend;++it){ if (!polyfracmod(*it,m,modulo,*it)) return false; } return true; } // replace p by lcm of p and g static bool vlcm(modpoly & p,const gen & g,int modulo){ if (g.type!=_VECT) return true; const modpoly & q=*g._VECTptr; environment env; env.modulo=modulo; env.moduloon=true; modpoly quo,rem,pgcd; // first check if g divides p if (!DivRem(p,q,&env,quo,rem,false)) return false; // setsizeerr(); if (rem.empty()) return true; // nope, compute gcd and divide p*q by gcd gcdmodpoly(p,q,&env,pgcd); if (is_undef(pgcd)) return false; DivRem(p,pgcd,&env,quo,rem); mulmodpoly(q,quo,&env,p); return true; } // rational recon and clear fraction, return size of lcm of denom static int polyfracmod(const vector< T_unsigned > & p0,const vector & m,int modulo,vector< T_unsigned > & p,int maxdeg){ p=p0; vector< T_unsigned >::iterator kt=p.begin(),ktend=p.end(); for (;kt!=ktend;++kt){ if (!polyfracmod(kt->g,m,modulo)) return 0; } // find lcm of deno vecteur ppcm(1,1); for (kt=p.begin();kt!=ktend;++kt){ if ((int)ppcm.size()>maxdeg+1) return 0; iterateur it=kt->g.begin(),itend=kt->g.end(); for (;it!=itend;++it){ if (it->type==_FRAC){ if (!vlcm(ppcm,it->_FRACptr->den,modulo)) return 0; if ((int)ppcm.size()>maxdeg+1) return 0; } else { if (it->type==_POLY){ vector< monomial >::const_iterator jt=it->_POLYptr->coord.begin(),jtend=it->_POLYptr->coord.end(); for (;jt!=jtend;++jt){ if (jt->value.type==_FRAC){ if (!vlcm(ppcm,jt->value._FRACptr->den,modulo)) return 0; if ((int)ppcm.size()>maxdeg+1) return 0; } } } } } } // clear deno gen den(ppcm,_POLY1__VECT); environment env; env.modulo=modulo; env.moduloon=true; if (ppcm.size()>1){ gen adjust=invmod(ppcm.front(),modulo); if (is_undef(adjust)) return 0; for (kt=p.begin();kt!=ktend;++kt){ iterateur it=kt->g.begin(),itend=kt->g.end(); for (;it!=itend;++it){ if (it->type==_FRAC){ if (it->_FRACptr->den.type!=_VECT) return 0; // setsizeerr(); modpoly quot,rest; if (!DivRem(ppcm,*it->_FRACptr->den._VECTptr,&env,quot,rest,false) || !rest.empty()) return 0; // setsizeerr(); *it = smod(adjust*(it->_FRACptr->num * gen(quot,_POLY1__VECT)),modulo); } else { if (it->type==_POLY){ vector< monomial >::iterator jt=it->_POLYptr->coord.begin(),jtend=it->_POLYptr->coord.end(); for (;jt!=jtend;++jt){ if (jt->value.type==_FRAC){ modpoly quot,rest; if (!DivRem(ppcm,*jt->value._FRACptr->den._VECTptr,&env,quot,rest,false) || !rest.empty()) return 0; // setsizeerr(); jt->value = smod(adjust*(jt->value._FRACptr->num * gen(quot,_POLY1__VECT)),modulo); } else jt->value = smod(adjust*(den * jt->value),modulo) ; } } else *it = smod(adjust*den * (*it),modulo); } } } } return int(ppcm.size()); } static bool gentoint(const vector< T_unsigned > & p0,vector< T_unsigned,hashgcd_U> > & p){ vector< T_unsigned >::const_iterator it=p0.begin(),itend=p0.end(); p.clear(); p.reserve(itend-it); for (;it!=itend;++it){ vecteur::const_iterator jt=it->g.begin(),jtend=it->g.end(); p.push_back(T_unsigned,hashgcd_U>(vector(0),it->u)); vector & res=p.back().g; for (;jt!=jtend;++jt){ if (jt->type!=_INT_) return false; res.push_back(jt->val); } } return true; } static void inttogen(const vector< T_unsigned,hashgcd_U> > & p0,vector< T_unsigned > & p){ vector< T_unsigned,hashgcd_U> >::const_iterator it=p0.begin(),itend=p0.end(); p.clear(); p.reserve(itend-it); for (;it!=itend;++it){ vector::const_iterator jt=it->g.begin(),jtend=it->g.end(); p.push_back(T_unsigned(vecteur(0),it->u)); vecteur & res=p.back().g; for (;jt!=jtend;++jt){ res.push_back(*jt); } } } // eval last variable of polynomials in g at x modulo modulo, answer in res // polynomial coeff must be integers static bool horner_back_ext(const gen & g,int x,int modulo,gen & res){ if (g.type!=_POLY){ if (g.type==_FRAC){ gen r1,r2; if (!horner_back_ext(g._FRACptr->num,x,modulo,r1) || !horner_back_ext(g._FRACptr->den,x,modulo,r2)) return false; res=r1/r2; return true; } if (g.type==_EXT){ gen r1,r2; if (!horner_back_ext(*g._EXTptr,x,modulo,r1) || !horner_back_ext(*(g._EXTptr+1),x,modulo,r2)) return false; res=algebraic_EXTension(r1,r2); return true; } if (g.type!=_VECT){ res=g; return true; } res=gen(*g._VECTptr,g.subtype); if (res.type!=_VECT) return true; gen tmp; iterateur it=res._VECTptr->begin(),itend=res._VECTptr->end(); for (;it!=itend;++it){ if (!horner_back_ext(*it,x,modulo,tmp)) return false; *it=tmp; } return true; } const polynome & p = *g._POLYptr; int dim=p.dim; polynome r(dim-1); vector< monomial >::const_iterator it=p.coord.begin(),itend=p.coord.end(); if (it==itend){ res=0; return true; } index_t current(it->index.begin(),it->index.begin()+dim-1),nouveau(dim-1); int expo=it->index.back(),newexpo; if (it->value.type!=_INT_) return false; // setsizeerr(); int val=it->value.val; for (++it;it!=itend;++it){ // IMPROVE: for dense poly: skip to it+expo check if same outer exponents if (it->value.type!=_INT_) return false; // setsizeerr(); index_t::iterator nit=nouveau.begin(); index_t::const_iterator jt=it->index.begin(),jtend=it->index.begin()+dim-1; for (;jt!=jtend;++nit,++jt){ *nit=*jt; } if (nouveau==current){ val=(powmod(x,expo-(newexpo=it->index.back()),modulo)*longlong(val)+it->value.val)%modulo; expo=newexpo; continue; } if (expo) val=(powmod(x,expo,modulo)*longlong(val))%modulo; r.coord.push_back(monomial(val,current)); swap(current,nouveau); val=it->value.val; expo=it->index.back(); } if (expo) val=(powmod(x,expo,modulo)*longlong(val))%modulo; r.coord.push_back(monomial(val,current)); res=r; return true; } // return false if the lcoeff is 0 static bool horner_back_ext2(const vector< T_unsigned > & g,int x,int modulo,vector< T_unsigned > & res){ gen tmp; res.clear(); vector< T_unsigned >::const_iterator it=g.begin(),itend=g.end(); bool b=true; for (;it!=itend;++it){ if (!horner_back_ext(it->g,x,modulo,tmp)) return false; if (is_zero(tmp)){ if (it==res.begin()) b=false; continue; } if (tmp.type==_VECT) res.push_back(T_unsigned(*tmp._VECTptr,it->u)); else res.push_back(T_unsigned(vecteur(1,tmp),it->u)); } return b; } struct ext_gcd_t { gen pi_t_minus_tk; vector< T_unsigned > lagrange; // the vecteur represents dependency wrt the extension variable (z) // each coeff depends on t1..tk-1 variable as a _POLY // the depency wrt tk is represented by coeffs=vecteur of type poly1[] }; static gen hornermod(const gen & g,int tk,int modulo){ if (g.type==_POLY){ polynome res(g._POLYptr->dim); vector< monomial >::const_iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); for (;it!=itend;++it){ gen tmp=hornermod(it->value,tk,modulo); if (!is_zero(tmp)) res.coord.push_back(monomial(tmp,it->index)); } return res; } if (g.type!=_VECT) return g; int res=0; const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); for (;it!=itend;++it){ res = ( res*longlong(tk)+it->val)%modulo; } return res; } static vecteur interp_tk(const vecteur & interp_ancien,const gen & ancienpi,const vecteur & interp_tk,int tk,int modulo){ int as=int(interp_ancien.size()),bs=int(interp_tk.size()),cs=giacmax(as,bs); vecteur res(cs); gen pi(hornermod(ancienpi,tk,modulo)); if (pi.type!=_INT_) return vecteur(1,gensizeerr(gettext("interp_tk"))); int invpi=invmod(pi.val,modulo); for (int i=0;i=cs) old=hornermod(interp_ancien[i+as-cs],tk,modulo); if (i+bs>=cs) nouv=interp_tk[i+bs-cs]; gen alpha=smod(invpi*(nouv-old),modulo); // multiply by ancienpi and add to interp_ancien if (alpha.type==_POLY) alpha=ancienpi * (*alpha._POLYptr); else alpha=ancienpi*alpha; if (i+as>=cs){ gen g =interp_ancien[i+as-cs]; if (g.type==_POLY && alpha.type!=_POLY) res[i]=smod(addpoly(*g._POLYptr,alpha),modulo); else { if (g.type!=_POLY && alpha.type==_POLY) res[i]=smod(addpoly(*alpha._POLYptr,g),modulo); else res[i]=smod(g+alpha,modulo); } } else res[i]=smod(alpha,modulo); } return res; } static void untrunc(gen & g,int innerdim){ if (g.type==_VECT){ polynome p(innerdim); index_t i(innerdim); const_iterateur jt=g._VECTptr->begin(),jtend=g._VECTptr->end(); if (jtend==jt){ g=0; return; } i.back()=int(jtend-jt)-1; for (;jt!=jtend;--i.back(),++jt){ if (!is_zero(*jt)) p.coord.push_back(monomial(*jt,i)); } g=p; return; } if (g.type!=_POLY) return; polynome p(innerdim); vector< monomial >::const_iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); for (;it!=itend;++it){ index_t i(it->index.iref()); i.push_back(0); if (it->value.type!=_VECT) p.coord.push_back(monomial(it->value,i)); else { const_iterateur jt=it->value._VECTptr->begin(),jtend=it->value._VECTptr->end(); if (jtend==jt) continue; i.back()=int(jtend-jt)-1; for (;jt!=jtend;--i.back(),++jt){ if (!is_zero(*jt)) p.coord.push_back(monomial(*jt,i)); } } } g=p; } static void untrunc(vector< T_unsigned > & p,int dim){ vector< T_unsigned >::iterator it=p.begin(),itend=p.end(); for (;it!=itend;++it){ iterateur jt=it->g.begin(),jtend=it->g.end(); for (;jt!=jtend;++jt) untrunc(*jt,dim); } } static int gcd_ext(const vector< T_unsigned > & p0,const vector< T_unsigned > & q0,const std::vector & vars,const vecteur & Pmin,int modulo, vector< T_unsigned > & D, vector< T_unsigned > & Pcof,vector< T_unsigned > & Qcof,bool compute_pcof,bool compute_qcof, int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return 0; #else // FIXME cofactors are not yet implemented gen extension; const_iterateur it=Pmin.begin(),itend=Pmin.end(); for (;it!=itend;++it){ if (it->type==_POLY) extension=*it; else { if (!it->is_integer()) return 0; } } if (extension.type==_POLY){ // return -1; // avoid inf loop while not implemented Modred pminmodulo(modulo,Pmin); // otherwise give ext variables values, subroutine P van Hoeij-Monagan // Algo. for Polynomial GCD computation over Algebraic Function Fields int d=1,n=1,innerdim=extension._POLYptr->dim; if (innerdim<1) return 0; // setsizeerr(); map,ext_gcd_t> m; // index is obtained by taking lcoeff with respect to outer var // (hashgcd_U) and inner var (index_t) vector< T_unsigned > p0tk,q0tk,dtk,pcoftk,qcoftk,nouveau,test,prevtest,tmprem; for (int tcount=0;;++tcount){ int tk=giac_rand(context0) % modulo; gen pmintk; if (!horner_back_ext(Pmin,tk,modulo,pmintk)) return 0; if (pmintk.type!=_VECT) continue; vecteur & Pmintk = *pmintk._VECTptr; if (is_zero(Pmintk.front())) continue; if (!horner_back_ext2(p0,tk,modulo,p0tk)) continue; if (!horner_back_ext2(q0,tk,modulo,q0tk)) continue; int res=gcd_ext(p0tk,q0tk,vars,Pmintk,modulo,dtk,pcoftk,qcoftk,false,false,nthreads); if (res==-1) return -1; if (res!=1){ ++d; if (d>n) return 0; else continue; } // dtk lcoeff should not depend on extension variable if (dtk.empty()) return 0; // setsizeerr(); vecteur & lcoeffdtk = dtk.front().g; if (lcoeffdtk.size()!=1) return 0; // setsizeerr(); index_t innerlcoeff(innerdim-1); if (lcoeffdtk.front().type==_POLY && !lcoeffdtk.front()._POLYptr->coord.empty()) innerlcoeff=lcoeffdtk.front()._POLYptr->coord.front().index.iref(); // adjust m[pair(dtk.front().u,innerlcoeff)] pair pa(dtk.front().u,innerlcoeff); map,ext_gcd_t>::iterator jt=m.find(pa),jtend=m.end(); if (jt==jtend){ m[pa].pi_t_minus_tk=gen(makevecteur(1,-tk),_POLY1__VECT); m[pa].lagrange=dtk; } else { gen ancienpi=m[pa].pi_t_minus_tk; vector< T_unsigned >::iterator jt=m[pa].lagrange.begin(),jtend=m[pa].lagrange.end(),kt=dtk.begin(),ktend=dtk.end(); nouveau.clear(); for (;;){ if (jt==jtend && kt==ktend) break; if (jt==jtend || jt->u>kt->u){ nouveau.push_back(T_unsigned(interp_tk(vecteur(0),ancienpi,kt->g,tk,modulo),kt->u)); ++kt; continue; } if (kt==ktend || jt->uu){ nouveau.push_back(T_unsigned(interp_tk(jt->g,ancienpi,vecteur(0),tk,modulo),jt->u)); ++jt; continue; } nouveau.push_back(T_unsigned(interp_tk(jt->g,ancienpi,kt->g,tk,modulo),kt->u)); ++jt; ++kt; } swap(m[pa].lagrange,nouveau); m[pa].pi_t_minus_tk=smod(ancienpi*gen(makevecteur(1,-tk),_POLY1__VECT),modulo); } // rational reconstruction for m[pa].lagrange modulo m[pa].pi_tk_minus_t // must be done on each component of the vecteur // which should be either a poly1[] or a _POLY with poly1[] coeffs // Clear fractions dependency wrt tk variable vector mo; vecteur2vector_int(*m[pa].pi_t_minus_tk._VECTptr,modulo,mo); if (!polyfracmod(m[pa].lagrange,mo,modulo,test,tcount)) continue; // now replace inner tk dependency as a poly1[] to an usual polynome untrunc(test,innerdim); // ?CHECK? clean lcoeff so that it is 1 wrt extension variable z // trial division, if success return 1 else continue // ?FIXME? use pseudo-division test if (debug_infolevel) CERR << CLOCK() << " algmodgcd hashdivrem " << test.size() <,hashgcd_U> > p_orig,q_orig,d,pcof,qcof; if (!gentoint(p0,p_orig) || !gentoint(q0,q_orig)) return 0; // ok, ext of dim 0 vector pmin; vecteur2vector_int(Pmin,modulo,pmin); int res=mod_gcd_ext(p_orig,q_orig,vars,pmin,modulo,d,pcof,qcof,compute_pcof,compute_qcof,nthreads); // convert back inttogen(d,D); if (compute_pcof) inttogen(pcof,Pcof); if (compute_qcof) inttogen(qcof,Qcof); return res; #endif // NO_TEMPLATE_MULTGCD } vector operator / (const vector & v,const vector & b){ #ifndef NO_STDEXCEPT setsizeerr(gettext("vector operator /")); #endif return v; } vector operator % (const vector & v,const vector & b){ #ifndef NO_STDEXCEPT setsizeerr(gettext("vector operator %")); #endif return v; } bool operator > (const vector & v,double q){ #ifndef NO_STDEXCEPT setsizeerr(gettext("vector operator >")); #endif return false; } static std::vector trim(const std::vector & res){ // trim res std::vector::const_iterator ita=res.begin(),itaend=res.end(); for (;ita!=itaend;++ita){ if (*ita) break; } return vector(ita,itaend); } std::vector operator + (const std::vector & a, const std::vector & b){ std::vector::const_iterator ita=a.begin(),itaend=a.end(),itb=b.begin(),itbend=b.end(); unsigned s=unsigned(itaend-ita),t=unsigned(itbend-itb); if (s>=t){ std::vector res(a); std::vector::iterator itres=res.begin()+(s-t); for (;itb!=itbend;++itb,++itres) *itres += (*itb); if (res.empty() || res.front()) return res; return trim(res); } std::vector res(b); std::vector::iterator itres=res.begin()+(t-s); for (;ita!=itaend;++ita,++itres) *itres += (*ita); return res; } std::vector operator - (const std::vector & a, const std::vector & b){ std::vector::const_iterator ita=a.begin(),itaend=a.end(),itb=b.begin(),itbend=b.end(); unsigned s=unsigned(itaend-ita),t=unsigned(itbend-itb); if (s>=t){ std::vector res(a); std::vector::iterator itres=res.begin()+(s-t); for (;itb!=itbend;++itb,++itres) *itres -= (*itb); if (res.empty() || res.front()) return res; // trim res return trim(res); } std::vector res(b); std::vector::iterator itres=res.begin(); for (;t>s;--t,++itres) *itres = -*itres; for (;ita!=itaend;++ita,++itres) *itres = *ita - *itres; return res; } std::vector operator - (const std::vector & a){ std::vector res(a); std::vector::iterator ita=res.begin(),itaend=res.end(); for (;ita!=itaend;++ita) *ita = -*ita; return res; } std::vector operator % (const std::vector & a,int modulo){ std::vector res(a); std::vector::iterator ita=res.begin(),itaend=res.end(); for (;ita!=itaend;++ita) *ita %= modulo; if (res.empty() || res.front()) return res; // trim res return trim(res); } static void distmult_ext(const vector< T_unsigned,hashgcd_U> > & p,const vector< vector > & v,vector< T_unsigned,hashgcd_U> > & pv,hashgcd_U var,const vector & pmin,int modulo){ if (&pv==&p){ vector< T_unsigned,hashgcd_U> > tmp; distmult_ext(p,v,tmp,var,pmin,modulo); swap(pv,tmp); return; } pv.clear(); vector< T_unsigned,hashgcd_U> >::const_iterator it=p.begin(),itend=p.end(); vector< vector >::const_iterator jtbeg=v.begin(),jtend=v.end(),jt; int vs=int(jtend-jtbeg),j; pv.reserve((itend-it)*vs); for (;it!=itend;++it){ for (jt=jtbeg,j=1;jt!=jtend;++j,++jt){ if (!is_zero(*jt)){ vector tmp; mulext(it->g,*jt,pmin,modulo,tmp); pv.push_back( T_unsigned< vector,hashgcd_U >(tmp,it->u+(vs-j)*var)); } } } } static bool horner(const vector< T_unsigned,hashgcd_U> > & p,int x,const std::vector & vars,vector< T_unsigned,hashgcd_U> > & px,int modulo,int maxdeg){ px.clear(); hashgcd_U var=vars[vars.size()-2]; hashgcd_U var2=vars.back(); vector< T_unsigned,hashgcd_U> >::const_iterator it=p.begin(),itend=p.end(),it1; vector::const_iterator jtend=vars.end()-1; hashgcd_U ucur,uend; if (maxdeg>=0){ uend=(maxdeg+1)*vars.front(); // dichotomy to find start position int pos1=0,pos2=int(itend-it),pos; for (;pos2-pos1>1;){ pos=(pos1+pos2)/2; if ((it+pos)->uu>=uend) ++it; } if (x==0){ for (;it!=itend;){ ucur=it->u; uend=(ucur/var)*var; if (ucur==uend){ const vector & g=smod(it->g,modulo); if (!is_zero(g)) px.push_back(T_unsigned,hashgcd_U>(g,uend)); ++it; continue; } register int nterms = (ucur-uend)/var2; if (ntermsu==uend){ it += nterms; const vector & g=smod(it->g,modulo); if (!is_zero(g)) px.push_back(T_unsigned,hashgcd_U>(g,uend)); ++it; continue; } for (++it;it!=itend;++it){ if (it->u<=uend){ if (it->u==uend){ const vector & g=smod(it->g,modulo); if (!is_zero(g)) px.push_back(T_unsigned,hashgcd_U>(g,uend)); ++it; } break; } } } return true; } for (;it!=itend;){ ucur=it->u; uend=(ucur/var)*var; if (ucur==uend){ px.push_back(*it); ++it; continue; } vector g(0); int nterms=(ucur-uend)/var2+1; // Check if the next group of monomials is dense wrt to xn it1=it+nterms; if (//false && ntermsu==uend ){ for (;it!=it1;++it){ // g = (g*x+it->g)%modulo; mulmod(g,x,modulo); addmod(g,it->g,modulo); } g=smod(g,modulo); if (!is_zero(g)) px.push_back(T_unsigned,hashgcd_U>(g,uend)); continue; } for (;it!=itend;++it){ const hashgcd_U & u=it->u; if (u,hashgcd_U>(g,uend)); } break; } mulmod(g,powmod(x,(ucur-u)/var2,modulo),modulo); addmod(g,it->g,modulo); ucur=u; } // end for if (it==itend){ if (!is_zero(g)){ mulmod(g,powmod(x,(ucur-uend)/var2,modulo),modulo); g=smod(g,modulo); if (!is_zero(g)) px.push_back(T_unsigned,hashgcd_U>(g,uend)); } } } return true; } #ifndef NO_TEMPLATE_MULTGCD static void * do_recursive_gcd_ext_call(void * ptr_){ #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted) return 0; gcd_call_param< vector > * ptr = (gcd_call_param< vector > *) ptr_; vector< vector > & Delta = *ptr->Delta; vector< vector > & lcoeffp = *ptr->lcoeffp; vector< vector > & lcoeffq = *ptr->lcoeffq; vector & alphav = * ptr->alphav; // vector< vector< vector > > & pv = *ptr->pv; // vector< vector< vector > > & qv = *ptr->qv; vector< vector > dim2palpha; vector< vector > dim2qalpha; // vector< vector< vector > > & dim2gcdv = *ptr->dim2gcdv; // vector< vector< vector > > & dim2pcofactorv = *ptr->dim2pcofactorv; // vector< vector< vector > > & dim2qcofactorv = *ptr->dim2qcofactorv; vector< T_unsigned< vector ,hashgcd_U> > & p = * ptr->p; vector< T_unsigned< vector ,hashgcd_U> > & q = * ptr->q; vector< vector< T_unsigned< vector ,hashgcd_U> > > & gcdv = * ptr->gcdv; vector< vector< T_unsigned< vector ,hashgcd_U> > > & pcofactorv = * ptr->pcofactorv; vector< vector< T_unsigned< vector ,hashgcd_U> > > & qcofactorv = * ptr->qcofactorv; index_t & pdeg = * ptr->pdeg; index_t & qdeg = * ptr->qdeg; vector< T_unsigned< vector ,hashgcd_U> > palpha,qalpha; index_t pdegalpha,qdegalpha; int vpos=ptr->vpos; int alpha1 = alphav[vpos]; int modulo = ptr->modulo; const vector & pmin = *ptr->pminptr; const vector & vars = * ptr->vars; vector & vars_truncated = * ptr->vars_truncated; // index_t & shift_vars = *ptr->shift_vars; index_t & shift_vars_truncated = *ptr->shift_vars_truncated; bool compute_cof = ptr->compute_cof; bool compute_qcofactor = ptr->compute_qcofactor; // bool dim2 = ptr->dim2; int nthreads = ptr->nthreads; // Eval p and q at xn=alpha if (!horner(p,alpha1,vars,palpha,modulo,-1)) return 0; degree(palpha,shift_vars_truncated,pdegalpha); pdegalpha.push_back(pdeg.back()); if (pdegalpha!=pdeg) return 0; if (!horner(q,alpha1,vars,qalpha,modulo,-1)) return 0; degree(qalpha,shift_vars_truncated,qdegalpha); qdegalpha.push_back(qdeg.back()); if (qdegalpha!=qdeg) return 0; vector< T_unsigned< vector ,hashgcd_U> > & g=gcdv[vpos]; vector< T_unsigned< vector ,hashgcd_U> > & gp=pcofactorv[vpos]; vector< T_unsigned< vector ,hashgcd_U> > & gq=qcofactorv[vpos]; if ( (ptr->ext_gcd_ok=mod_gcd_ext(palpha,qalpha,vars_truncated,pmin,modulo,g,gp,gq,compute_cof,compute_qcofactor,nthreads))!=1){ g.clear(); return ptr; } vector tmp; // adjust lcoeff of gcd to be the same as Delta // smallmult(smod(hornermod(Delta,alpha1,modulo)*longlong(invmod(g.front().g,modulo)),modulo),g,g,modulo); if (!invmodext(g.front().g,pmin,modulo,tmp)){ ptr->ext_gcd_ok=-1; return ptr; } mulext(tmp,hornermod_ext(Delta,alpha1,modulo),pmin,modulo); smallmult(tmp,g,g,modred(modulo,pmin)); if (compute_cof){ // adjust gp lcoeff // smallmult(smod(hornermod(lcoeffp,alpha1,modulo)*longlong(invmod(gp.front().g,modulo)),modulo),gp,gp,modulo); if (!invmodext(gp.front().g,pmin,modulo,tmp)){ ptr->ext_gcd_ok=-1; return ptr; } mulext(tmp,hornermod_ext(lcoeffp,alpha1,modulo),pmin,modulo); smallmult(tmp,gp,gp,modred(modulo,pmin)); if (compute_qcofactor){ // adjust gq cofactor // smallmult(smod(hornermod(lcoeffq,alpha1,modulo)*longlong(invmod(gq.front().g,modulo)),modulo),gq,gq,modulo); if (!invmodext(gq.front().g,pmin,modulo,tmp)){ ptr->ext_gcd_ok=-1; return ptr; } mulext(tmp,hornermod_ext(lcoeffq,alpha1,modulo),pmin,modulo); smallmult(tmp,gq,gq,modred(modulo,pmin)); } } return ptr; } #endif //NO_TEMPLATE_MULTGCD static int mod_gcd_ext(const vector< T_unsigned,hashgcd_U> > & p_orig,const vector< T_unsigned,hashgcd_U> > & q_orig,const std::vector & vars,const vector & pmin,int modulo, vector< T_unsigned,hashgcd_U> > & d, vector< T_unsigned,hashgcd_U> > & pcof,vector< T_unsigned,hashgcd_U> > & qcof,bool compute_pcofactor,bool compute_qcofactor, int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return 0; #else //nthreads=1; // FIXME, !=1 segfaults on compaq mini hashgcd_U var=vars.front(); int dim=int(vars.size()); if (dim==1){ vector< vector > P,Q,D,cof,tmp; vector tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8; convert(p_orig,var,P); convert(q_orig,var,Q); if (!gcdsmallmodpoly_ext(P,Q,pmin,modulo,D)) return 0; // should compute pcof and qcof if compute_p/qcofactor true if (compute_pcofactor){ if (!divrem_(P,D,pmin,modulo,&cof,true,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8)) return 0; convert_back(cof,var,pcof); } if (compute_qcofactor){ if (!divrem_(Q,D,pmin,modulo,&cof,true,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8)) return 0; convert_back(cof,var,qcof); } convert_back(D,var,d); return 1; } modred pminmodulo(modulo,pmin); // multi-dim gcd where pmin does not depend on vars -> Brown's recursive P algorithm // normalization: leading coeff in Z = 1 std::vector vars_truncated(vars); vars_truncated.pop_back(); hashgcd_U varxn=vars_truncated.back(),var2=vars.back(); index_t shift_vars,shift_vars_truncated; if (!find_shift(vars,shift_vars)) return 0; // setsizeerr(); shift_vars_truncated=shift_vars; shift_vars_truncated.pop_back(); short int shiftxn=shift_vars_truncated.back(),shift2=shift_vars.back(); // Make p and q primitive as polynomials in x1,...,xn-1 // with coeff polynomial in xn vector< T_unsigned,hashgcd_U> > p(p_orig),q(q_orig),pcont,qcont,dcont,tmp,pcofactor,qcofactor; vector< vector > pcontxn,qcontxn,dcontxn,pcofcontxn,qcofcontxn; if (debug_infolevel>20-dim) CERR << "gcdmod_ext threads " << nthreads << " content begin " << "dim " << dim << " " << CLOCK() << endl; if (!pp_mod(p,&pmin,modulo,varxn,var2,pcontxn)) return 0; if (!pp_mod(q,&pmin,modulo,varxn,var2,qcontxn)) return 0; if (!gcdsmallmodpoly_ext(pcontxn,qcontxn,pmin,modulo,dcontxn)) return 0; if (debug_infolevel>20-dim) CERR << "gcdmod content in " << "dim " << dim << " " << CLOCK() << endl; // Make p and q primitive as polynomial in xn with coeff in x1...xn-1 if (!pp_mod(p,&pmin,modulo,vars,pcont,nthreads)) return 0; if (!pp_mod(q,&pmin,modulo,vars,qcont,nthreads)) return 0; mod_gcd_ext(pcont,qcont,vars_truncated,pmin,modulo,dcont,pcofactor,qcofactor,false,false,nthreads); // don't use pv and qv here! // multiply pcofactor and qcofactor by the initial contents dep. on xn if (debug_infolevel>20-dim) CERR << "gcdmod content end " << "dim " << dim << " " << CLOCK() << endl; distmult_ext(dcont,dcontxn,dcont,var2,pmin,modulo); // ready for gcd computation by interpolation with respect to xn // first find degree of gcd with respect to xn int pxndeg=degree_xn(p,shiftxn,shift2),qxndeg=degree_xn(q,shiftxn,shift2),gcddeg=0; vector< vector > pb(pxndeg+1),qb(qxndeg+1),db; vector b(dim-1),bnext(dim-1); index_t vzero; // coeff of vzero correspond to zero or non zero int nzero=1; // Number of zero coeffs for (int essai=0;essai<2;){ if (debug_infolevel>20-dim) CERR << "gcdmod degree? " << essai << " dim " << dim << " " << CLOCK() << endl; if (!horner(p,b,vars,pb,modulo) || !horner(q,b,vars,qb,modulo)) return false; for (;;){ for (int i=0;i tmp(invmod(pcofactor.front().g,pminmodulo)); mulext(tmp,p_orig.front().g,pmin,modulo); smallmult(tmp,pcofactor,pcofactor,pminmodulo); } if (compute_qcofactor){ smallmult(q,qcofactor,qcofactor,pminmodulo,0); vector tmp(invmod(qcofactor.front().g,pminmodulo)); mulext(tmp,q_orig.front().g,pmin,modulo); smallmult(tmp,qcofactor,qcofactor,pminmodulo); } return true; } if (!essai){ // 1st gcd test gcddeg=dbdeg; nzero=find_nonzero(db,vzero); ++essai; continue; } // 2nd try if (dbdeg > lcoeffp,lcoeffq,lcoeffg,Delta,tmpcont; hashgcd_U lcoeffpu,lcoeffqu; if (debug_infolevel>20-dim) CERR << "gcdmod lcoeff begin " << "dim " << dim << " " << CLOCK() << endl; lcoeffpu=lcoeff(p,varxn,var2,lcoeffp); lcoeffqu=lcoeff(q,varxn,var2,lcoeffq); if (!gcdsmallmodpoly_ext(lcoeffp,lcoeffq,pmin,modulo,Delta)) return 0; if (debug_infolevel>20-dim){ CERR << "lcoeff p, q, gcd" << lcoeffp << "," << lcoeffq << "," << Delta << endl; CERR << "gcdmod lcoeff end " << "dim " << dim << " " << CLOCK() << endl; } // estimate time for full lift or division try // size=p.size()+q.size() // sumdeg=pxndeg+qxndeg // %age=gcddeg/min(pxndeg,qxndeg) // %age^dim*(1-%age)^dim*size^2 estimates the time for division try // gcddeg*size estimates the time for lifting to gcddeg // sumdeg*size estimates the time for full lifting // if sumdeg<(gcddeg+%age^dim*(1-%age)^dim*size) do full lifting int Deltadeg = int(Delta.size())-1,liftdeg=(compute_qcofactor?qxndeg:pxndeg)+Deltadeg; int gcddeg_plus_delta=gcddeg+Deltadeg; int liftdeg0=giacmax(liftdeg-gcddeg,gcddeg_plus_delta); // once liftdeg0 is reached we can replace g/gp/gq computation // by a check that d*dp=dxn*lcoeff(d*dp)/Delta at alpha // and d*dq=dxn*lcoeff(d*dq)/lcoeff(qxn) at alpha int sumdeg = pxndeg+qxndeg; double percentage = double(gcddeg)/giacmin(pxndeg,qxndeg); int sumsize = int(p.size()+q.size()); // ? add a malus factor for division double gcdlift=gcddeg+std::pow(percentage,dim)*std::pow(1-percentage,dim)*sumsize; bool compute_cof = false; // dim==2 || sumdeg20-dim) CERR << "gcdmod ext degree begin " << "dim " << dim << " " << CLOCK() << " compute_cof " << compute_cof << "(" << sumdeg/gcdlift << ")" << endl; int ptotaldeg=degree(p,shift_vars,pdeg); int qtotaldeg=degree(q,shift_vars,qdeg); if (debug_infolevel>20-dim){ CERR << "pdeg " << pdeg << " " << ptotaldeg << endl; CERR << "qdeg " << qdeg << " " << qtotaldeg << endl; } if (debug_infolevel>20-dim) CERR << "gcdmod degree end " << "dim " << dim << " " << CLOCK() << endl; int spdeg=0,sqdeg=0; for (int i=0;i20-dim) CERR << "gcdmod find alpha dim " << dim << " " << CLOCK() << endl; if (debug_infolevel>20-dim) CERR << " p " << p << " q " << q << endl; vector< T_unsigned,hashgcd_U> > palpha,qalpha,dp,dq; // d, dp and dq are the current interpolated values of gcd and cofactors vector alphav; vector< vector< T_unsigned,hashgcd_U> > > gcdv,pcofactorv,qcofactorv; gcd_call_param< vector > gcd_par; gcd_par.Delta=Δ gcd_par.lcoeffp=&lcoeffp; gcd_par.lcoeffq=&lcoeffq; gcd_par.alphav=&alphav; gcd_par.pv=0; gcd_par.qv=0; gcd_par.dim2gcdv=0; gcd_par.dim2pcofactorv=0; gcd_par.dim2qcofactorv=0; gcd_par.p=&p; gcd_par.q=&q; gcd_par.gcdv=&gcdv; gcd_par.pcofactorv=&pcofactorv; gcd_par.qcofactorv=&qcofactorv; gcd_par.pdeg=&pdeg; gcd_par.qdeg=&qdeg; gcd_par.vars=&vars; gcd_par.vars_truncated=&vars_truncated; gcd_par.shift_vars=&shift_vars; gcd_par.shift_vars_truncated=&shift_vars_truncated; gcd_par.compute_cof=compute_cof; gcd_par.compute_qcofactor=compute_qcofactor; gcd_par.dim2=false; gcd_par.pminptr=&pmin; gcd_par.modulo=modulo; // Warning: leaving nthreads > 1 is a bad idea if too many allocations // with lock happen if (dim>3 && sumsize > modgcd_cachesize ){ gcd_par.nthreads=nthreads; nthreads=1; } else gcd_par.nthreads=1; if (nthreads>gcddeg_plus_delta) nthreads=gcddeg_plus_delta+1; if (nthreads>1){ int todo=compute_cof?(liftdeg0+1):(gcddeg_plus_delta+1); double nth=todo/double(nthreads); // if (nthreads>=4 && (todo%nthreads==0)) nth = (todo+1)/double(nthreads); // keep one proc for a bad prime nth=std::ceil(nth); nthreads=int(std::ceil(todo/nth)); if (debug_infolevel>20-dim) CERR << "Using " << nthreads << " threads " << nth << " " << todo/nth << endl; } // return -1; // avoid inf loop while not implemented /* ******************************* BEGIN LOOP ******************************* */ int gcd_ext_ok=1; for (alpha=-1;;){ if (gcd_ext_ok!=1) return -1; // First check if we are ready to interpolate if (!compute_cof && e>gcddeg_plus_delta){ if (debug_infolevel>20-dim){ CERR << "gcdmod before interp " << dim << " clock= " << CLOCK() << gcdv << endl; } interpolate(alphav,gcdv,d,var2,modulo); vector< T_unsigned,hashgcd_U> > pquo,qquo,tmprem,pD(d); pp_mod(pD,&pmin,modulo,varxn,var2,tmpcont); if (debug_infolevel>20-dim){ CERR << "gcdmod pp1mod dim " << dim << " clock= " << CLOCK() << " d " << d << endl; CERR << "gcdmod alphav " << alphav << endl << "gcdv " << gcdv << endl << "gcdmod content " << tmpcont << endl; } // This removes the polynomial in xn that we multiplied by // (it was necessary to know the lcoeff of the interpolated poly) if (debug_infolevel>20-dim) CERR << "gcdmod check dim " << dim << " " << CLOCK() << endl; // Now, gcd divides pD for gcddeg+1 values of x1 // degree(pD)<=degree(gcd) // ?CHECK? should we pseudo-divide? if (hashdivrem(p,pD,pquo,tmprem,vars,pminmodulo,0,true)==1 && tmprem.empty()){ if (hashdivrem(q,pD,qquo,tmprem,vars,pminmodulo,0,true)==1 && tmprem.empty()){ smallmult(pD,dcont,d,pminmodulo,0); smallmult(invmod(d.front().g,pminmodulo),d,d,pminmodulo); if (compute_pcofactor){ smallmult(pcofactor,pquo,pcofactor,pminmodulo,0); vector tmp(invmod(pcofactor.front().g,pminmodulo)); mulext(tmp,p_orig.front().g,pmin,modulo); smallmult(tmp,pcofactor,pcofactor,pminmodulo); } if (compute_qcofactor){ smallmult(qcofactor,qquo,qcofactor,pminmodulo,0); vector tmp(invmod(qcofactor.front().g,pminmodulo)); mulext(tmp,q_orig.front().g,pmin,modulo); smallmult(tmp,qcofactor,qcofactor,pminmodulo); } if (debug_infolevel>20-dim) CERR << "gcdmod found dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); return true; } // end if hashdivrem(q,...) } // end if hashdivrem(p,...) if (debug_infolevel>20-dim) CERR << "Gcdmod bad guess " << endl; // restart gcdv.clear(); alphav.clear(); pcofactorv.clear(); qcofactorv.clear(); e=0; } // end if (e>gcddeg+delta) if (compute_cof && e>liftdeg0 ){ // interpolate d and dp interpolate(alphav,gcdv,d,var2,modulo); if (debug_infolevel>20-dim) CERR << "end interpolate gcd " << CLOCK() << endl; interpolate(alphav,pcofactorv,dp,var2,modulo); if (debug_infolevel>20-dim) CERR << "end interpolate p cof " << CLOCK() << endl; // check that d(alpha)*dp(alpha)=palpha with lcoeff adjusted // for e<=liftdeg vector< T_unsigned,hashgcd_U> > g,gp; for (++alpha;e<=liftdeg;++e,++alpha){ alpha1=alpha%2?-(alpha+1)/2:alpha/2; while (is_zero(hornermod_ext(lcoeffp,alpha1,modulo))){ ++alpha; alpha1=alpha%2?-(alpha+1)/2:alpha/2; } #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || alpha>=modulo) return false; int maxtotaldeg=ptotaldeg+1-e; if (!horner(p,alpha1,vars,palpha,modulo,maxtotaldeg)) return false; if (debug_infolevel>20-dim) CERR << "gcdmod horner d " << alpha << " dim " << dim << " " << CLOCK() << endl; if (!horner(d,alpha1,vars,g,modulo,maxtotaldeg)) return false; if (debug_infolevel>20-dim) CERR << "gcdmod horner dp " << alpha << " dim " << dim << " " << CLOCK() << endl; if (!horner(dp,alpha1,vars,gp,modulo,maxtotaldeg)) return false; vector tmp; mulext(gp.back().g,g.back().g,pmin,modulo,tmp); invmod(tmp,pminmodulo); mulext(tmp,palpha.back().g,pmin,modulo); smallmult(tmp,gp,gp,pminmodulo); // FIXME is_p_a_times_b if (true // !is_p_a_times_b(palpha,gp,g,vars,modulo,maxtotaldeg) ){ // Bad guess, go find some new gcd and interpolate e=liftdeg0+1; break; } } // end for ( e loop ) if (e>liftdeg){ // enough evaluation point // divide d,dp,dq by their content in xn pp_mod(d,&pmin,modulo,varxn,var2,tmpcont); pp_mod(dp,&pmin,modulo,varxn,var2,tmpcont); // check xn degrees of d+dp=degree(pxn), d+dq=degree(qxn) int dxndeg=degree_xn(d,shiftxn,shift2),dpxndeg=degree_xn(dp,shiftxn,shift2); // int dqxndeg=degree_xn(dq,shiftxn,shift2); if ( dxndeg+dpxndeg==pdeg.back() ){ smallmult(d,dcont,d,pminmodulo,0); if (compute_pcofactor){ smallmult(dp,pcofactor,pcofactor,pminmodulo,0); vector tmp(invmod(pcofactor.front().g,pminmodulo)); mulext(tmp,p_orig.front().g,pmin,modulo); smallmult(tmp,pcofactor,pcofactor,pminmodulo); } if (compute_qcofactor){ interpolate(alphav,qcofactorv,dq,var2,modulo); pp_mod(dq,&pmin,modulo,varxn,var2,tmpcont); smallmult(dq,qcofactor,qcofactor,pminmodulo,0); vector tmp(invmod(qcofactor.front().g,pminmodulo)); mulext(tmp,q_orig.front().g,pmin,modulo); smallmult(tmp,qcofactor,qcofactor,pminmodulo); } if (debug_infolevel>20-dim) CERR << "gcdmod end dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); // divtest=false; return true; } // failure, restart gcdv.clear(); alphav.clear(); pcofactorv.clear(); qcofactorv.clear(); e=0; } // end if (e>lifdeg) } // end if (compute_cof && e in liftdeg0..liftdeg) // *************************************************** // // Not ready to interpolate, try to eval new points // *************************************************** // for (int thread=0;thread20-dim) CERR << "gcdmod eval alpha1=" << alpha1 << " dim " << dim << " " << CLOCK() << endl; break; } // end for (;;) // alpha is probably admissible // (we will test later if degree of palpha/qalpha wrt x1...xn-1 is max) // prepare room for gcd and cofactors if (debug_infolevel>20-dim) CERR << "dim " << dim << " palpha " << palpha << " qalpha " << qalpha << endl ; alphav.push_back(alpha1); gcdv.push_back(vector< T_unsigned< vector,hashgcd_U> >(0)); pcofactorv.push_back(vector< T_unsigned< vector,hashgcd_U> >(0)); qcofactorv.push_back(vector< T_unsigned< vector,hashgcd_U> >(0)); } // end for (int thread=0;thread > > gcd_call_param_v(nthreads,gcd_par); #ifdef HAVE_PTHREAD_H pthread_t tab[nthreads-1]; #endif for (int thread=0;thread > *) ptr)->ext_gcd_ok; } else { if (gcd_ext_ok==1) gcd_ext_ok=gcd_call_param_v[thread].ext_gcd_ok; } #else if (gcd_ext_ok==1) gcd_ext_ok=gcd_call_param_v[thread].ext_gcd_ok; #endif if (gcd_ext_ok!=1) continue; // improve: kill other threads } if (gcd_ext_ok!=1) continue; for (int thread=0;thread2 && gcddeg-nzero==e){ // We have enough evaluations, let's try SPMOD // Build the matrix, each line has coeffs / vzero vector< vector > m,minverse; for (int j=0;j<=e;++j){ index_t::const_reverse_iterator it=vzero.rbegin(),itend=vzero.rend(); vector line; for (int p=alphav[j],pp=1;it!=itend;++it,pp=smod(p*pp,modulo)){ if (*it) line.push_back(pp); } reverse(line.begin(),line.end()); m.push_back(line); } // assume gcd is the vector of non zero coeffs of the gcd in x^n // we have the relation // m*gcd=gcdv // invert m (if invertible) longlong det_mod_p; if (smallmodinv(m,minverse,modulo,det_mod_p) && det_mod_p){ // hence gcd=minverse*gcdv, where the i-th component of gcd // must be "multiplied" by xn^degree_corresponding_vzero[i] vector< vector< T_unsigned > > minversegcd(e+1); for (int j=0;j<=e;++j){ for (int k=0;k<=e;++k){ vector< T_unsigned > tmp(gcdv[k]); smallmult(minverse[j][k],tmp,tmp,modulo); smalladd(minversegcd[j],tmp,minversegcd[j],modulo); // CERR << minversegcd[j] << endl; } } vector< T_unsigned > trygcd,pquo,qquo,tmprem; index_t::const_iterator it=vzero.begin(),itend=vzero.end(); int deg=itend-it-1; for (int j=0;it!=itend;++it,--deg){ if (!*it) continue; smallshift(minversegcd[j],deg*var2,minversegcd[j]); smalladd(trygcd,minversegcd[j],trygcd,modulo); ++j; } // Check if trygcd is the gcd! pp_mod(trygcd,0,modulo,varxn,var2,tmpcont); if (hashdivrem(p,trygcd,pquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ if (hashdivrem(q,trygcd,qquo,tmprem,vars,modulo,0,false)==1 && tmprem.empty()){ smallmult(trygcd,dcont,d,modulo,0); smallmult(invmod(d.front().g,modulo),d,d,modulo); if (compute_pcofactor){ smallmult(pcofactor,pquo,pcofactor,modulo,0); smallmult(smod(longlong(p_orig.front().g)*invmod(pcofactor.front().g,modulo),modulo),pcofactor,pcofactor); } if (compute_qcofactor){ smallmult(qcofactor,qquo,qcofactor,modulo,0); smallmult(smod(longlong(q_orig.front().g)*invmod(qcofactor.front().g,modulo),modulo),qcofactor,qcofactor); } if (debug_infolevel>20-dim) CERR << "gcdmod found dim " << dim << " " << CLOCK() << endl; if (pqswap) swap(pcofactor,qcofactor); return true; } // end q divisible by trygcd } // end p divisible by trygcd } // end m invertible } // end if (dim>2 && ) } // end SPMOD */ if (debug_infolevel>20-dim) CERR << "gcdmod interp dim " << dim << " " << CLOCK() << endl; ++e; continue; } // end gdeg==delta if (comp==0 || comp==-1){ if (debug_infolevel>20-dim) CERR << "Bad reduction " << alphav[vpos] << endl; // bad reduction: all indices of gdeg are >= to delta and gdeg!=delta alphav.erase(alphav.begin()+vpos); gcdv.erase(gcdv.begin()+vpos); pcofactorv.erase(pcofactorv.begin()+vpos); qcofactorv.erase(qcofactorv.begin()+vpos); if (comp==0) continue; } // previous alpha where bad reduction if (debug_infolevel>20-dim && vpos) CERR << "Bads reductions " << alphav[vpos-1] << endl; alphav.erase(alphav.begin(),alphav.begin()+vpos); gcdv.erase(gcdv.begin(),gcdv.begin()+vpos); pcofactorv.erase(pcofactorv.begin(),pcofactorv.begin()+vpos); qcofactorv.erase(qcofactorv.begin(),qcofactorv.begin()+vpos); if (comp==-1){ e=0; continue; } // restart everything with this value of alpha // this will happen (almost all the time) at first iteration delta=gdeg; e=1; continue; } // end for (int thread=0;thread > & p_orig,const gen & modulo,vector< T_unsigned > & p,const gen & pimod){ vector< T_unsigned > q; vector< T_unsigned< vecteur,hashgcd_U> >::const_iterator it=p_orig.begin(),itend=p_orig.end(); vector< T_unsigned >::const_iterator jt=p.begin(),jtend=p.end(); q.reserve(jtend-jt); for (;it!=itend && jt!=jtend;){ if (it->u==jt->u){ const vecteur & tmp=it->g; q.push_back(T_unsigned(ichinrem_ext(tmp,jt->g,modulo,pimod),it->u)); ++it; ++jt; continue; } if (it->uu){ q.push_back(T_unsigned(ichinrem_ext(vecteur(0),jt->g,modulo,pimod),jt->u)); ++jt; } else { const vecteur & tmp=it->g; q.push_back(T_unsigned(ichinrem_ext(tmp,0,modulo,pimod),it->u)); ++it; } } for (;it!=itend;++it){ const vecteur & tmp=it->g; q.push_back(T_unsigned(ichinrem_ext(tmp,0,modulo,pimod),it->u)); } for (;jt!=jtend;++jt) q.push_back(T_unsigned(ichinrem_ext(vecteur(0),jt->g,modulo,pimod),jt->u)); swap(p,q); } // add pmin extension to all vecteurs inside p static void make_ext(vector< T_unsigned > &p,const gen & pmin){ vector< T_unsigned >::iterator jt=p.begin(),jtend=p.end(); for (;jt!=jtend;++jt){ if (jt->g.type==_VECT) jt->g=algebraic_EXTension(jt->g,pmin); } } static bool fracmod(const vector< T_unsigned > &p,const gen & modulo,vector< T_unsigned > & q){ q=p; gen tmp; vector< T_unsigned >::iterator jt=q.begin(),jtend=q.end(); for (;jt!=jtend;++jt){ if (!fracmod(jt->g,modulo,tmp)) return false; jt->g=tmp; } return true; } // should be called from gausspol.cc in gcdheu after monomial packing like in modpoly.cc gcd_modular bool gcd_ext(const vector< T_unsigned > & p_orig,const vector< T_unsigned > & q_orig,vector< T_unsigned > & d, vector< T_unsigned > & pcofactor, vector< T_unsigned > & qcofactor,const std::vector & vars, bool compute_cofactors,int nthreads){ #ifdef NO_TEMPLATE_MULTGCD return false; #else // find extension compute_cofactors=false; gen coefft; int tp=is_integer(p_orig,coefft),tq=is_integer(q_orig,coefft); if (!tp || !tq) return false; if (tp<3 && tq<3) return false; if (tp==3 && tq==2) tp=4; if (tq==3 && tp==2) tp=4; if (tq>tp) tp=tq; index_t shift_vars; if (!find_shift(vars,shift_vars)) return false; coefft=*(coefft._EXTptr+1); if (coefft.type!=_VECT || coefft._VECTptr->empty() || !is_integer_vecteur(*coefft._VECTptr)) return false; vecteur & pminv=*coefft._VECTptr; gen pmin0=pminv.front(); gen m=536871000,pimod=1; gen lcoeffp=p_orig.front().g,lcoeffq=q_orig.front().g; d.clear(); pcofactor.clear(); qcofactor.clear(); vector< T_unsigned > p,q,g,pcof,qcof; index_t pdeg,qdeg,pdegmod,qdegmod,gdegmod,gdegmod2,gdeg; degree(p_orig,shift_vars,pdeg); degree(q_orig,shift_vars,qdeg); gdeg=min(pdeg,qdeg); if (tp==4){ // return false; vector< T_unsigned > p1,q1,g1,pcof1,qcof1,p2,q2,g2,pcof2,qcof2; // tp==4, use primes == 1 % 4, find gcds replacing i by both // sqrts of -1 modulo the primes, ichinrem for (int n=0;;n++){ m=nextprime(m+1); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || m.type!=_INT_) return false; int modulo=m.val; // computing gcd in Z[i]: use a modulo =1[4], find gcd for both roots of -1 mod modulo if (modulo%4==3) continue; // min poly mod modulo, coeffs should not contain i==sqrt(-1) if (is_zero(smod(pmin0,modulo))) continue; if (debug_infolevel>(int)shift_vars.size()) CERR << "Modular algebraic extension gcd " << modulo << endl; vecteur pmin; pmin.reserve(pminv.size()); const_iterateur it=pminv.begin(),itend=pminv.end(); for (;it!=itend;++it){ pmin.push_back(smod(*it,modulo)); } // find square root of -1 int i=modsqrtminus1(modulo); gen lp1=complex_smod_ext(lcoeffp,i,modulo),lp2=complex_smod_ext(lcoeffp,-i,modulo); gen lq1=complex_smod_ext(lcoeffq,i,modulo),lq2=complex_smod_ext(lcoeffq,-i,modulo); if (is_zero(lp1) || is_zero(lq1) || is_zero(lp2) || is_zero(lq2)) continue; if (!complex_smod_ext(p_orig,i,modulo,p1)) return false; // setsizeerr(); degree(p1,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; if (!complex_smod_ext(q_orig,i,modulo,q1)) return false; // setsizeerr(); degree(q1,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; if (!complex_smod_ext(p_orig,-i,modulo,p2)) return false; // setsizeerr(); degree(p2,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; if (!complex_smod_ext(q_orig,-i,modulo,q2)) return false; // setsizeerr(); degree(q2,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; int res=gcd_ext(p1,q1,vars,pmin,modulo,g1,pcof1,qcof1,compute_cofactors,compute_cofactors,nthreads); // ?FIXME? check that g1, pcof1, qcof1 are normalized #ifdef TIMEOUT control_c(); #endif if (res==-1 || ctrl_c || interrupted) return false; if (!res) continue; res=gcd_ext(p2,q2,vars,pmin,modulo,g2,pcof2,qcof2,compute_cofactors,compute_cofactors,nthreads); #ifdef TIMEOUT control_c(); #endif if (res==-1 || ctrl_c || interrupted) return false; if (!res) continue; // CERR << p1 << " " << q1 << " " << g1 << endl << p2 << " " << q2 << " " << g2 << endl ; // IMPROVE? use a map for modulo -> g1,g2,etc. // and do chinese remaindering for poly having the same lcoeff degree(g1,shift_vars,gdegmod); degree(g2,shift_vars,gdegmod2); if (gdegmod2!=gdegmod) continue; if (gdegmod.front()>gdeg.front()) // unlucky prime continue; int cmp=compare(gdegmod,gdeg); if (cmp==-1){ // full restart d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; continue; } if (cmp!=-2){ // restart with this gcd gdeg=gdegmod; d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; } complex_ichinrem_ext(g1,g2,i,modulo,d,pimod); if (compute_cofactors){ complex_ichinrem_ext(pcof1,pcof2,i,modulo,pcofactor,pimod); complex_ichinrem_ext(qcof1,qcof2,i,modulo,qcofactor,pimod); } pimod = modulo * pimod; // rational reconstruction and division test vector< T_unsigned > dtest,pcofactortest,qcofactortest,pquo,qquo,rem; fracmod(d,pimod,dtest); lcmdeno(dtest); gen tmp; ppz(dtest,tmp,true); make_ext(dtest,coefft); // make extensions if (compute_cofactors){ fracmod(pcofactor,pimod,pcofactortest); fracmod(qcofactor,pimod,qcofactortest); lcmdeno(pcofactortest); ppz(pcofactortest,tmp,true); lcmdeno(qcofactortest); ppz(qcofactortest,tmp,true); // lcoeff(p)*dtest*pcofactortest =? p*lcoeff(dtest)*lcoeff(pcofactortest) // this is true mod pimod, check sizes // If AxB=C mod m,Pmin(theta) // then A*B-C=mD mod Pmin(theta) // The size of coeffs of A*B-C is <= // <= |A|*|B|* min(size(A),size(B)) * (degpmin+1) +|C| // but division by Pmin(theta) may enlarge the remainder // by a multiplication by ((degpmin+1)*|Pmin|)^(degpmin-1) // remove denominators and integer content gen lcoeffdpcof=dtest.front().g*pcofactortest.front().g; simplify(lcoeffp,lcoeffdpcof); gen lcoeffdqcof=dtest.front().g*qcofactortest.front().g; simplify(lcoeffq,lcoeffdqcof); gen maxp=max(p_orig,context0),maxq=max(q_orig,context0),maxpcof=max(pcofactortest,context0),maxqcof=max(qcofactortest,context0),maxd=max(dtest,context0); gen maxpmin=_max(coefft,context0),degpmin=int(pminv.size()-1); if (is_undef(maxpmin)) return false; gen multpmin=pow((degpmin+1)*maxpmin,degpmin,context0); gen test1=lcoeffdpcof*maxp+lcoeffp*maxd*maxpcof*gen(int(std::min(dtest.size(),pcofactortest.size())))*(degpmin+1)*multpmin; gen test2=lcoeffdqcof*maxq+lcoeffq*maxd*maxqcof*gen(int(std::min(dtest.size(),qcofactortest.size())))*(degpmin+1)*multpmin; if (is_strictly_greater(pimod,test1,context0) && is_strictly_greater(pimod,test2,context0)){ // leave d unchanged but adjust pcofactor and qcofactor swap(d,dtest); smallmult(lcoeffp/lcoeffdpcof,pcofactortest,pcofactor); smallmult(lcoeffq/lcoeffdqcof,qcofactortest,qcofactor); return true; } } if (hashdivrem(p_orig,dtest,pquo,rem,vars,0 /* reduce */,0/*qmax*/,true)==1 && rem.empty()){ if (hashdivrem(q_orig,dtest,qquo,rem,vars,0 /* reduce */,0/*qmax*/,true)==1 && rem.empty()){ // FIXME normalize d swap(pcofactor,pquo); swap(qcofactor,qquo); swap(d,dtest); return true; } } } } // begin loop on modulo vector< T_unsigned > dtestold; for (int n=0;;n++){ m=nextprime(m+1); #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || m.type!=_INT_) return false; int modulo=m.val; if (debug_infolevel>(int)shift_vars.size()) CERR << "Modular algebraic extension gcd " << modulo << endl; // min poly mod modulo if (is_zero(smod(pmin0,modulo))) continue; vecteur pmin; pmin.reserve(pminv.size()); const_iterateur it=pminv.begin(),itend=pminv.end(); for (;it!=itend;++it){ pmin.push_back(smod(*it,modulo)); } // reduce mod m if (!smod_ext(p_orig,m.val,p) || !smod_ext(q_orig,m.val,q)) return false; // check degrees of p and q mod modulo degree(p,shift_vars,pdegmod); if (pdegmod!=pdeg) continue; degree(q,shift_vars,qdegmod); if (qdegmod!=qdeg) continue; // if one lcoeff of p/q is not invertible mod modulo // then it won't be when variables of p or q are evaluated // hence checking in dim 1 is sufficient if (debug_infolevel>0) CERR << CLOCK()*1e-6 << " begin gcd_ext mod " << modulo <<"," << pmin << endl; int res=gcd_ext(p,q,vars,pmin,modulo,g,pcof,qcof,compute_cofactors,compute_cofactors,nthreads); if (debug_infolevel>0) CERR << CLOCK()*1e-6 << " end gcd_ext mod " << modulo <<"," << pmin << endl; #ifdef TIMEOUT control_c(); #endif if (ctrl_c || interrupted || res==-1) return false; if (!res) continue; // ?FIXME? check that g is normalized to have lcoeff==1 // IMPROVE? use a map for modulo -> lcoeff -> g, pcof, qcof, pimodulo // and chinese remainder with modulos having the same lcoeff degree(g,shift_vars,gdegmod); if (gdegmod.front()>gdeg.front()) // unlucky prime continue; // Check with respect to other variables // N.B.: it could be faster to store all g,pcof,qcof // to avoid a full restart each time an unlucky content is reached // Here we must have sufficiently successive primes // that do not have unlucky content int cmp=compare(gdegmod,gdeg); if (cmp==-1){ // full restart d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; continue; } if (cmp!=-2){ // restart with this gcd gdeg=gdegmod; d.clear(); pcofactor.clear(); qcofactor.clear(); pimod=1; } // Same degrees // chinese remainder of g and d ichinrem_ext(g,modulo,d,pimod); if (compute_cofactors){ ichinrem_ext(pcof,modulo,pcofactor,pimod); ichinrem_ext(qcof,modulo,qcofactor,pimod); CERR << "g " << d << endl; CERR << "pcof " << pcofactor << endl; CERR << "qcof " << qcofactor << endl; } pimod = modulo * pimod; // rational reconstruction and division test vector< T_unsigned > dtest,pcofactortest,qcofactortest,pquo,qquo,rem; fracmod(d,pimod,dtest); lcmdeno(dtest); gen tmp; ppz(dtest,tmp,true); make_ext(dtest,coefft); // make extensions if (0 && compute_cofactors){ // DISABLED: I don't know how equality mod m could be translated // (absolute values of algebraic extension) fracmod(pcofactor,pimod,pcofactortest); fracmod(qcofactor,pimod,qcofactortest); lcmdeno(pcofactortest); ppz(pcofactortest,tmp,true); lcmdeno(qcofactortest); ppz(qcofactortest,tmp,true); // lcoeff(p)*dtest*pcofactortest =? p*lcoeff(dtest)*lcoeff(pcofactortest) // this is true mod pimod, check sizes // If AxB=C mod m,Pmin(theta) // then A*B-C=mD mod Pmin(theta) // The size of coeffs of A*B-C is <= // <= |A|*|B|* min(size(A),size(B)) * (degpmin+1) +|C| // but division by Pmin(theta) may enlarge the remainder // by a multiplication by ((degpmin+1)*|Pmin|)^(degpmin-1) // remove denominators and integer content gen lcoeffdpcof=dtest.front().g*pcofactortest.front().g; simplify(lcoeffp,lcoeffdpcof); gen lcoeffdqcof=dtest.front().g*qcofactortest.front().g; simplify(lcoeffq,lcoeffdqcof); gen maxp=max(p_orig,context0); gen maxq=max(q_orig,context0); gen maxpcof=max(pcofactortest,context0); gen maxqcof=max(qcofactortest,context0); gen maxd=max(dtest,context0); gen maxpmin=_max(coefft,context0),degpmin=int(pminv.size()-1); if (is_undef(maxpmin)) return false; gen multpmin=pow((degpmin+1)*maxpmin,degpmin,context0); gen test1=lcoeffdpcof*maxp+lcoeffp*maxd*maxpcof*gen(int(std::min(dtest.size(),pcofactortest.size())))*(degpmin+1)*multpmin; gen test2=lcoeffdqcof*maxq+lcoeffq*maxd*maxqcof*gen(int(std::min(dtest.size(),qcofactortest.size())))*(degpmin+1)*multpmin; if (is_strictly_greater(pimod,test1,context0) && is_strictly_greater(pimod,test2,context0)){ // leave d unchanged but adjust pcofactor and qcofactor swap(d,dtest); smallmult(lcoeffp/lcoeffdpcof,pcofactortest,pcofactor); smallmult(lcoeffq/lcoeffdqcof,qcofactortest,qcofactor); return true; } } // This division might take a very long time if not successfull // Make it only when we are sure (might be improved) if (dtest==dtestold){ if (debug_infolevel>0) CERR << CLOCK()*1e-6 << " gcd_ext checking gcd" << endl; if (hashdivrem(p_orig,dtest,pquo,rem,vars,0 /* reduce */,0/*qmax*/,true)==1 && rem.empty()){ // CERR << "p " << p_orig << endl << "dtest " << dtest << endl << "quo" << pquo << endl; if (debug_infolevel>0) CERR << CLOCK()*1e-6 << " gcd_ext checking gcd 1st test ok" << endl; if (hashdivrem(q_orig,dtest,qquo,rem,vars,0 /* reduce */,0/*qmax*/,true)==1 && rem.empty()){ // CERR << "q " << q_orig << endl << "dtest " << dtest << endl << "quo" << qquo << endl; if (debug_infolevel>0) CERR << CLOCK()*1e-6 << " gcd_ext checking gcd 2nd test ok" << endl; swap(pcofactor,pquo); swap(qcofactor,qquo); swap(d,dtest); return true; } } } dtestold=dtest; } #endif // NO_TEMPLATE_MULTGCD } #ifndef NO_NAMESPACE_GIAC } // namespace giac #endif // ndef NO_NAMESPACE_GIAC