From sources-request at octave dot org Mon Nov  1 11:30:00 2004
Subject: unit simplex
From: Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov>
To: sources at octave dot org
Date: Mon, 1 Nov 2004 12:20:07 -0500

Just in case anybody needs n+1 points equally distributed
on an n-dimensional sphere.

- Paul

% V = unitsimplex(n)
%
% Return a regular n-dimensional simplex with points on
% the n-sphere of radius 1.

% This program is public domain

function z = unitsimplex(n)
  % Find a regular simplex with a vertex at the origin
  % and unit length sides.  Use a recursive algorithm
  % which adds vertices one at a time, starting from
  % the center of the simplex in the previous dimension
  % and choosing the coordinate in the next dimension
  % such the the distance to the origin is unit length.
  % By symmetry, since we started from the center of
  % the simplex the distance to all other vertices will
  % also be unit length.
  z = zeros(n+1,n);
  z(1,1) = 1;   % the first vertex is at (1,0,...,0)
  for i=2:n
    % The center of the next higher dimensional simplex
    % is above the center of the last simplex, but with
    % the final coordinate as the average of all previous 
    % coordinates, which is easy to calculate since all but
    % but one of these coordinates is 0.  Since z(i,:) = 0,
    % the following is equivalent to:
    %    z(i,:) = mean(z(1:i,:))
    z(i,:) = z(i-1,:);
    z(i,i-1) /= i; 
    z(i,i) = sqrt(1-sumsq(z(i,:)));
  end
  % the final vertex is at (0,0,...,0)

  % Find the center of the simplex and subtract it.
  % Use the same algorithm we use above for the center.
  c = z(n,:);
  c(n) /= n+1;
  for i=1:n+1, z(i,:) -= c; end

  % Scale the simplex into a sphere of unit radius.  Since
  % the last vertex we added will only have one coordinate
  % away from the origin (all other coordinates were at the
  % center of the simplex), make that coordinate have unit
  % length.
  z /= z(n,n);
  

%!# test that all distances are the same
%!test
%! z = unitsimplex(5);
%! d = zeros(rows(z));
%! for i=1:rows(z)
%!   for j=i+1:rows(z)
%!     d(i,j) = d(j,i) = sqrt(sumsq(z(i,:)-z(j,:)));
%!   end
%! end
%! assert(d,d(2,1)*(1-eye(rows(d))),100*eps)

%!# test that we are sitting in a unit sphere
%!assert(sumsq(unitsimplex(5)'),ones(1,6),100*eps);

%!# test that we are centered in the sphere
%!assert(mean(unitsimplex(5)),zeros(1,5),100*eps);