Unit simplex in N dimension, as a triangular matrix. 1st vertex is the origin.

unit_simplex.py

import numpy

'''
Unit simplex in N dimension, as a triangular matrix,
only positive coeffs.
1st vertex is the origin, not in the returned matrix

I don't know if this algorithm have been discovered
before. I find it elegant, and I believe it can be
very accurate and stable with arbitrary precision numbers
'''
def get_unit_simplex(N, dtype = float):
    S = numpy.zeros((N, N), dtype = dtype)
    
    # Iterative construction
    S[0, 0] = 1
    for i in range(1, N):
        S[i, :i - 1] = S[i - 1, :i - 1]
        S[i, i - 1] = S[i - 1, i - 1] / ((i + 1) ** 2)
        S[i, i] = 1 - S[i, :].sum()
        
    # Job done
    return numpy.sqrt(S)