# How much pizza and how much frozen yogurt? ...with Gröbner bases

In a recent blog post I tried to get yacas to solve a system of polynomial equations. Unfortunately it could not do that, so I solved it numerically instead. Now it is possible – together with many other systems of polynomial equations thanks to fixing a small error in yacas. It has now been fixed, also in Ryacas (development version), so hurry up and update Ryacas to the latest version 0.9.0.9122.

First I will show how to do it directly in Ryacas (and afterwards show what happens behind the curtains).

library(Ryacas)

## Solving a system of equations

We define the system from the recent blog post in R:

L <- "x^2 * (y/4) - a*(3*x + 3*y/2 - 45)"
eq1 <- paste0("D(x) ", L) %>% yac_str()
eq2 <- paste0("D(y) ", L) %>% yac_str()
eq3 <- paste0("D(a) ", L) %>% yac_str()
eqs <- c(eq1, eq2, eq3)
eqs
## [1] "(x*y)/2-3*a"      "x^2/4-(3*a)/2"    "45-(3*x+(3*y)/2)"

As known from the theory of Gröbner bases (see e.g. the book Ideals, Varieties, and Algorithms by Cox, Little, Oshea), instead of finding roots in this entangled system, we can find the Gröbner bases and find roots for those instead:

eq_str <- paste0("{", paste0(eqs, " == 0", collapse = ", "), "}") %>%
y_fn("Solve", "{x, y, a}")
eq_str
## [1] "Solve({(x*y)/2-3*a == 0, x^2/4-(3*a)/2 == 0, 45-(3*x+(3*y)/2) == 0}, {x, y, a})"
sol_all <- eq_str %>% yac_str()
sol_all
## [1] "{{x==0,y==30,a==0},{x==10,y==10,a==50/3}}"

As seen, we get two solutions. The solution $$a = 0$$ is not feasible in our case (read the recent blog post) so we just pick the second:

sol_elible_vars <- sol_all %>% y_fn("Nth", 2) %>% yac_str()
sol_elible_vars
## [1] "{x==10,y==10,a==50/3}"

We have now finally arrived at the answer: $(x, y, a) = (10, 10, 50/3) .$

It is also possible to postprocess the answer a bit if that is needed:

# Remove variable names if wanted:
sol_elible <- sol_elible_vars %>% y_rmvars() %>% yac_str()
sol_elible
## [1] "{10,10,50/3}"
# And maybe convert to numbers:
sol_elible %>% yac_expr()
## expression(c(10, 10, 50/3))
# ...maybe even evaluated:
sol_elible %>% yac_expr() %>% eval()
## [1] 10.00000 10.00000 16.66667

## Enter Gröbner bases

We could actually have solved the system manually earlier using Gröbner bases.

As known from the theory of Gröbner bases (but I will not start trying to explain it because I do not know them nearly as well as I should), instead of finding roots in this entangled system, we can find the Gröbner bases and find roots for those instead:

fs <- paste0("{", paste0(eqs, collapse = ", "), "}")
fs
## [1] "{(x*y)/2-3*a, x^2/4-(3*a)/2, 45-(3*x+(3*y)/2)}"
gs <- fs %>% y_fn("Groebner") %>% yac_str() %>% as_r()
gs
## [1] "(x*y)/2-3*a"            "x^2/4+((-3)*a)/2"
## [3] "(-3)*x+((-3)*y)/2+45"   "9*y*a-90*a"
## [5] "15*y-y^2/2-6*a"         "(-6)*a^2+100*a"
## [7] "((-15)*y)/2+(-9)*a+225"

As seen the system is now “triangulated”. So we proceed:

yac_str("Solve((-6)*a^2+100*a == 0, a)")
## [1] "{a==0,a==50/3}"

As mentioned, the solution $$a = 0$$ is not feasible in our case (read the recent blog post), so we use $$a = 50/3$$ and obtain:

yac_str("WithValue(a, 50/3, Solve(9*y*a-90*a == 0, y))")
## [1] "{y==10}"

So $$y = 10$$. And we get

yac_str("WithValue({a, y}, {50/3, 10}, Solve((x*y)/2-3*a == 0, x))")
## [1] "{x==10}"

We have now finally arrived at the answer: $(x, y, a) = (10, 10, 50/3) .$

##### Mikkel Meyer Andersen
###### Assoc. Professor of Applied Statistics

My research interests include applied statistics and computational statistics.