Alternative algorithm to compute the symbolic powers of a prime ideal in positive characteristic

Given a prime ideal P in a regular ring of positive characteristic, symbPowerPrimePosChar computes its symbolic powers. Unfortunately, this algorithm is slower than others.

i1 : R = ZZ/7[x,y,z]

o1 = R

o1 : PolynomialRing

i2 : P = ker map(ZZ/7[t],R,{t^3,t^4,t^5})

             2         2     2   3
o2 = ideal (y  - x*z, x y - z , x  - y*z)

o2 : Ideal of R

i3 : J = symbPowerPrimePosChar(P,2)

             4       2     2 2   2 3    3       2 2      3   3 2    4     3 
o3 = 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 )

o3 : Ideal of R

The symbolic powers of P do not coincide with its powers.

i4 : J == P^2

o4 = false

We can also test it a bit faster, without computing the symbolic powers of $P$.

i5 : isSymbolicEqualOrdinary(P,2)

o5 = false