# Returns the radius of an (n-1)-sphere defined by n+1 points in R^n
# L = set of n+1 points, each an n-dimensional vector
def sphere_radius(L):
    n = len(L)-1
    matL = [[0]*(n+2)]
    for pt in L:
        tempL = [vector(pt)*vector(pt)]
        for pos in range(n):
            tempL.append(pt[pos])
        tempL.append(1)
        matL.append(tempL)
    M = matrix(CDF,n+2,n+2,matL)
    return sqrt(reduce(lambda x,y: x+y, map(lambda z: minor(M,0,z-1)**2/(4*minor(M,0,0)**2),range(1,n+2)))+(-1)**n*minor(M,0,n+1)/minor(M,0,0))

# Returns the radius of a (k-1)-sphere defined by k+1 points in R^n
# L = set of k+1 points, each an n-dimensional vector
def sphere_radius_general(L):
    k = len(L)-1
    n = len(L[0])

    # Shift vectors so one is at origin
    L1 = []
    for vec in L[1:]:
        L1.append(vec-L[0])

    # Complete basis of R^n
    L2 = []
    for ind in range(n-k):
        L2.append([0]*ind+[1]+[0]*(n-ind-1))
    M = matrix(CDF,n,n,L1+L2)
    Q,R = M.transpose().QR()

    # Make vectors for passing to other function
    L3 = [vector(CDF,[0]*k)]
    Qinv = Q.inverse()
    for vec in L1:
        L3.append((Qinv*vec)[:k-n])
    return sphere_radius(L3)

# Returns the (i,j)-minor (determinant when ith row, jth col removed) of input matrix mat
def minor(mat,i,j):
    return mat.delete_rows([i]).delete_columns([j]).det()

# Test
v1 = vector(CDF,[3,2,1])
v2 = vector(CDF,[0,-1,3.3])
v3 = vector(CDF,[5,6.8,-9.1])
A = matrix(CDF,3,3,[v2-v1,v3-v1,[1,0,0]])
Q,R=A.transpose().QR()
a1 = (Q.inverse()*(v2-v1))[:-1]
a2 = (Q.inverse()*(v3-v1))[:-1]
print sphere_radius([vector(CDF,[0,0]),a1,a2])
print sphere_radius_general([v1,v2,v3])