Add mul_ function for polynomial multiplications. Close #49

check/features.frm

@@ -184,6 +184,22 @@ assert result("F6") == result("F61")
 assert result("F71") =~ expr("f(nosquare,q2)*functions(p1,p2,q2,N1_?,N1_?)")
 assert result("F72") =~ expr("f(N1_?,q2)*functions(p1,p2,q2,N1_?)*nosquare.nosquare")
 *--#] CoToTensor :
+*--#[ Issue49 :
+* Add mul_ function for polynomial multiplications
+Symbols x,y,z;
+#$p = (1+x+y+z)^4;
+#$q = $p+1;
+#$r = mul_($p,$q);
+L r1 = $r;
+L r2 = $p^2 + $p;
+.sort
+Drop;
+L Zero = r1 - r2;
+P;
+.end
+assert succeeded?
+assert result("Zero") =~ expr("0")
+*--#] Issue49 : 
 *--#[ Issue72 :
 * "Setups: PATHVALUE not yet implemented"
 #:incdir foo

sources/declare.h

@@ -1499,6 +1499,7 @@ extern WORD  *poly_gcd(PHEAD WORD *, WORD *, WORD);
 extern WORD  *poly_div(PHEAD WORD *, WORD *, WORD);
 extern WORD  *poly_rem(PHEAD WORD *, WORD *, WORD);
 extern WORD  *poly_inverse(PHEAD WORD *, WORD *);
+extern WORD  *poly_mul(PHEAD WORD *, WORD *);
 extern WORD *poly_ratfun_add(PHEAD WORD *, WORD *);
 extern int   poly_ratfun_normalize(PHEAD WORD *);
 extern int   poly_factorize_argument(PHEAD WORD *, WORD *);

sources/ftypes.h

@@ -457,8 +457,9 @@ typedef int (*TFUN1)();
 #define PUTFIRST 102
 #define PERMUTATIONS 103
 #define PARTITIONS 104
+#define MULFUNCTION 105
 
-#define MAXBUILTINFUNCTION 104
+#define MAXBUILTINFUNCTION 105
 #define FIRSTUSERFUNCTION 150
 
 /*

sources/inivar.h

@@ -230,6 +230,7 @@ static struct fixedfun {
 	,{"putfirst_"   ,1 ,0             ,0               ,0} /* PUTFIRST */
 	,{"perm_"       ,1 ,0             ,0               ,0} /* PERMUTATIONS */
 	,{"partitions_" ,1 ,0             ,0               ,0} /* PARTITIONS */
+	,{"mul_"        ,0 ,0             ,0               ,0} /* MULFUNCTION */
 };
 
 FIXEDSET fixedsets[] = {

sources/polywrap.cc

@@ -1752,6 +1752,45 @@ WORD *poly_inverse(PHEAD WORD *arga, WORD *argb) {
 
 /*
   	#] poly_inverse : 
+  	#[ poly_mul :
+*/
+
+WORD *poly_mul(PHEAD WORD *a, WORD *b) {
+
+#ifdef DEBUG
+	cout << "*** [" << thetime() << "]  CALL : poly_mul" << endl;
+#endif
+
+	// Extract variables
+	vector<WORD *> e;
+	e.push_back(a);
+	e.push_back(b);
+	poly::get_variables(BHEAD e, false, false);  // TODO: any performance effect by sort_vars=true?
+
+	// Convert to polynomials
+	poly pa(poly::argument_to_poly(BHEAD a, false, true));
+	poly pb(poly::argument_to_poly(BHEAD b, false, true));
+
+	// Check for modulus calculus
+	WORD modp = poly_determine_modulus(BHEAD true, true, "polynomial multiplication");
+	pa.setmod(modp, 1);
+
+	// multiplication
+	pa *= pb;
+
+	// convert to Form notation
+	WORD *res = (WORD *)Malloc1((pa.size_of_form_notation() + 1) * sizeof(WORD), "poly_mul");
+	poly::poly_to_argument(pa, res, false);
+
+	// clean up and reset modulo calculation
+	poly_free_poly_vars(BHEAD "AN.poly_vars_mul");
+	AN.ncmod = AC.ncmod;
+
+	return res;
+}
+
+/*
+  	#] poly_mul : 
   	#[ poly_free_poly_vars :
 */
 

sources/proces.c

@@ -1368,7 +1368,8 @@ Important: we may not have enough spots here
 					return(1);
 				}
               }
