diff --git a/check/fixes.frm b/check/fixes.frm index 94c7198c..b64a609e 100644 --- a/check/fixes.frm +++ b/check/fixes.frm @@ -5613,3 +5613,129 @@ Print; assert succeeded? assert result("test") =~ expr("0") *--#] PullReq860_3 : +*--#[ PolyRatFunExpandModulus : +* Reduce expansion coefficients before sorting, including uncombined terms. +Modulus 13; +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,3); +Local F = rat(1,1+ep/2); +Local G = rat(1+ep,2+ep); +Local H = rat(1+ep,2); +Local N = rat(1,2); +Local S = rat(ep,2); +Local D = rat(1,2*ep); +Local L = rat(ep^2,2*ep+ep^2); +Local P2 = 6*rat(7); +* Continue beyond the quadratic coefficient that vanishes modulo 13. +Local R = rat(1,1+4*ep+3*ep^2); +Print +s; +.end +Modulus PlusMin,13; +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,3); +Local F = rat(1,1+ep/2); +Local G = rat(1+ep,2+ep); +Local H = rat(1+ep,2); +Local N = rat(1,2); +Local S = rat(ep,2); +Local D = rat(1,2*ep); +Local L = rat(ep^2,2*ep+ep^2); +Local P2 = 6*rat(7); +Local R = rat(1,1+4*ep+3*ep^2); +Print +s; +.end +assert succeeded? +modes = ["Modulus 13", "PlusMin 13"] +[ + ["F", "+ rat(1 + 6*ep + 10*ep^2 + 8*ep^3)", "+ rat(1 + 6*ep - 3*ep^2 - 5*ep^3)"], + ["G", "+ rat(7 + 10*ep + 8*ep^2 + 9*ep^3)", "+ rat(-6 - 3*ep - 5*ep^2 - 4*ep^3)"], + ["H", "+ rat(7 + 7*ep)", "+ rat(-6 - 6*ep)"], + ["N", "+ rat(7)", "+ rat(-6)"], + ["S", "+ rat(7*ep)", "+ rat(-6*ep)"], + ["D", "+ rat(7*ep^-1)", "+ rat(-6*ep^-1)"], + ["L", "+ rat(7*ep + 3*ep^2 + 5*ep^3 + 4*ep^4)", "+ rat(-6*ep + 3*ep^2 + 5*ep^3 + 4*ep^4)"], + ["P2", "+ rat(3)", "+ rat(3)"], + ["R", "+ rat(1 + 9*ep + 12*ep^3)", "+ rat(1 - 4*ep - ep^3)"] +].each do |name, *values| + values.each_with_index do |value, index| + assert result(name, index) =~ expr(value), "#{modes[index]}: #{name}" + end +end +*--#] PolyRatFunExpandModulus : +*--#[ PolyFunModulusCoefficient : +* Compact numeric arguments need the same reduction as general arguments. +Modulus 13; +CFunction pf; +PolyFun pf; +Local F = 6*pf(7); +Local G = 1/2*pf(7); +Local H = -6*pf(7); +Print +s; +.end +Modulus PlusMin,13; +CFunction pf; +PolyFun pf; +Local F = 6*pf(7); +Local G = 1/2*pf(7); +Local H = -6*pf(7); +Print +s; +.end +* A nonzero product can vanish when the modulus is composite. +Modulus 15; +CFunction pf; +PolyFun pf; +Local F = 3*pf(5); +Print; +.end +assert succeeded? +modes = ["Modulus 13", "PlusMin 13"] +[ + ["F", "+ pf(3)", "+ pf(3)"], + ["G", "+ pf(10)", "+ pf(-3)"], + ["H", "+ pf(10)", "+ pf(-3)"] +].each do |name, *values| + values.each_with_index do |value, index| + assert result(name, index) =~ expr(value), "#{modes[index]}: #{name}" + end +end +assert result("F", 2) =~ expr("0"), "Modulus 15: F" +*--#] PolyFunModulusCoefficient : +*--#[ DollarFactorizeModulusContent : +* Content extraction must reduce coefficients before sorting. +Modulus AlsoDollars,13; +Symbol x; +Local F = 1; +$p = 2*x+4*x^2; +FactDollar $p; +Print " N = %$;",$p[0]; +Print " F1 = %$;",$p[1]; +Print " F2 = %$;",$p[2]; +Print " F3 = %$;",$p[3]; +ModuleOption noparallel; +.end +Modulus PlusMin,AlsoDollars,13; +Symbol x; +Local F = 1; +$p = 2*x+4*x^2; +FactDollar $p; +Print " N = %$;",$p[0]; +Print " F1 = %$;",$p[1]; +Print " F2 = %$;",$p[2]; +Print " F3 = %$;",$p[3]; +ModuleOption noparallel; +.end +assert succeeded? +modes = ["Modulus 13", "PlusMin 13"] +[ + ["N", "3", "3"], + ["F1", "x", "x"], + ["F2", "x + 7", "x - 6"], + ["F3", "4", "4"] +].each do |name, *values| + values.each_with_index do |value, index| + assert result(name, index) =~ expr(value), "#{modes[index]}: #{name}" + end +end +*--#] DollarFactorizeModulusContent : diff --git a/doc/manual/statements.tex b/doc/manual/statements.tex index b6820631..cef4680a 100644 --- a/doc/manual/statements.tex +++ b/doc/manual/statements.tex @@ -3526,11 +3526,14 @@ \section{modulus} The options are \begin{description} \item[NoFunctions] Modulus calculus is not performed inside function -arguments. +arguments, except for the active PolyFun (see \ref{substapolyfun}) and the +coefficients of the active PolyRatFun in expand mode (see +\ref{substapolyratfun}). \item[AlsoFunctions] Modulus calculus is also performed inside function arguments. \item[CoefficientsOnly] Modulus calculus is neither performed inside function -arguments nor on powers of symbols. +arguments nor on powers of symbols. The exceptions listed under NoFunctions +also apply. \item[PlusMin] The values of numbers are reduced to the range $(-value+1)/2$ to $(value-1)/2$. \item[Positive] The values of numbers are reduced to the range $0$ to @@ -4632,6 +4635,11 @@ \section{polyratfun} addition of two coefficients is a simple and straightforward operation that does not need the expensive polynomial GCD computations. +\noindent When Modulus (see \ref{substamodulus}) is active, the coefficients +of this expansion are reduced modulo the specified value. This also applies +with the default NoFunctions option: the active PolyRatFun in expand mode +is exempt from this restriction, so AlsoFunctions is not required. + \noindent Of course one can program such expansions externally and maybe better suited for the problem at hand, but using this option of the polyratfun is much faster and gives fewer chances of mistakes. diff --git a/sources/proces.c b/sources/proces.c index b0fc8f74..b75501e1 100644 --- a/sources/proces.c +++ b/sources/proces.c @@ -5172,6 +5172,8 @@ int PrepPoly(PHEAD WORD *term,WORD par) if ( m >= AT.WorkSpace && m < AT.WorkTop ) AT.WorkPointer = m; if ( Normalize(BHEAD v) ) Terminate(-1); + /* Modulo reduction of the argument */ + if ( AN.ncmod != 0 && Modulus(v) ) Terminate(-1); AT.WorkPointer = oldworkpointer; m = w; if ( *v == 4 && v[2] == 1 && (v[1]&MAXPOSITIVE) == v[1] ) { diff --git a/sources/ratio.c b/sources/ratio.c index e1f29daa..5ca724cf 100644 --- a/sources/ratio.c +++ b/sources/ratio.c @@ -1340,6 +1340,16 @@ WORD *MultiplyWithTerm(PHEAD WORD *in, WORD *term, WORD par) *termout = tt - termout; // in or term can contain non-symbols: we call this in TakeContent Normalize(BHEAD termout); + /* Modulo reduction before sorting */ + if ( AN.ncmod != 0 ) { + if ( *termout != 0 && Modulus(termout) ) { + LowerSortLevel(); + AR.CompareRoutine = (COMPAREDUMMY)oldcompareroutine; + AR.SortType = oldsorttype; + goto CalledFrom; + } + if ( *termout == 0 ) { in += *in; continue; } + } StoreTerm(BHEAD termout); in += *in; } @@ -3221,6 +3231,10 @@ NormArg:; rr = r + *r; i = *r; rrr = outarg; NCOPY(rrr,r,i); Normalize(BHEAD outarg); + if ( AN.ncmod != 0 && *outarg != 0 && Modulus(outarg) ) { + LowerSortLevel(); AT.TrimPower = 0; + error = 3; goto onerror; + } if ( *outarg > 0 ) StoreTerm(BHEAD outarg); } r = fun+FUNHEAD+ARGHEAD; @@ -3301,6 +3315,10 @@ NormArg:; } } else { error = 2; goto onerror; } /* should not happen! */ + if ( AN.ncmod != 0 && *rr != 0 ) { + if ( Modulus(rr) ) { error = 3; goto onerror; } + r = rr + *rr; *r = 0; + } } else { /* Multi-term numerator. */ m = arg1+ARGHEAD; @@ -3332,6 +3350,12 @@ NormArg:; m++; while ( m < mm ) *r++ = *m++; *rrr = r-rrr; Normalize(BHEAD rrr); + if ( AN.ncmod != 0 ) { + if ( *rrr != 0 && Modulus(rrr) ) { + LowerSortLevel(); error = 3; goto onerror; + } + if ( *rrr == 0 ) continue; + } StoreTerm(BHEAD rrr); } EndSort(BHEAD rr,1); @@ -3413,6 +3437,12 @@ NormArg:; if ( r < AT.WorkTop && r >= AT.WorkSpace ) AT.WorkPointer = r; Normalize(BHEAD rrr); + if ( AN.ncmod != 0 ) { + if ( *rrr != 0 && Modulus(rrr) ) { + LowerSortLevel(); error = 3; goto onerror; + } + if ( *rrr == 0 ) continue; + } if ( ABS(rrr[*rrr-1]) == *rrr-1 ) { if ( AR.PolyFunPow >= 0 ) { StoreTerm(BHEAD rrr); @@ -3425,6 +3455,11 @@ NormArg:; } EndSort(BHEAD rr,1); } + /* Fast numerators bypass sorting, but still need reduction. */ + if ( *arg1 < 0 && AN.ncmod != 0 && *rr != 0 ) { + if ( Modulus(rr) ) { error = 3; goto onerror; } + rr[*rr] = 0; + } r = rr; while ( *r ) r += *r; i = r-rr; r = fun + FUNHEAD + ARGHEAD; @@ -3590,6 +3625,9 @@ NormArg:; *thecopy = rr - thecopy; AT.WorkPointer = rr; Normalize(BHEAD thecopy); + if ( AN.ncmod != 0 && *thecopy != 0 && Modulus(thecopy) ) { + LowerSortLevel(); error = 3; goto onerror; + } if ( *thecopy > 0 ) StoreTerm(BHEAD thecopy); AT.WorkPointer = thecopy; } @@ -3725,8 +3763,11 @@ int InvPoly(PHEAD WORD *inpoly, WORD maxpow, WORD sym) } } /* - Copy c to the proper location + Reduce before reusing c, and represent modular zeroes by null pointers. + Copy c to the proper location. */ + if ( AN.ncmod != 0 && lenc != 0 && + TakeModulus((UWORD *)c,&lenc,AC.cmod,AN.ncmod,UNPACK) ) goto calcerror; if ( lenc == 0 ) AT.pWorkSpace[outpointers+j] = 0; else { AT.pWorkSpace[outpointers+j] = w;