diff --git a/check/features.frm b/check/features.frm index 55f8317e..07c48f6d 100644 --- a/check/features.frm +++ b/check/features.frm @@ -3340,3 +3340,156 @@ Warning: ignoring out-of-place setup command: EOF assert result("test") =~ expr("10") *--#] warnings_2 : +*--#[ PolyRatFunExpandFixed : +#- +CFunction rat; +Symbol ep,a; +PolyRatFun rat(expand,ep,4,fixed); + +Local F1 = rat(1,ep^2-ep^4); +Local F2 = rat(ep^2,1-ep); +Local F3 = rat(1+ep^4+ep^5); +Local F4 = rat(a^3+a^2*ep^4+a*ep^5,16); +Local F5 = rat(1+a^3*ep^5+ep^6,ep); +Local F6 = rat(1,1-ep)^2; +Local F7 = rat(0,1-ep); +Local F8 = rat(ep^5,1-ep); +Local F9 = rat(1,1-ep); +Local F10 = 2; +Local F11 = 2/5; +Local F12 = 12345678901234567890; +Print F1,...,F8; +.sort + +* Multiply expanded functions and pull scalar coefficients into rat. +Drop F1,...,F8; +Multiply 2/3*rat(1,1-ep); +Print; +.sort + +Multiply rat(1,1-ep); +Print F9; +.end +#- +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,4,fixed); + +Local L = rat(ep^2+ep^3+ep^4); +.sort + +Multiply rat(ep^-1+1+ep); +Print; +.end +#- +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,0,fixed); + +Local Z1 = rat(ep); +Local Z2 = rat(ep,2); +Local Z3 = rat(ep^2,ep); +Local Z4 = rat(1,ep-ep^2); +Print; +.end +#- +CFunction rat; +Symbol ep,a; +PolyRatFun rat(expand,ep,-3,fixed); + +Local N1 = rat(1,ep^5-ep^6); +Local N2 = rat(a^3*ep^-5+2*a*ep^-3+ep^-2+a^2+ep,16); +Local N3 = rat(ep^-4+a*ep^-2+a^2+ep,ep); +Local N4 = rat(ep^-5+2*ep^-3+ep^-2+1); +Local N5 = rat(1,ep^5-ep^6)^2; +Local N6 = rat(1,ep^2); +Local N7 = rat(1,ep^3); +Local N8 = rat(1,16)*rat(1,16)+rat(ep,ep)*rat(ep,ep); +Local N9 = rat(1,ep)+rat(a,ep)+rat(ep,ep)+rat(1,16) + +rat(a,16)+rat(ep,16)+rat(1)+rat(a)+rat(ep); +.sort + +* These products use the already truncated N1, unlike N5 above. +Local N10 = N1*N1; +Local N11 = 4294967297/4294967299*N1*N1; +Print; +.end +#- +CFunction rat; +Symbol ep,a; +PolyRatFun rat(expand,ep,-1,fixed); + +Local S1 = rat(1,ep); +Local S2 = rat(a,ep); +Print; +.end +assert succeeded? +assert result("F1") =~ expr("rat(1 + ep^-2 + ep^2 + ep^4)") +assert result("F2") =~ expr("rat(ep^2 + ep^3 + ep^4)") +assert result("F3") =~ expr("rat(1 + ep^4)") +assert result("F4") =~ expr("rat(1/16*a^3 + 1/16*ep^4*a^2)") +assert result("F5") =~ expr("rat(ep^-1 + ep^4*a^3)") +assert result("F6") =~ expr("rat(1 + 2*ep + 3*ep^2 + 4*ep^3 + 5*ep^4)") +assert result("F7") =~ expr("0") +assert result("F8") =~ expr("0") +assert result("F9", 0) =~ expr("rat(2/3 + 4/3*ep + 2*ep^2 + 8/3*ep^3 + 10/3*ep^4)") +assert result("F9", 1) =~ expr("rat(2/3 + 2*ep + 4*ep^2 + 20/3*ep^3 + 10*ep^4)") +assert result("F10") =~ expr("rat(4/3 + 4/3*ep + 4/3*ep^2 + 4/3*ep^3 + 4/3*ep^4)") +assert result("F11") =~ expr("rat(4/15 + 4/15*ep + 4/15*ep^2 + 4/15*ep^3 + 4/15*ep^4)") +assert result("F12") =~ expr(" + rat(8230452600823045260 + 8230452600823045260*ep + + 8230452600823045260*ep^2 + 8230452600823045260*ep^3 + + 8230452600823045260*ep^4) +") +assert result("L") =~ expr("rat(ep + 2*ep^2 + 3*ep^3 + 2*ep^4)") +assert result("Z1") =~ expr("0") +assert result("Z2") =~ expr("0") +assert result("Z3") =~ expr("0") +assert result("Z4") =~ expr("rat(1 + ep^-1)") +assert result("N1") =~ expr("rat(ep^-5 + ep^-4 + ep^-3)") +assert result("N2") =~ expr("rat(1/16*ep^-5*a^3 + 1/8*ep^-3*a)") +assert result("N3") =~ expr("rat(ep^-5 + ep^-3*a)") +assert result("N4") =~ expr("rat(ep^-5 + 2*ep^-3)") +assert result("N5") =~ expr("rat(ep^-10 + 2*ep^-9 + 3*ep^-8 + 4*ep^-7 + 5*ep^-6 + 6*ep^-5 + 7*ep^-4 + 8*ep^-3)") +assert result("N6") =~ expr("0") +assert result("N7") =~ expr("rat(ep^-3)") +assert result("N8") =~ expr("0") +assert result("N9") =~ expr("0") +assert result("N10") =~ expr("rat(ep^-10 + 2*ep^-9 + 3*ep^-8 + 2*ep^-7 + ep^-6)") +assert result("N11") =~ expr(" + rat(4294967297/4294967299*ep^-10 + 8589934594/4294967299*ep^-9 + + 12884901891/4294967299*ep^-8 + 8589934594/4294967299*ep^-7 + + 4294967297/4294967299*ep^-6) +") +assert result("S1") =~ expr("rat(ep^-1)") +assert result("S2") =~ expr("rat(ep^-1*a)") +*--#] PolyRatFunExpandFixed : +*--#[ PolyRatFunExpandNegativeRelative : +#- +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,-3,relative); +.end +assert compile_error?("Illegal option field in PolyRatFun statement.") +*--#] PolyRatFunExpandNegativeRelative : +*--#[ PolyRatFunOptions : +#- +CFunction rat; +Symbol ep; +PolyRatFun rat(expand,ep,4,relative); + +Local F = rat(1,ep^2-ep^3); +Print; +.end +#- +CFunction rat; +Symbol ep; +PolyRatFun rat(finddivergence,ep); + +Local G = rat(13+ep,ep^2-ep^3); +Print; +.end +assert succeeded? +assert result("F") =~ expr("rat(1 + ep^-2 + ep^-1 + ep + ep^2)") +assert result("G") =~ expr("rat(1,ep^2)") +*--#] PolyRatFunOptions : diff --git a/check/fixes.frm b/check/fixes.frm index 94c7198c..b2882617 100644 --- a/check/fixes.frm +++ b/check/fixes.frm @@ -5613,3 +5613,54 @@ Print; assert succeeded? assert result("test") =~ expr("0") *--#] PullReq860_3 : +*--#[ PolyRatFunExpandDenominators : +#- +CFunction rat; +Symbol ep,a,b; +PolyRatFun rat(expand,ep,5); + +* Relative depth starts at the lowest power, including for simple denominators. +Local F1 = rat(1++...+,ep); +Local F2 = rat(1++...+,ep-ep^2); +Local F3 = rat(a^3,ep); +Local F4 = rat(a^3*ep,ep-ep^3); +Local F5 = rat(a^3,16); +Local F6 = rat(a,-16); +Local F7 = rat(a^3+a^2*ep+a*b*ep^3+b^2*ep^6,16); +Local F8 = rat(a^3*ep^-2+a*b+b^2*ep^3+a^2*ep^4,16); +Local F9 = rat(a^3*ep^3+a*b*ep^4+b^2*ep^9,ep); +Local F10 = rat(a^3*ep^2+a*b*ep^4+b^2*ep^9,16); +Print; +.end +#- +CFunction rat; +Symbol a,b,ep; +PolyRatFun rat(expand,ep,0); + +* Missing ep has power zero; other symbols are coefficients in either order. +Local Z1 = rat(1+ep+ep^2,ep); +Local Z2 = rat(a^3+ep+b^2*ep^2,ep); +Local Z3 = rat(a^3*ep^3+a*b*ep^4,ep); +Local Z4 = rat(a^3*ep^-2+a*b+ep,16); +Local Z5 = rat(a^3+b^2,16); +Local Z6 = rat(a^3*ep^2+a*b*ep^4,16); +Print; +.end +assert succeeded? +assert result("F1") =~ expr("rat(1 + ep^-1 + ep + ep^2 + ep^3 + ep^4)") +assert result("F2") =~ expr("rat(2 + ep^-1 + 3*ep + 4*ep^2 + 5*ep^3 + 6*ep^4)") +assert result("F3") =~ expr("rat(ep^-1*a^3)") +assert result("F4") =~ expr("rat(a^3 + ep^2*a^3 + ep^4*a^3)") +assert result("F5") =~ expr("rat(1/16*a^3)") +assert result("F6") =~ expr("rat(-1/16*a)") +assert result("F7") =~ expr("rat(1/16*a^3 + 1/16*ep*a^2 + 1/16*ep^3*a*b)") +assert result("F8") =~ expr("rat(1/16*ep^-2*a^3 + 1/16*a*b + 1/16*ep^3*b^2)") +assert result("F9") =~ expr("rat(ep^2*a^3 + ep^3*a*b)") +assert result("F10") =~ expr("rat(1/16*ep^2*a^3 + 1/16*ep^4*a*b)") +assert result("Z1") =~ expr("rat(ep^-1)") +assert result("Z2") =~ expr("rat(a^3*ep^-1)") +assert result("Z3") =~ expr("rat(a^3*ep^2)") +assert result("Z4") =~ expr("rat(1/16*a^3*ep^-2)") +assert result("Z5") =~ expr("rat(1/16*b^2 + 1/16*a^3)") +assert result("Z6") =~ expr("rat(1/16*a^3*ep^2)") +*--#] PolyRatFunExpandDenominators : diff --git a/doc/manual/statements.tex b/doc/manual/statements.tex index b6820631..2745ac80 100644 --- a/doc/manual/statements.tex +++ b/doc/manual/statements.tex @@ -4547,6 +4547,7 @@ \section{polyratfun} Type & Declaration statement\\ Syntax & polyratfun {\tt<}name of function{\tt>}; \\ & polyratfun {\tt<}name of function{\tt>},{\tt<}name of function{\tt>}; \\ + & polyratfun {\tt<}name of function{\tt>}({\tt<}options{\tt>}); \\ & polyratfun; \\ See also & polyfun (\ref{substapolyfun}), moduleoption (\ref{substamoduleoption}) @@ -4618,19 +4619,44 @@ \section{polyratfun} the end of the calculation one can see how many poles in x were maximally present, and hence how far one has to expand in x. Because the contents of the polyratfun are extremely simple, the expensive rational arithmetic is -completely absent and things should go rather fast. +completely absent and things should go rather fast. +The keyword \texttt{finddivergence} is a synonym for +\texttt{divergence}, so the following declaration has the same effect: +\begin{verbatim} + PolyRatFun rat(finddivergence,x); +\end{verbatim} \noindent In the second option one can specify how far one should expand: \begin{verbatim} PolyRatFun rat(expand,x,power); + PolyRatFun rat(expand,x,power,relative); +\end{verbatim} +These declarations are equivalent: \texttt{relative} is the default. The +denominator can only be a polynomial in x. Its inverse will be expanded and +multiplied by the numerator. The expansion depth \texttt{power} is a relative, +nonnegative depth, meaning that if the rational function starts at $x^p$, terms +through $x^{p+\mathrm{power}}$ are retained and higher powers are discarded. +For example, with a leading power of $x^{-2}$, a depth of 4 retains terms +through $x^2$. + +\noindent To specify an absolute maximum exponent, use +\begin{verbatim} + PolyRatFun rat(expand,x,power,fixed); +\end{verbatim} +In this case all powers of x greater than \texttt{power} are discarded, +independently of the leading power. The exponent \texttt{power} may also be +negative. For example, with +\begin{verbatim} + PolyRatFun rat(expand,x,-3,fixed); \end{verbatim} -In this case the denominator can only be a polynomial in the variable x. It -will be expanded and multiplied by the numerator and eventually all terms -with powers of x that are greater than 'power' will be discarded. The -remaining incidence of the function rat will then have only one argument, -like the polyfun (see \ref{substapolyfun}). The advantage is that now the -addition of two coefficients is a simple and straightforward operation that -does not need the expensive polynomial GCD computations. +the function \texttt{rat(1,x\^{}5-x\^{}6)} becomes +\texttt{rat(x\^{}-5+x\^{}-4+x\^{}-3)}. + +\noindent In either expansion mode, the remaining incidence of the function +rat has only one argument, like the polyfun (see \ref{substapolyfun}). The +advantage is that now the addition of two coefficients is a simple and +straightforward operation that does not need the expensive polynomial GCD +computations. \noindent Of course one can program such expansions externally and maybe better suited for the problem at hand, but using this option of the diff --git a/sources/compcomm.c b/sources/compcomm.c index 4c281613..ae998ecd 100644 --- a/sources/compcomm.c +++ b/sources/compcomm.c @@ -5515,7 +5515,7 @@ int CoPolyRatFun(UBYTE *s) symbols[AC.lPolyFunVar].maxpower = MAXPOWER; } else if ( StrICmp(s,(UBYTE *)"expand") == 0 ) { - WORD x = 0, etype = 2; + WORD x = 0, etype = 2, sign = 1; if ( c != ',' ) { MesPrint("&Illegal option field in PolyRatFun statement."); return(1); @@ -5532,6 +5532,10 @@ int CoPolyRatFun(UBYTE *s) } *t = c; while ( *t == ',' || *t == ' ' || *t == '\t' ) t++; + if ( *t == '-' ) { + sign = -1; + t++; + } if ( *t > '9' || *t < '0' ) { MesPrint("&Illegal option field in PolyRatFun statement."); return(1); @@ -5560,9 +5564,14 @@ int CoPolyRatFun(UBYTE *s) return(1); } } + // Relative expansion depth cannot have a negative power. + if ( etype == 2 && sign == -1 ) { + MesPrint("&Illegal option field in PolyRatFun statement."); + return(1); + } AR.PolyFunExp = AC.lPolyFunExp = etype; AR.PolyFunVar = AC.lPolyFunVar; - AR.PolyFunPow = AC.lPolyFunPow = x; + AR.PolyFunPow = AC.lPolyFunPow = sign*x; symbols[AC.lPolyFunVar].minpower = -MAXPOWER; symbols[AC.lPolyFunVar].maxpower = MAXPOWER; } diff --git a/sources/proces.c b/sources/proces.c index b0fc8f74..82ea6292 100644 --- a/sources/proces.c +++ b/sources/proces.c @@ -5011,7 +5011,7 @@ int PrepPoly(PHEAD WORD *term,WORD par) But because this is an undocumented feature for very special purposes, we don't do anything about it. (30-aug-2011) */ - if ( AR.PolyFunType == 2 && AR.PolyFunExp != 2 ) { + if ( AR.PolyFunType == 2 && AR.PolyFunExp != 2 && AR.PolyFunExp != 3 ) { WORD oldtype = AR.SortType; AR.SortType = SORTHIGHFIRST; if ( poly_ratfun_normalize(BHEAD term) != 0 ) Terminate(-1); @@ -5035,7 +5035,8 @@ int PrepPoly(PHEAD WORD *term,WORD par) i = ABS(m[-1]); if ( par > 0 ) { if ( count == 0 ) return(0); - else if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 && AR.PolyFunExp == 2) ) + else if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 + && (AR.PolyFunExp == 2 || AR.PolyFunExp == 3)) ) goto DoOne; else if ( AR.PolyFunType == 2 ) goto DoTwo; @@ -5049,7 +5050,8 @@ int PrepPoly(PHEAD WORD *term,WORD par) poly = t = tstop; if ( i == 3 && m[-2] == 1 && (m[-3]&MAXPOSITIVE) == m[-3] ) { *m++ = AR.PolyFun; - if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 && AR.PolyFunExp == 2) ) { + if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 + && (AR.PolyFunExp == 2 || AR.PolyFunExp == 3)) ) { *m++ = FUNHEAD+2; FILLFUN(m) *m++ = -SNUMBER; @@ -5069,7 +5071,8 @@ int PrepPoly(PHEAD WORD *term,WORD par) else { WORD *vm; r = tstop; - if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 && AR.PolyFunExp == 2) ) { + if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 + && (AR.PolyFunExp == 2 || AR.PolyFunExp == 3)) ) { *m++ = AR.PolyFun; *m++ = FUNHEAD+ARGHEAD+i+1; FILLFUN(m) @@ -5123,7 +5126,8 @@ int PrepPoly(PHEAD WORD *term,WORD par) #] Create a PolyFun : */ } - else if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 && AR.PolyFunExp == 2) ) { + else if ( AR.PolyFunType == 1 || (AR.PolyFunType == 2 + && (AR.PolyFunExp == 2 || AR.PolyFunExp == 3)) ) { DoOne:; /* #[ One argument : @@ -5472,8 +5476,8 @@ int PolyFunMul(PHEAD WORD *term) return(0); } ReStart: - if ( AR.PolyFunType == 2 && ( ( AR.PolyFunExp != 2 ) - || ( AR.PolyFunExp == 2 && AN.PolyNormFlag > 1 ) ) ) { + if ( AR.PolyFunType == 2 && ( ( AR.PolyFunExp != 2 && AR.PolyFunExp != 3 ) + || AN.PolyNormFlag > 1 ) ) { WORD count1 = 0, count2 = 0, count3; WORD oldtype = AR.SortType; t = term + 1; t1 = term + *term; t1 -= ABS(t1[-1]); @@ -5671,6 +5675,7 @@ int PolyFunMul(PHEAD WORD *term) */ w = AT.WorkPointer; NewSort(BHEAD0); + if ( AR.PolyFunType == 2 && AR.PolyFunExp == 3 ) AT.TrimPower = 1; if ( AR.PolyFunType == 2 && AR.PolyFunExp == 2 ) { AT.TrimPower = 1; /* @@ -5825,6 +5830,7 @@ PolyCall:; MLOCK(ErrorMessageLock); PolyCall2:; AR.PolyFunPow = OldPolyFunPow; + if ( AR.PolyFunExp == 3 ) AT.TrimPower = 0; MesCall("PolyFunMul"); MUNLOCK(ErrorMessageLock); SETERROR(-1) diff --git a/sources/ratio.c b/sources/ratio.c index e1f29daa..29f06b41 100644 --- a/sources/ratio.c +++ b/sources/ratio.c @@ -3221,7 +3221,17 @@ NormArg:; rr = r + *r; i = *r; rrr = outarg; NCOPY(rrr,r,i); Normalize(BHEAD outarg); - if ( *outarg > 0 ) StoreTerm(BHEAD outarg); + if ( *outarg > 0 ) { + if ( AR.PolyFunExp == 3 && AR.PolyFunPow < 0 ) { + if ( ABS(outarg[*outarg-1]) == *outarg-1 + || outarg[1] != SYMBOL ) continue; + for ( i = 3; i < outarg[2]; i += 2 ) { + if ( outarg[i] == AR.PolyFunVar ) break; + } + if ( i >= outarg[2] ) continue; + } + StoreTerm(BHEAD outarg); + } } r = fun+FUNHEAD+ARGHEAD; EndSort(BHEAD r,1); @@ -3257,7 +3267,15 @@ NormArg:; } } else { error = 5; goto onerror; } if ( arg2 == NULL ) { - if ( *arg1 < 0 ) { fun[2] = CLEANFLAG; goto done; } + if ( *arg1 < 0 ) { + if ( AR.PolyFunExp == 3 && ( ( AR.PolyFunPow < 0 && (fun[2]&MUSTCLEANPRF) != 0 + && ( *arg1 == -SNUMBER || *arg1 == -SYMBOL ) ) + || ( AR.PolyFunPow == 0 && *arg1 == -SYMBOL + && arg1[1] == AR.PolyFunVar ) ) ) { + *arg1 = -SNUMBER; arg1[1] = 0; + } + fun[2] = CLEANFLAG; goto done; + } if ( fun[2] == CLEANFLAG ) goto done; goto NormArg; /* Note: should not come here */ } @@ -3265,6 +3283,19 @@ NormArg:; Produce the output argument in outarg */ if ( outarg == NULL ) outarg = TermMalloc("ExpandRat")+ARGHEAD; + if ( AR.PolyFunExp == 3 && AR.PolyFunPow < 0 && AN.PolyNormFlag != 0 + && (fun[2]&MUSTCLEANPRF) == 0 + && ( *arg1 == -SNUMBER || ( *arg1 > 0 && *arg1 == ARGHEAD+arg1[ARGHEAD] + && arg1[ARGHEAD] == ABS(arg1[*arg1-1])+1 ) ) + && ( *arg2 == -SNUMBER || ( *arg2 > 0 && *arg2 == ARGHEAD+arg2[ARGHEAD] + && arg2[ARGHEAD] == ABS(arg2[*arg2-1])+1 ) ) ) { +/* + With other polyratfun factors, rational normalization inserts a clean + numeric factor containing the term's coefficient. Keep it until it is + multiplied into the series. +*/ + AR.PolyFunPow = 0; + } if ( *arg2 <= 0 ) { /* @@ -3301,6 +3332,8 @@ NormArg:; } } else { error = 2; goto onerror; } /* should not happen! */ + if ( AR.PolyFunExp == 3 && AR.PolyFunPow < + ( *arg1 == -SYMBOL && arg1[1] == AR.PolyFunVar ? 0 : -1 ) ) *rr = 0; } else { /* Multi-term numerator. */ m = arg1+ARGHEAD; @@ -3316,8 +3349,12 @@ NormArg:; if ( minpow2 > 0 ) minpow2 = 0; } else { - if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) { - if ( rrm[4] < minpow2 ) minpow2 = rrm[4]; + if ( rrm[1] == SYMBOL ) { + for ( i = 3; i < rrm[2]; i += 2 ) { + if ( rrm[i] == AR.PolyFunVar ) break; + } + eppow = i < rrm[2] ? rrm[i+1] : 0; + if ( eppow < minpow2 ) minpow2 = eppow; } else { error = 5; goto onerror; } } @@ -3325,6 +3362,7 @@ NormArg:; } AR.PolyFunPow += minpow2-1; } + AT.TrimPower = 1; while ( m < arg2 ) { r = outarg; rrr = r++; mm = m + *m; @@ -3332,8 +3370,19 @@ NormArg:; m++; while ( m < mm ) *r++ = *m++; *rrr = r-rrr; Normalize(BHEAD rrr); + if ( *rrr == 0 ) continue; + if ( AR.PolyFunPow < 0 ) { + // Terms independent of AR.PolyFunVar were not discarded by Normalize + // We remove these now. + if ( ABS(rrr[*rrr-1]) == *rrr-1 ) continue; + for ( i = 3; i < rrr[2]; i += 2 ) { + if ( rrr[i] == AR.PolyFunVar ) break; + } + if ( i >= rrr[2] ) continue; + } StoreTerm(BHEAD rrr); } + AT.TrimPower = 0; EndSort(BHEAD rr,1); r = rr; while ( *r ) r += *r; } @@ -3358,7 +3407,9 @@ NormArg:; rr = r = outarg; if ( arg2[1] == 0 ) { error = 1; goto onerror; } if ( *arg1 == -SNUMBER ) { /* Things may not be normalized */ - if ( arg1[1] == 0 ) { *r++ = 0; } + if ( arg1[1] == 0 || ( AR.PolyFunExp == 3 && AR.PolyFunPow < 0 ) ) { + *r++ = 0; + } else { co1[0] = ABS(arg1[1]); co1[1] = 1; co2[0] = 1; co2[1] = ABS(arg2[1]); @@ -3377,6 +3428,8 @@ NormArg:; if ( arg2[1] < 0 ) *r++ = -3; else *r++ = 3; *r = 0; + if ( AR.PolyFunExp == 3 && ( AR.PolyFunPow < 0 + || ( AR.PolyFunPow == 0 && arg1[1] == AR.PolyFunVar ) ) ) *rr = 0; } else if ( *arg1 < 0 ) { error = 2; goto onerror; } else { /* Multi-term numerator. */ @@ -3393,8 +3446,12 @@ NormArg:; if ( minpow2 > 0 ) minpow2 = 0; } else { - if ( rrm[1] == SYMBOL && rrm[2] == 4 && rrm[3] == AR.PolyFunVar ) { - if ( rrm[4] < minpow2 ) minpow2 = rrm[4]; + if ( rrm[1] == SYMBOL ) { + for ( i = 3; i < rrm[2]; i += 2 ) { + if ( rrm[i] == AR.PolyFunVar ) break; + } + eppow = i < rrm[2] ? rrm[i+1] : 0; + if ( eppow < minpow2 ) minpow2 = eppow; } else { error = 5; goto onerror; } } @@ -3402,6 +3459,7 @@ NormArg:; } AR.PolyFunPow += minpow2; } + AT.TrimPower = 1; while ( m < arg2 ) { r = rr; rrr = r++; mm = m + *m; @@ -3413,16 +3471,20 @@ NormArg:; if ( r < AT.WorkTop && r >= AT.WorkSpace ) AT.WorkPointer = r; Normalize(BHEAD rrr); + if ( *rrr == 0 ) continue; if ( ABS(rrr[*rrr-1]) == *rrr-1 ) { - if ( AR.PolyFunPow >= 0 ) { - StoreTerm(BHEAD rrr); - } + if ( AR.PolyFunPow < 0 ) continue; } - else if ( rrr[1] == SYMBOL && rrr[2] == 4 && - rrr[3] == AR.PolyFunVar && rrr[4] <= AR.PolyFunPow ) { - StoreTerm(BHEAD rrr); + else if ( rrr[1] != SYMBOL ) continue; + else if ( AR.PolyFunPow < 0 ) { + for ( i = 3; i < rrr[2]; i += 2 ) { + if ( rrr[i] == AR.PolyFunVar ) break; + } + if ( i >= rrr[2] ) continue; } + StoreTerm(BHEAD rrr); } + AT.TrimPower = 0; EndSort(BHEAD rr,1); } r = rr; while ( *r ) r += *r; @@ -3541,7 +3603,13 @@ NormArg:; Hence there is no need to free memory afterwards. */ if ( AR.PolyFunExp == 3 ) { + if ( minpow > AR.PolyFunPow ) { + fun[1] = FUNHEAD+2; fun[2] = CLEANFLAG; + fun[FUNHEAD] = -SNUMBER; fun[FUNHEAD+1] = 0; + goto done; + } ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow-minpow,AR.PolyFunVar); + minpow = 0; /* The multiplication cutoff is absolute. */ } else { ipoly = InvPoly(BHEAD denominator,AR.PolyFunPow,AR.PolyFunVar);