// MagicSeries.cpp : console application : Walter Trump 2005-01-17

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define min(a, b)  (((a) < (b)) ? (a) : (b)) 
#define max(a, b)  (((a) > (b)) ? (a) : (b)) 

int main()
{
  const int nMax=25;					// order n (Maximum)
  const int uMax=nMax*nMax;				// integer u (Maximum)
  const int sMax=(uMax+1)*nMax/2;		// sum s (Maximum)
  int Mag[nMax+1], sMn[nMax+1], sMx[nMax+1];
  double (*Num)[uMax+1];		
  Num = new double [sMax+1][uMax+1];	// main array

  printf(" Number of magic series for squares \n\n");
  double t=(double)clock()/CLK_TCK;
  
  Mag[nMax]=sMax;
  sMx[nMax]=sMax;
  sMn[nMax]=sMax;
  for(int n=nMax-1; n>=1; n--)			// magic constants and borders for loops
  {
	  Mag[n]=(n*(n*n+1))/2;
	  sMx[n]=(n*(2*sMx[n+1]-n+1))/(2*n+2);
	  sMn[n]=max((n*(n+1))/2 , sMn[n+1]-uMax);
	  sMn[n]=min(sMn[n] , Mag[n]);
  }

  for(int s=sMx[1]; s>=1; s--)			// basis of iteration
  {
    for(int u=s; u<=uMax; u++)
      Num[s][u]= 1;
  }

  for(n=2; n<=nMax; n++)				// iteration
  {
    int Sn=(n*(n-1))/2;
	for(int s=sMx[n]; s>=sMn[n]; s--)
    {
	  double Sum=0;
	  int uMn=max(n , 1+(s+Sn-1)/n);
	  int uMx=min(s-Sn, uMax);
	  for(int u=uMn; u<=uMx; u++)
	  {
        Sum += Num[s-u][u-1];
		Num[s][u]= Sum;
	  }
	  for(u=uMx+1; u<=uMax; u++)
		Num[s][u]= Sum;
    }
	printf(" Order %2.0i:  %1.14E \n", n, Num[Mag[n]][n*n]);
  }

  t=((double)clock()/CLK_TCK)-t;
  printf("\n Time: %.3f ms \n",t);

  delete [] Num;
  printf("\n Press Enter !  ");
  getchar();
  return 0;
}