function [ intval, errval, intcls ] = sphrlr( s, d, f, n )
%
%   [ intval, errval, intcls ] = sphrlr( s, d, f, n )
%  SPHRLR computes a randomized fully symmetric rule approximation 
%   to an integral over the surface of the unit sphere: norm(x) = 1. 
%**************  parameters for sphrlr  ********************************
%***** input parameters for sphrlr( s, d, f, n )
%  s       integer number of variables
%  d       integer degree parameter
%  f       user defined real function integrand f(x).
%  n       number of random rotations for the integration rule
%****** output parameters
%  intval  intval will be an approximation of polynomial degree 2d + 1.
%  errval  error estimate for intval: 3 x standard error
%  intcls  integer total number of f values needed for intval
%****** example use
% >> ft = inline('x(1)^2*x(2)^4*cos(x(4))*exp(-x(3))','x');
% >> for d = 2 : 6, [ iv ev pts ] = sphrlr(4,d,ft,8); disp([d iv ev pts]), end
% 
%     Sphrlr uses a transformation of Sylvester's interpolation formula
%     for simplices, as described in "Fully Symmetric Interpolatory
%     Rules for Multiple Integrals over Hyper-Spherical Surfaces",
%       J. Comp. Appl. Math. 157 (2003), pp. 187-195
%   Software and paper written by 
%     Alan Genz
%     Department of Mathematics
%     Washington State University
%     Pullman, Washington 99164-3113  USA
%     email: alangenz@wsu.edu
%
%***********************************************************************
%
%***  begin loop for each d
%      for each d find all distinct partitions m with mod(m) = d
%
  intval = 0; errval = 0; intcls = 0;
  for i = 1 : n
    q = ranrf(s); intvli = 0; m = 0*[1:s]; m(1) = d; 
    while m(1) <= d 
      %     
      %***  calculate weight for partition m and fully symmetric sums
      %
      wt = weight( s, m, d ); 
      if abs(wt) > 1d-10, [ fs sumcls ] = fulsmr( s, m, d, f, q ); 
	intvli = intvli + wt*fs; intcls = intcls + sumcls;
      end
      m = nxprtr( s, m );
    end
    df = ( intvli - intval )/i; 
    intval = intval + df; errval = ( i - 2 )*errval/i + df^2; 
  end
  srf = srfsph(s); errval = 3*srf*sqrt(errval); intval = srf*intval;
%
%  end function sphrlr
%
function v = srfsph( s )
%
% compute surface content for s-sphere
%
   v = 2; if mod(s,2) == 0, v = v*pi; end
   for i = s-2 : -2 : 1, v = 2*v*pi/i; end
%
% end function srfsph
%
function q = ranrf( n )
%
% RANRF(n) generates random nxn orthogonal matrix as product of reflectors 
%
  q = eye(n); x = zeros(n,1);
  for j = n - 1 : -1 : 1, x(j:n) = randn(n-j+1,1); s = sign(x(j));
    % determine reflector and apply reflector to q
    sg = s*norm(x); x(j) = x(j) + sg; gm = 1/( sg*x(j) );  
    if s < 0, q(j:n,j:n) = q(j:n,j:n) - x(j:n)*( gm*x(j:n)'*q(j:n,j:n) );
    else, q(j:n,j:n) =  x(j:n)*( gm*x(j:n)'*q(j:n,j:n) ) - q(j:n,j:n); end
  end
%
% end function ranrf
%
function ws = weight( s, m, d )
%
%***  calculate weight for partition m
%
  ws = 0; k = 0*[1:d];
  while 1
    is = 0; ks = 0; mc = 1; kn = 1; pd = 1;
    for i = 1 : d, is = is + 1;
      if k(i) == 0, pd = pd*( 1 - is );
      else, pd = d*kn*pd/( s + ks ); ks = ks + 2; kn = kn + 2;
      end 
      if is == m(mc), mc = mc + 1; is = 0; kn = 1; end
    end
    ws = ws + pd;
    for i = 1 : d
      k(i) = k(i) + 1; if k(i) < 2, break, end, k(i) = 0;
    end
    if i == d & k(d) == 0, break, end
  end
  for i = 1 : s, for j = 1 : m(i), ws = ws/( m(i) - j + 1 ); end, end
%
function [ fs, sumcls ] = fulsmr( s, m, d, f, q )
%
%***  compute fully symmetric basic rule sum
%
  sumcls = 0; fs = 0; fin = 0;
  %
  %*******  compute centrally symmetric sum for permutation m
  %
  while fin == 0, x = -sqrt( m/d ); 
    [ sym, symcls ] = symsmr( s, x, f, q );
    fs = fs + sym; sumcls = sumcls + symcls;
    %     
    %*******  find next distinct permutation of m and loop back
    %          to compute next centrally symmetric sum
    %
    [ m, fin ] = nxperm( s, m );
  end
%     
%*****end loop for permutations of m and associated sums
%
%  end function fulsmr
%
function [ sym, sumcls ] = symsmr( s, zi, f, q )
%
%***  function to compute symmetric basic rule sum
%
  sumcls = 0; sym = 0; z = zi;
  %
  %*******  compute centrally symmetric sum for zi
  %
  while 1
    sumcls = sumcls + 1; sym = sym + ( feval( f, z*q' ) - sym )/sumcls;
    for i = 1 : s, z(i) = -z(i); if z(i) > 0, break, end, end
    if z == zi, break, end
  end
%
%  end function symsmr
%
%
function [ m, fin ] = nxperm( s, mi )
%
%***  compute the next permutation of mi
%
  m = mi; fin = 0;
  for i = 2 : s
    if m(i-1) > m(i), mi = m(i); ixchng = i - 1;
      if i > 2
	for l = 1 : ixchng/2
	  ml = m(l); imnusl = i - l; m(l) = m(imnusl); m(imnusl) = ml;
          if ml  <= mi, ixchng = ixchng - 1; end
	  if m(l) > mi, lxchng = l; end
	end
	if m(ixchng) <= mi, ixchng = lxchng; end
      end 
      m(i) = m(ixchng); m(ixchng) = mi; return
    end 
  end
  %
  %*** restore original order to m.
  %
  for i = 1 : s/2, mi = m(i); m(i) = m(s-i+1); m(s-i+1) = mi; end
  fin = 1;
%
%  end function nxperm
%
function p = nxprtr( s, pi )
%
%***  compute the next s partition of pi
%
  p = pi; psum = p(1);
  for i = 2 : s, psum = psum + p(i);
    if p(1) <= p(i) + 1, p(i) = 0;
    else
      p(1) = psum - (i-1)*( p(i) + 1 );
      for l = 2 : i, p(l) = p(i) + 1; end
      return
    end
  end
  p(1) = psum + 1;
%
%  end function nxtptr
%