[R] [External Help] Multivariate Polynomials in R

[Leonard Mada]

Dear R Users,

I have written some R code for multivariate polynomials in R. I am 
looking forward for some help in redesigning and improving the code.

Although this code was not planned initially to be released as a 
package, the functionality has become quite versatile over time. I will 
provide some examples below. If anyone is interested in multivariate 
polynomials and has some spare time, or has some students interested to 
learn some interesting math, feel free to contact me.

The immediate focus should be on:
1) Writing/improving the automatic tests;
2) Redesigning the code (and build an R package);

As the code has grown in size, I am very cautious to change anything, 
until proper tests are written. I have started to write some test code 
(the link to the GitHub page is below), but I am not yet very confident 
how to properly write the tests and also lack the time as well. I will 
appreciate any expertise and help on this topic.

Ultimately, I hope to be able to focus more on the math topics. I will 
post a separate call for some of these topics.

CODE DETAILS

The source files are on GitHub:
https://github.com/discoleo/R/blob/master/Math/Polynomials.Helper.R
- (all) files named Polynomials.Helper.XXX.R are needed; (~ 25 files, 
including the test files);
- if requested, I can also upload a zip file with all these source files;
- the code started as some helper scripts (which is why all those files 
are mixed with other files);

The multivariate polynomials are stored as data.frames and R's 
aggregate() function is the workhorse: but it proved sufficiently fast 
and the code works well even with polynomials with > 10,000 monomials. I 
have some older Java code which used a TreeMap (sorted map), but I do 
not maintain that code anymore. I was very reserved initially regarding 
the efficiency of the data frame; but it worked well! And it proved very 
useful for sub-setting specific monomials!

I have attached some concrete examples below.

Sincerely,

Leonard


### Examples

source("Polynomials.Helper.R")
# - requires also the other Helper scripts;
# - not strictly needed (but are loaded automatically):
#   library(polynom)
#   library(pracma)

### Example 1:
n = 2; # Power "n" will be evaluated automatically
p1 = toPoly.pm("x^n*y + b*z - R")
p2 = toPoly.pm("y^n*z + b*x - R")
p3 = toPoly.pm("z^n*x + b*y - R")

pR = solve.lpm(p1, p2, p3, xn=c("z", "y"))
str(pR) # 124 monomials
# tweaking manually can improve the results;
pR = solve.lpm(p1, p2, p3, xn=c("y", "z"))
str(pR)
# pR[[2]]$Rez: 19 monomials: much better;

pR2 = div.pm(pR[[2]]$Rez, "x^3 + b*x - R", "x")
# Divisible!
str(pR2)
# Order 12 polynomial in x (24 monomials);

### Note:
# - the P[12] contains the distinct roots:
#   it is the minimal order polynomial;
# - the trivial solution (x^3 + b*x = R) was removed;
# - this is the naive way to solve this system (but good as Demo);

# print the coefficients of x;
# (will be used inside the function coeff.S3Ht below)
pR2 = pR2$Rez;
pR2$coeff = - pR2$coeff; # positive leading coeff;
toCoeff(pR2, "x")

### Quick Check
solve.S3Ht = function(R, b) {
     coeff = coeff.S3Ht(R, b);
     x = roots(coeff); # library(pracma)
     # Note: pracma uses leading to free coeff;
     z = b*x^11 - R*x^10 - 2*R^2*b*x^5 + 2*R^2*b^2*x^3 + R*b^4*x^2 - R*b^5;
     z = z / (- R^2*x^6 - R*b^2*x^5 + 3*R*b^3*x^3 - b^6);
     y = (R - z^2*x) / b;
     sol = cbind(x, y, z);
     return(sol);
}
coeff.S3Ht = function(R, b) {
     coeff = c(b^2, - 2*R*b, R^2 - b^3, 3*R*b^2,
         - 3*R^2*b + b^4, R^3 - 4*R*b^3,
         2*R^2*b^2 - b^5, 5*R*b^4,
         R^4 - R^2*b^3 + b^6, - 3*R^3*b^2 - 3*R*b^5,
         - R^4*b + 3*R^2*b^4 - b^7,
         2*R^3*b^3 - R*b^6,
         - R^2*b^5 + b^8);
     return(coeff);
}

R = 5; b = -2;
sol = solve.S3Ht(R, b)
# all 12 sets of solutions:
x = sol[,1]; y = sol[,2]; z = sol[,3];

### Test:
x^2*y + b*z
y^2*z + b*x
z^2*x + b*y

id = 1;
eval.pm(p1, list(x=x[id], y=y[id], z=z[id], b=b, R=R))


##############

### Example 2:

n = 5
p1 = toPoly.pm("(x + a)^n + (y + a)^n - R1")
p2 = toPoly.pm("(x + b)*(y + b) - R2")

# Very Naive way to solve:
pR = solve.pm(p1, p2, "y")
str(pR)
table(pR$Rez$x)
# Order 10 with 109 monomials;
# [very naive!]