# Julia Code for sampling real cubics

[9]:

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


## Computing number of real lines

[13]:

## 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