Deriving Heron's formula with Sympy using Groebner basis

Question

I'm trying to implement with Sympy the procedure described in this lecture, where Groebner basis is used to derive Heron's formula, but I'm not able to implement the last step, the computation of the correct Groebner basis, surely because of my lack of understanding of the subject. Here is the code I wrote:

import sympy as sy
sy.init_printing(use_unicode=True, wrap_line=False)

def cross(a,b): # 2D cross product
    return a[0]*b[1]-a[1]*b[0]

n = 3 # number of vertices of triangle
# defining 3n+1 symbols
P = sy.Matrix(sy.MatrixSymbol('P', n, 2)) # matrix of coordinates of vertices
s = sy.symbols(f's0:{n}') # lengths of polygon' sides
A,R,cx,cy = sy.symbols('A R cx cy') # area, radius, center coordinates
C = sy.Matrix([[cx,cy]])
P[0,0] = 0
P[0,1] = 0
P[1,1] = 0
# defining 2n+1 equations
eq_area = 2*A - sum(map(cross,*zip(*((P[i,:],P[j,:]) for i,j in zip(range(n),list(range(1,n))+[0]))))) # area of polygon
eqs_vonc = [R**2-((C-P.row(r))*sy.transpose(C-P.row(r)))[0] for r in range(P.rows)] # vertices on circumference
eqs_sides = [s[i]**2-((P[i,:]-P[j,:])*sy.transpose(P[i,:]-P[j,:]))[0] for i,j in zip(range(n),list(range(1,n))+[0])] # sides lenghts
eqs = [eq_area]+eqs_sides+eqs_vonc
# compute Groebner basis
G = sy.groebner(eqs,A,R,*s,*C) # just tried

How should the last step be implemented in order to obtain the Heron's formula?

Answer 1

With Groebner bases you have to be very careful with symbol ordering and exactly which symbols are considered to be unknowns vs parameters. The call to Mathematica's GroebnerBasis asks to find a basis in some symbols while eliminating others. SymPy's groebner function does not have an explicit option for eliminating others but it can be done simply placing those symbols earlier in the list of symbols passed to groebner (earlier symbols are eliminated before later ones in lex ordering). Then we can just discard the equations in the basis that contain the symbols that we wanted to eliminate.

Here we want to eliminate cx, cy (*C) and P[1,0], P[2,0] and P[2,1] first so we put them at the beginning of the symbol list. After that we have the symbols whose equations we are interested in which are R and s0, s1, s2 (*s):

In [130]: G = groebner(eqs, [*C, P[1,0], P[2,0], P[2,1], R, *s])

In [133]: for eq in G: print(eq)
4*A*cx - s0**2*P[2, 1]
16*A**2*cy - 16*A**2*P[2, 1] + s0**2*s1**2*P[2, 1] + s0**2*s2**2*P[2, 1] - s1**4*P[2, 1] + 2*s1**2*s2**2*P[2, 1] - s2**4*P[2, 1]
2*A*P[1, 0] - s0**2*P[2, 1]
4*A*P[2, 0] - s0**2*P[2, 1] + s1**2*P[2, 1] - s2**2*P[2, 1]
-4*A**2 + s0**2*P[2, 1]**2
4*A**2*s0**2 - 8*A**2*s1**2 - 8*A**2*s2**2 + 16*A**2*P[2, 1]**2 + s1**4*P[2, 1]**2 - 2*s1**2*s2**2*P[2, 1]**2 + s2**4*P[2, 1]**2
16*A**2*R**2 - s0**2*s1**2*s2**2
16*A**2 + s0**4 - 2*s0**2*s1**2 - 2*s0**2*s2**2 + s1**4 - 2*s1**2*s2**2 + s2**4

The first few equations have symbols that we wanted to eliminate so we can skip them. Only the final two equations have everything eliminated. Either of those could be solved for A**2 in terms of some of the other symbols:

In [134]: solve(G[-1], A**2)[0].factor()
Out[134]: 
-(s₀ - s₁ - s₂)⋅(s₀ - s₁ + s₂)⋅(s₀ + s₁ - s₂)⋅(s₀ + s₁ + s₂) 
─────────────────────────────────────────────────────────────
                              16                             

In [135]: solve(G[-2], A**2)[0].factor()
Out[135]: 
  2   2   2
s₀ ⋅s₁ ⋅s₂ 
───────────
       2   
   16⋅R 

I don't know if either of those is what you are looking for.