Sullivant's algorithm for primes in a polynomial ring

Given ideals I and J in a polynomial ring, we compute their join I*J:

i1 : S = QQ[x,y,z];

i2 : I = ideal(x^3,x^2*y^2,x*z^3,y^4,y^2*z);

o2 : Ideal of S

i3 : J = joinIdeals(I,I)

             5   5    7   4 3   3 5   2 3 3   3 5
o3 = ideal (x , y z, y , x y , x z , x y z , x y )

o3 : Ideal of S

Following Seth Sullivant's "Combinatorial symbolic powers", J. Algebra 319 (2008), no. 1, 115–142, we can compute symbolic powers of prime ideals using join:

i4 : A = QQ[x,y,z];

i5 : symbolicPowerJoin(ideal(x,y),2)

             2        2
o5 = ideal (y , x*y, x )

o5 : Ideal of A

i6 : f = map(QQ[t],A,{t^3,t^4,t^5})

                      3   4   5
o6 = map (QQ[t], A, {t , t , t })

o6 : RingMap QQ[t] <-- A

i7 : P = ker f;

o7 : Ideal of A

i8 : symbolicPowerJoin(P,2)

             4       2     2 2   2 3    3       2 2      3   3 2    4     3 
o8 = ideal (y  - 2x*y z + x z , x y  - x y*z - y z  + x*z , x y  - x z - y z
     ------------------------------------------------------------------------
            2   5      3     2       3
     + x*y*z , x  + x*y  - 3x y*z + z )

o8 : Ideal of A