?? zz_px1.c
字號:
void ProbMinPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m) { long n = F.n; if (m < 1 || m > n) Error("ProbMinPoly: bad args"); long i; vec_ZZ_p R(INIT_SIZE, n); for (i = 0; i < n; i++) random(R[i]); DoMinPolyMod(h, g, F, m, R);}void MinPolyMod(ZZ_pX& hh, const ZZ_pX& g, const ZZ_pXModulus& F, long m) { ZZ_pX h, h1; long n = F.n; if (m < 1 || m > n) Error("MinPoly: bad args"); /* probabilistically compute min-poly */ ProbMinPolyMod(h, g, F, m); if (deg(h) == m) { hh = h; return; } CompMod(h1, h, g, F); if (IsZero(h1)) { hh = h; return; } /* not completely successful...must iterate */ long i; ZZ_pX h2, h3; ZZ_pXMultiplier H1; vec_ZZ_p R(INIT_SIZE, n); for (;;) { R.SetLength(n); for (i = 0; i < n; i++) random(R[i]); build(H1, h1, F); UpdateMap(R, R, H1, F); DoMinPolyMod(h2, g, F, m-deg(h), R); mul(h, h, h2); if (deg(h) == m) { hh = h; return; } CompMod(h3, h2, g, F); MulMod(h1, h3, H1, F); if (IsZero(h1)) { hh = h; return; } }}void IrredPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m) { vec_ZZ_p R(INIT_SIZE, 1); if (m < 1 || m > F.n) Error("IrredPoly: bad args"); set(R[0]); DoMinPolyMod(h, g, F, m, R);}void diff(ZZ_pX& x, const ZZ_pX& a){ long n = deg(a); long i; if (n <= 0) { clear(x); return; } if (&x != &a) x.rep.SetLength(n); for (i = 0; i <= n-1; i++) { mul(x.rep[i], a.rep[i+1], i+1); } if (&x == &a) x.rep.SetLength(n); x.normalize();}void MakeMonic(ZZ_pX& x){ if (IsZero(x)) return; if (IsOne(LeadCoeff(x))) return; ZZ_p t; inv(t, LeadCoeff(x)); mul(x, x, t);} void PlainMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n){ ZZ_pX y; mul(y, a, b); trunc(x, y, n);}void FFTMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n){ if (IsZero(a) || IsZero(b)) { clear(x); return; } long d = deg(a) + deg(b); if (n > d + 1) n = d + 1; long k = NextPowerOfTwo(d + 1); FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k); ToFFTRep(R1, a, k); ToFFTRep(R2, b, k); mul(R1, R1, R2); FromFFTRep(x, R1, 0, n-1);}void MulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n){ if (n < 0) Error("MulTrunc: bad args"); if (deg(a) <= NTL_ZZ_pX_FFT_CROSSOVER || deg(b) <= NTL_ZZ_pX_FFT_CROSSOVER) PlainMulTrunc(x, a, b, n); else FFTMulTrunc(x, a, b, n);} void PlainSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n){ ZZ_pX y; sqr(y, a); trunc(x, y, n);}void FFTSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n){ if (IsZero(a)) { clear(x); return; } long d = 2*deg(a); if (n > d + 1) n = d + 1; long k = NextPowerOfTwo(d + 1); FFTRep R1(INIT_SIZE, k); ToFFTRep(R1, a, k); mul(R1, R1, R1); FromFFTRep(x, R1, 0, n-1);}void SqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n){ if (n < 0) Error("SqrTrunc: bad args"); if (deg(a) <= NTL_ZZ_pX_FFT_CROSSOVER) PlainSqrTrunc(x, a, n); else FFTSqrTrunc(x, a, n);}void FastTraceVec(vec_ZZ_p& S, const ZZ_pX& f){ long n = deg(f); if (n <= 0) Error("FastTraceVec: bad args"); if (n == 0) { S.SetLength(0); return; } if (n == 1) { S.SetLength(1); set(S[0]); return; } long i; ZZ_pX f1; f1.rep.SetLength(n-1); for (i = 0; i <= n-2; i++) f1.rep[i] = f.rep[n-i]; f1.normalize(); ZZ_pX f2; f2.rep.SetLength(n-1); for (i = 0; i <= n-2; i++) mul(f2.rep[i], f.rep[n-1-i], i+1); f2.normalize(); ZZ_pX f3; InvTrunc(f3, f1, n-1); MulTrunc(f3, f3, f2, n-1); S.SetLength(n); S[0] = n; for (i = 1; i < n; i++) negate(S[i], coeff(f3, i-1));}void PlainTraceVec(vec_ZZ_p& S, const ZZ_pX& ff){ if (deg(ff) <= 0) Error("TraceVec: bad args"); ZZ_pX f; f = ff; MakeMonic(f); long n = deg(f); S.SetLength(n); if (n == 0) return; long k, i; ZZ acc, t; ZZ_p t1; S[0] = n; for (k = 1; k < n; k++) { mul(acc, rep(f.rep[n-k]), k); for (i = 1; i < k; i++) { mul(t, rep(f.rep[n-i]), rep(S[k-i])); add(acc, acc, t); } conv(t1, acc); negate(S[k], t1); }}void TraceVec(vec_ZZ_p& S, const ZZ_pX& f){ if (deg(f) <= NTL_ZZ_pX_TRACE_CROSSOVER) PlainTraceVec(S, f); else FastTraceVec(S, f);}void ComputeTraceVec(const ZZ_pXModulus& F){ vec_ZZ_p& S = *((vec_ZZ_p *) &F.tracevec); if (S.length() > 0) return; if (!F.UseFFT) { PlainTraceVec(S, F.f); return; } long i; long n = F.n; FFTRep R; ZZ_pX P, g; g.rep.SetLength(n-1); for (i = 1; i < n; i++) mul(g.rep[n-i-1], F.f.rep[n-i], i); g.normalize(); ToFFTRep(R, g, F.l); mul(R, R, F.HRep); FromFFTRep(P, R, n-2, 2*n-4); S.SetLength(n); S[0] = n; for (i = 1; i < n; i++) negate(S[i], coeff(P, n-1-i));}void TraceMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pXModulus& F){ long n = F.n; if (deg(a) >= n) Error("trace: bad args"); if (F.tracevec.length() == 0) ComputeTraceVec(F); InnerProduct(x, a.rep, F.tracevec);}void TraceMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pX& f){ if (deg(a) >= deg(f) || deg(f) <= 0) Error("trace: bad args"); project(x, TraceVec(f), a);}void PlainResultant(ZZ_p& rres, const ZZ_pX& a, const ZZ_pX& b){ ZZ_p res; if (IsZero(a) || IsZero(b)) clear(res); else if (deg(a) == 0 && deg(b) == 0) set(res); else { long d0, d1, d2; ZZ_p lc; set(res); long n = max(deg(a),deg(b)) + 1; ZZ_pX u(INIT_SIZE, n), v(INIT_SIZE, n); ZZVec tmp(n, ZZ_pInfo->ExtendedModulusSize); u = a; v = b; for (;;) { d0 = deg(u); d1 = deg(v); lc = LeadCoeff(v); PlainRem(u, u, v, tmp); swap(u, v); d2 = deg(v); if (d2 >= 0) { power(lc, lc, d0-d2); mul(res, res, lc); if (d0 & d1 & 1) negate(res, res); } else { if (d1 == 0) { power(lc, lc, d0); mul(res, res, lc); } else clear(res); break; } } rres = res; }}void ResIterHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red, vec_ZZ_p& cvec, vec_long& dvec){ M_out(0,0).SetMaxLength(d_red); M_out(0,1).SetMaxLength(d_red); M_out(1,0).SetMaxLength(d_red); M_out(1,1).SetMaxLength(d_red); set(M_out(0,0)); clear(M_out(0,1)); clear(M_out(1,0)); set(M_out(1,1)); long goal = deg(U) - d_red; if (deg(V) <= goal) return; ZZVec tmp(deg(U)+1, ZZ_pInfo->ExtendedModulusSize); ZZ_pX Q, t(INIT_SIZE, d_red); while (deg(V) > goal) { append(cvec, LeadCoeff(V)); append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V)); PlainDivRem(Q, U, U, V, tmp); swap(U, V); mul(t, Q, M_out(1,0)); sub(t, M_out(0,0), t); M_out(0,0) = M_out(1,0); M_out(1,0) = t; mul(t, Q, M_out(1,1)); sub(t, M_out(0,1), t); M_out(0,1) = M_out(1,1); M_out(1,1) = t; }} void ResHalfGCD(ZZ_pXMatrix& M_out, const ZZ_pX& U, const ZZ_pX& V, long d_red, vec_ZZ_p& cvec, vec_long& dvec){ if (IsZero(V) || deg(V) <= deg(U) - d_red) { set(M_out(0,0)); clear(M_out(0,1)); clear(M_out(1,0)); set(M_out(1,1)); return; } long n = deg(U) - 2*d_red + 2; if (n < 0) n = 0; ZZ_pX U1, V1; RightShift(U1, U, n); RightShift(V1, V, n); if (d_red <= NTL_ZZ_pX_HalfGCD_CROSSOVER) { ResIterHalfGCD(M_out, U1, V1, d_red, cvec, dvec); return; } long d1 = (d_red + 1)/2; if (d1 < 1) d1 = 1; if (d1 >= d_red) d1 = d_red - 1; ZZ_pXMatrix M1; ResHalfGCD(M1, U1, V1, d1, cvec, dvec); mul(U1, V1, M1); long d2 = deg(V1) - deg(U) + n + d_red; if (IsZero(V1) || d2 <= 0) { M_out = M1; return; } ZZ_pX Q; ZZ_pXMatrix M2; append(cvec, LeadCoeff(V1)); append(dvec, dvec[dvec.length()-1]-deg(U1)+deg(V1)); DivRem(Q, U1, U1, V1); swap(U1, V1); ResHalfGCD(M2, U1, V1, d2, cvec, dvec); ZZ_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1); mul(t, Q, M1(1,0)); sub(t, M1(0,0), t); swap(M1(0,0), M1(1,0)); swap(M1(1,0), t); t.kill(); t.SetMaxLength(deg(M1(1,1))+deg(Q)+1); mul(t, Q, M1(1,1)); sub(t, M1(0,1), t); swap(M1(0,1), M1(1,1)); swap(M1(1,1), t); t.kill(); mul(M_out, M2, M1); }void ResHalfGCD(ZZ_pX& U, ZZ_pX& V, vec_ZZ_p& cvec, vec_long& dvec){ long d_red = (deg(U)+1)/2; if (IsZero(V) || deg(V) <= deg(U) - d_red) { return; } long du = deg(U); long d1 = (d_red + 1)/2; if (d1 < 1) d1 = 1; if (d1 >= d_red) d1 = d_red - 1; ZZ_pXMatrix M1; ResHalfGCD(M1, U, V, d1, cvec, dvec); mul(U, V, M1); long d2 = deg(V) - du + d_red; if (IsZero(V) || d2 <= 0) { return; } M1(0,0).kill(); M1(0,1).kill(); M1(1,0).kill(); M1(1,1).kill(); ZZ_pX Q; append(cvec, LeadCoeff(V)); append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V)); DivRem(Q, U, U, V); swap(U, V); ResHalfGCD(M1, U, V, d2, cvec, dvec); mul(U, V, M1); }void resultant(ZZ_p& rres, const ZZ_pX& u, const ZZ_pX& v){ if (deg(u) <= NTL_ZZ_pX_GCD_CROSSOVER || deg(v) <= NTL_ZZ_pX_GCD_CROSSOVER) { PlainResultant(rres, u, v); return; } ZZ_pX u1, v1; u1 = u; v1 = v; ZZ_p res, t; set(res); if (deg(u1) == deg(v1)) { rem(u1, u1, v1); swap(u1, v1); if (IsZero(v1)) { clear(rres); return; } power(t, LeadCoeff(u1), deg(u1) - deg(v1)); mul(res, res, t); if (deg(u1) & 1) negate(res, res); } else if (deg(u1) < deg(v1)) { swap(u1, v1); if (deg(u1) & deg(v1) & 1) negate(res, res); } // deg(u1) > deg(v1) && v1 != 0 vec_ZZ_p cvec; vec_long dvec; cvec.SetMaxLength(deg(v1)+2); dvec.SetMaxLength(deg(v1)+2); append(cvec, LeadCoeff(u1)); append(dvec, deg(u1)); while (deg(u1) > NTL_ZZ_pX_GCD_CROSSOVER && !IsZero(v1)) { ResHalfGCD(u1, v1, cvec, dvec); if (!IsZero(v1)) { append(cvec, LeadCoeff(v1)); append(dvec, deg(v1)); rem(u1, u1, v1); swap(u1, v1); } } if (IsZero(v1) && deg(u1) > 0) { clear(rres); return; } long i, l; l = dvec.length(); if (deg(u1) == 0) { // we went all the way... for (i = 0; i <= l-3; i++) { power(t, cvec[i+1], dvec[i]-dvec[i+2]); mul(res, res, t); if (dvec[i] & dvec[i+1] & 1) negate(res, res); } power(t, cvec[l-1], dvec[l-2]); mul(res, res, t); } else { for (i = 0; i <= l-3; i++) { power(t, cvec[i+1], dvec[i]-dvec[i+2]); mul(res, res, t); if (dvec[i] & dvec[i+1] & 1) negate(res, res); } power(t, cvec[l-1], dvec[l-2]-deg(v1)); mul(res, res, t); if (dvec[l-2] & dvec[l-1] & 1) negate(res, res); PlainResultant(t, u1, v1); mul(res, res, t); } rres = res;}void NormMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pX& f){ if (deg(f) <= 0 || deg(a) >= deg(f)) Error("norm: bad args"); if (IsZero(a)) { clear(x); return; } ZZ_p t; resultant(t, f, a); if (!IsOne(LeadCoeff(f))) { ZZ_p t1; power(t1, LeadCoeff(f), deg(a)); inv(t1, t1); mul(t, t, t1); } x = t;}NTL_END_IMPL
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -