Julia Code for sampling real cubics

Loading Packages

using Pkg;
using Distributions, Random;
using HomotopyContinuation;
Pkg.status()

Computing number of real lines

## Variables of the coordinate ring R[x,y,z,w] of P^3_R
## (real projective sapce of dim 3)
HomotopyContinuation.@var(x,y,z,w)



## Variables for the parametrization of a line
##  (s,t, as + bt, cs + dt)
## The other possible parametrization have 0 measure.
HomotopyContinuation.@var(a,b,c,d)



### Orthogonal basis of the space of spherical harmonics of degree 3
### together with the associated norms.
H3 = [  [x*y*z,1],
        [x*y*w,1],
        [x*z*w,1],
        [y*z*w,1],

        [x^3 - x*(y^2 + z^2 + w^2), 2*sqrt(3)],
        [y^3 - y*(x^2 + z^2 + w^2), 2*sqrt(3)],
        [z^3 - z*(x^2 + y^2 + w^2), 2*sqrt(3)],
        [w^3 - w*(x^2 + y^2 + z^2), 2*sqrt(3)],

        [x*(y^2 - z^2), 2],
        [y*(z^2 - w^2), 2],
        [z*(y^2 - x^2), 2],
        [w*(x^2 - y^2), 2],

        [x*(y^2 - w^2) - (1/2)*x*(y^2 - z^2), sqrt(3)],
        [y*(z^2 - x^2) - (1/2)*y*(z^2 - w^2), sqrt(3)],
        [z*(y^2 - w^2) - (1/2)*z*(y^2 - x^2), sqrt(3)],
        [w*(x^2 - z^2) - (1/2)*w*(x^2 - y^2), sqrt(3)]
    ];


### Orthogonal basis of the space of spherical harmonics of degree 1 multiplied with the euclidean norm.
### together with the associated norms in the space V.
xH1 = [
        [x*(x^2 + y^2 + z^2 + w^2), 4*sqrt(3)],
        [y*(x^2 + y^2 + z^2 + w^2), 4*sqrt(3)],
        [z*(x^2 + y^2 + z^2 + w^2), 4*sqrt(3)],
        [w*(x^2 + y^2 + z^2 + w^2), 4*sqrt(3)]
       ];


#####################################################################
######  Function that samples a real cubic from the distribution P_lambda
#####################################################################
function randomSmoothCubicSurface(lambda)

    randSurface = 0*x;

    for T in H3
        pol  = T[1]
        norm = T[2]
        randSurface =  randSurface + lambda*rand(Normal(0,1))* (pol/norm);

    end;

    for T in xH1
        pol  = T[1]
        norm = T[2]
        randSurface = randSurface + (1-lambda)*rand(Normal(0,1))*(pol/norm);

    end;

    return randSurface;

end;


####################################################################################################
######  Function that returns the lines counts in a sample of M surfaces under P_lambda
####################################################################################################
function realLineCounts(lambda,M)

    results = Dict([3=>0.0, 7=>0.0, 15=>0.0, 27=>0.0]);

    for i in 1:M

        bool = true

        while bool

            f = randomSmoothCubicSurface(lambda);

            g0 = subs(f, [x,y,z,w] => [a, c, 1, 0]);

            g1 = subs(f, [x,y,z,w] => [a+b, c+d, 1, 1]);

            g2 = subs(f, [x,y,z,w] => [a-b, c-d, 1, -1]);

            g3 = subs(f, [x,y,z,w] => [b, d, 0, 1]);


            S = System([g0,g1,g2,g3]);
            A = solve(S);

            r = size(real_solutions(A))[1];


            if !(r in [3,7,15,27])
                bool = true
            else
                bool = false;

                results[r] += 1;

            end;

        end;
    end;

    return results;
end;



M = 10000;


l = 1/100

results = realLineCounts(l, M);
println(results)

Tracking 81 paths... 100%|██████████████████████████████| Time: 0:00:18
  # paths tracked:                  81
  # non-singular solutions (real):  27 (3)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         27 (3)
Dict(15 => 0.0, 7 => 0.0, 27 => 0.0, 3 => 100.0)

########## Return the distribution \pi^{(\lambda)} for \lambda = i/100 as i ranges in [1,100]
### Warning: this takes a very long time, we parallelized the counting.

for i in 1:100
    resutls = realLineCounts(i/100, M);
    println(results)
end