-              else if ( *t == DIVFUNCTION || *t == REMFUNCTION || *t == INVERSEFUNCTION ) {
+              else if ( *t == DIVFUNCTION || *t == REMFUNCTION
+                        || *t == INVERSEFUNCTION || *t == MULFUNCTION ) {
                 WORD *tf;
 				int todo = 1, numargs = 0;
                 tf = t + FUNHEAD;
@@ -1381,6 +1382,7 @@ Important: we may not have enough spots here
 					if ( *t == DIVFUNCTION ) AN.TeInFun = -9;
 					else if ( *t == REMFUNCTION ) AN.TeInFun = -10;
 					else if ( *t == INVERSEFUNCTION ) AN.TeInFun = -11;
+					else if ( *t == MULFUNCTION ) AN.TeInFun = -14;
 					AN.TeSuOut = 0;
 					AR.TePos = -1;
 					return(1);
@@ -3700,6 +3702,9 @@ AutoGen:	i = *AT.TMout;
 				case -13:
 					if ( DoPartitions(BHEAD term,level) ) goto GenCall;
 					break;
+				case -14:
+					if ( DIVfunction(BHEAD term,level,3) < 0 ) goto GenCall;
+					break;
 			}
 		}
 		else {

sources/ratio.c

@@ -2760,9 +2760,11 @@ WORD *PolyDiv(PHEAD WORD *a,WORD *b,char *text)
 	Note that the output can be just a number or many terms.
 	In case par == 0 the output is [arg1/arg2]
 	In case par == 1 the output is [arg1%arg2]
+	In case par == 2 the output is [inverse of arg1 modulus arg2]
+	In case par == 3 the output is [arg1*arg2]
 */
 
-WORD divrem[3] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION };
+WORD divrem[4] = { DIVFUNCTION, REMFUNCTION, INVERSEFUNCTION, MULFUNCTION };
 
 int DIVfunction(PHEAD WORD *term,WORD level,int par)
 {
@@ -2772,7 +2774,7 @@ int DIVfunction(PHEAD WORD *term,WORD level,int par)
 	WORD *proper1, *proper2, *proper3 = 0;
 	int numargs = 0, type1, type2, actionflag1, actionflag2;
 	WORD startebuf = cbuf[AT.ebufnum].numrhs;
-	if ( par < 0 || par > 2 ) {
+	if ( par < 0 || par > 3 ) {
 		MLOCK(ErrorMessageLock);
 		MesPrint("Internal error. Illegal parameter %d in DIVfunction.",par);
 		MUNLOCK(ErrorMessageLock);
@@ -2862,6 +2864,7 @@ divzero:;
 	if ( par == 0 )      proper3 = poly_div(BHEAD proper1, proper2,0);
 	else if ( par == 1 ) proper3 = poly_rem(BHEAD proper1, proper2,0);
 	else if ( par == 2 ) proper3 = poly_inverse(BHEAD proper1, proper2);
+	else if ( par == 3 ) proper3 = poly_mul(BHEAD proper1, proper2);
 	if ( proper3 == 0 ) goto CalledFrom;
 	if ( actionflag1 || actionflag2 ) {
 		if ( ( arg3 = TakeExtraSymbols(BHEAD proper3,startebuf) ) == 0 ) goto CalledFrom;

sources/token.c

@@ -948,6 +948,7 @@ int simp2token(SBYTE *s)
 						|| n == (DIVFUNCTION-FUNCTION)
 						|| n == (REMFUNCTION-FUNCTION)
 						|| n == (INVERSEFUNCTION-FUNCTION)
+						|| n == (MULFUNCTION-FUNCTION)
 						|| n == (FACTORIN-FUNCTION)
 						|| n == (FIRSTTERM-FUNCTION)
 						|| n == (CONTENTTERM-FUNCTION) )
@@ -1090,7 +1091,8 @@ tcommon:				v++; while ( *v >= 0 ) v++;
 						if ( n == GCDFUNCTION-FUNCTION
 						|| n == DIVFUNCTION-FUNCTION
 						|| n == REMFUNCTION-FUNCTION
-						|| n == INVERSEFUNCTION-FUNCTION ) {
+						|| n == INVERSEFUNCTION-FUNCTION
+						|| n == MULFUNCTION-FUNCTION ) {
 							*t = TEMPTY; s++;
 						}
 						else *fill++ = *s++;