Add absolute PolyRatFun expansion and fix relative mode

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 :

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^1>+...+<ep^10>,ep);
+Local F2 = rat(1+<ep^1>+...+<ep^10>,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 :

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 

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;
 		}

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)