%   OCTAVE SOURCE CODE
%   Author: Carol Ouellette, (carol.ouellette@mathscitech.org)
%   License: LGPL (Lesser GPL)
%   Copyright (c) 2010, Carol Ouellette.

%   Matrix solution of the finite summation of integer powers
%   using Equation (12) from the paper: 
%   "Finite Summation of Integer Powers, Part 3", 
%   Assad Ebrahim & Carol Ouellette, Feb 2010
%   Paper at: http://mathscitech.org/papers/ebrahim-sum-powers-3.pdf
%   Source Code at: http://mathscitech.org/articles/downloads#source

% =====================================
%  USAGE INSTRUCTIONS: from within Octave, run the following commands:
%   addpath("c:\\tmp")   % or wherever this source file is
%   sumkp_matrix(p)     % where p is a positive integer
%        - this returns coefficients for the closed form formula;
%	- high order coefficients listed first, i.e. a_p+1 at left, down to a_1 at right
% ============================================

% Function listing coefficients from a_p+1 to a_1 as fractions

function c = sumkp_matrix(p)
	M=zeros(p+1); 		% initialize matrix M
	for i=1:p+1
		for j=i:p+1		% set upper triangular elements
			M(i,j)=nchoosek(j,i-1);  %  Equation (12)
		end;
	end;
	% show
	disp("M:");
	disp(M);

	b=zeros(p+1,1); % initialize column matrix b
	for i=1:p+1
		b(i)=nchoosek(p,i-1);  %  Equation (12)
	end
	% show
	disp("b:");
	disp(b);  % show
	
	c=inv(M)*b;  % Solve to obtain coefficients

	c=fliplr(c');		% List from highest index to lowest
	sml=abs(c)<1e-8;	% Numerical correction: set tiny coeffients to exactly zero
	c(find(sml==1))=0;
	
	disp("Solution Coefficients (High to Low order):")
	c=rats(c);  		% return coefficients as fractions

end;


%   Revision History:
%   4/04/2010 - 004 - AKE - revised to match derivation in Finite Summations paper
%   3/08/2010 - 003 - CJO - extension
%   2/04/2010 - 002 - CJO - obtaining rational coefficients
%   1/29/2010 - 001 - CJO - initial exploration */