/* MAXIMA SOURCE CODE
   Author: Assad Ebrahim, (assad.ebrahim@mathscitech.org)
   License: LGPL (Lesser GPL)
   Copyright (c) 2010, Assad Ebrahim.
   
   Matrix Solution of the Finite Sum of p-th Powers of Integers
   using Equation (13) 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 wxMaxima, run the following commands:
   
   load("c:\\tmp\\sumkp_matrix.mac"); [replace with path to downloaded file]
   SpN_mat(p);                                [positive integer p: 1,2,...]

   NOTES:  - SpN_mat(p) returns the closed form formula for the
           finite summation of p-th powers of integers in a variety of formats: 
           coefficients, polynomial (series), polynomial (factored), prettified
   */


/* **********************
   PUBLIC Function: Solves the p+1 by p+1 linear system (*):  
   Input: positive integer p: 1,2,... 
   Return: closed form formula for S_p(N) */

   SpN_mat(p):= block([a, eq],     /* give subroutine variables local scope */

      v : makelist(a[i], i, 0, p), /* create list of unknowns (0-indexed)   */
      
                                   /* create list of equations (0-indexed)  */
     eq : makelist(sum(binom(j+1,i)*a[j],j,i,p) = binom(p,i), i, 0, p),  

     /* find coefficients of solution polynomial by solving linear system */
          print("Linear System of Equations:"),
		  print(eq),
     sol : linsolve(eq, v),
     
     /* create polynomial: inner product of {N^i} with coefficients */
          print("Solution (Coefficient Values):"),
		  print(reverse(sol)),
     pol : makelist(N^(i+1),i,0,p), /* monomials {N^i} */
     
     /* closed form formula */
     cff : rhs(sol.pol),  /* inner product of {N^i} with coefficients */
     
     /* prettified outputs */
          print("Polynomial Solutions:"),
          print("(a) series:", cff),
          print("(b) over common denominator:", ratsimp(cff)),
          print("(c) factors, irreducible over Q:", factor(cff)),
    
     cff  /* return closed form formula  */
  )$

/* ***********************
   Package complete */
print("Loading completed successfully.");
print("Commands: SpN_mat(p)   [p a positive integer: 1,2,3,...]");


/*
   Revision History:
   2/28/2010 - 002 - AKE - revised layout
   2/22/2010 - 001 - AKE - initial  */