# Import the necessary modules/libraries import numpy as np from numpy import poly1d from scipy import linalg import math import random import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # compute j of a vector in C2 def quatj(u): if u.size!=2: print("Error, the vector is not in C2") return if linalg.norm(u)==0: print("Error, the vector in C2 is 0") return return np.array([-u[1].conjugate(),u[0].conjugate()]) # Compute the stereographic projection of a non-zero vector in R3 def stereo(x): if x.size!=3: print("Error, the vector is not in R3") return if linalg.norm(x)==0: print("Error, the vector in R3 is 0") return r=linalg.norm(x) if x[2]!=r: return np.array([1,complex(x[0],x[1])/(r-x[2])]) elif x[2]==r: return np.array([0,-1]) # Convert a non-zero vector in C2 into a degree less than or equal to 1 polynomial of one variable t def poly(u): if u.size!=2: print("Error, the vector is not in C2") return if linalg.norm(u)==0: print("Error, the vector in C2 is 0") return return poly1d(np.array([u[0],-u[1]]), variable="t") # Compute the 2 by 2 determinant of two vectors in u and v in C2 def det(u,v): mytemp = np.zeros((2,2),dtype=np.complex) mytemp[0,:] = u[:] mytemp[1,:] = v[:] return linalg.det(mytemp) # given an np array with the n points in R3 as rows, compute the associated determinant def MyDet(x): n=x.shape[0] if x.shape[1]!=3: print("Error, the number of columns of x must be 3") return s=np.zeros((n,2), dtype=np.complex) for i in range(n): s[i,:] = stereo(-x[i,:]) neg_s=np.zeros((n,2),dtype=np.complex) for i in range(n): neg_s[i,:] = quatj(s[i,:]) p=np.zeros((n,n,2),dtype=np.complex) for i in range(n-1): for j in range(i+1,n): p[i,j,:]=stereo(-x[i,:]+x[j,:]) p[j,i,:]=quatj(p[i,j,:]) q=np.zeros((n,n,2),dtype=np.complex) for i in range(n-1): for j in range(i+1,n): q[i,j,:]=stereo(-x[i,:]-x[j,:]) q[j,i,:]=quatj(q[i,j,:]) M=np.zeros((2*n,2*n),dtype=np.complex) for i in range(n): mypoly = poly1d([complex(1,0)],variable="t") for j in range(n): if (j!=i): mypoly=mypoly*poly(p[i,j,:]) for j in range(n): if (ij): mypoly=mypoly*poly(q[j,i,:]) mypoly=mypoly*poly(s[i,:]) mycoeffs = mypoly.coeffs[::-1] M[2*i,0:mycoeffs.size]=mycoeffs[:] for i in range(n): mypoly = poly1d([complex(1,0)],variable="t") for j in range(n): if (j!=i): mypoly=mypoly*poly(p[j,i,:]) for j in range(n): if (ij): mypoly=mypoly*poly(q[i,j,:]) mypoly=mypoly*poly(neg_s[i,:]) mycoeffs = mypoly.coeffs[::-1] M[2*i+1,0:mycoeffs.size]=mycoeffs[:] P = 1 for i in range(n): P*=det(s[i],neg_s[i]) for i in range(n-1): for j in range(i+1,n): P*=(det(p[i,j],p[j,i])*det(q[i,j],q[j,i]))**2 return linalg.det(M)/P # returns a pseudo-random number between -a and a def rand(a): return 2*a*random.random()-a # generate N configurations of n points having each coordinate # between -a and +a, put these N configurations in an x with 3 # indices, such that the distances between 0 and x_i are all greater or equal to a*eps # and the distances between x_i and x_j, and those between x_i and -x_j are all # greater or equal to a*eps. compute the corresponding N determinants and put them in # an np array called MyD # return x and MyD. # eps should be small enough so that there are configurations with those requirements, but not # too small, otherwise one would get rounding errors # as a rule of thumb, take a to be approximately between 1.5*n and 2*n, eps to be 0.1 for N <= 30 def test(n,N,a, eps): x=np.zeros((N,n,3),dtype=np.float) for k in range(N): for i in range(n): condition = True while condition: for j in range(3): x[k,i,j]=rand(a) condition = (linalg.norm(x[k,i,:])