// Optimal randomized ski rental - B.P. Weems, 1/17/2004
// See S. Phillips & J. Westbrook, "On-line Algorithms: Competitive
//     Analysis and Beyond", Jan 24, 1998.
// Solved using system of linear equations & Householder method
// See P. Brinch-Hansen, "Householder Reduction of Linear Equations"
//     ACM Computing Surveys 24:2, June 1992, 185-194.
#include <stdio.h>
#include <math.h>

// Householder method - notice that matrix is stored transposed

#define SIZE 2000

typedef float column[SIZE+1];
typedef column matrix[SIZE+1];

int n;
matrix a;
column b;

float product(i,a,b)
int i;
column a,b;
{
float ab;
int k;

ab=0.0;
for (k=i;k<=n;k++)
  ab+=a[k]*b[k];
return ab;
}

eliminate(i,ai,vi)
int i;
column ai,vi;
{
float anorm,dii,fi,wii;
int k;

anorm=sqrt(product(i,ai,ai));
dii=(ai[i]>0.0) ? (-anorm) : anorm;
wii=ai[i]-dii;
fi=sqrt(-2.0*wii*dii);
vi[i]=wii/fi;
ai[i]=dii;
for (k=i+1;k<=n;k++)
{
  vi[k]=ai[k]/fi;
  ai[k]=0.0;
}
}

transform(i,aj,vi)
int i;
column aj,vi;
{
float fi;
int k;

fi=2.0*product(i,vi,aj);
for (k=i;k<=n;k++)
  aj[k]-=fi*vi[k];
}

reduce(a,b)
matrix a;
column b;
{
column vi;
int i,j;

for (i=1;i<=n-1;i++)
{
  eliminate(i,a[i],vi);
  for (j=i+1;j<=n;j++)
    transform(i,a[j],vi);
  transform(i,b,vi);
}
}

main()
{
int i,j,c;
double sum,cdouble,alpha,pi,contrib;
column x;

printf("Enter c (cost to buy, cost to rent is $1)\n");
scanf("%d",&c);
cdouble=c;

// Compute results using closed form
alpha=1.0/(pow(1.0+1.0/(cdouble-1.0),cdouble)-1.0) + 1.0;
printf("alpha %f\n",alpha);

pi=(alpha-1.0)/cdouble;
for (j=1;j<=c;j++)
{
  pi*=cdouble/(cdouble-1.0);
  printf("pi[%d] %f\n",j,pi);
}

for (i=1;i<=c;i++)
{
  sum=0.0;
  pi=(alpha-1.0)/cdouble;
  for (j=1;j<=c;j++)
  {
    pi*=cdouble/(cdouble-1.0);
    contrib=pi*(i<j ? i : j-1+cdouble);
//    printf("contrib[%d] %f\n",j,contrib);
    sum+=contrib;
  }
//  printf("%d has expected c.r. %f\n",i,sum/i);
}

// Householder to get alpha and pi values 
n=c;
for (i=1;i<=c;i++)
{
  for (j=1;j<=c;j++)
  {
    contrib=i<j ? i : j-1+cdouble;
    a[j][i]=contrib; // matrix is store transposed
//    printf("%f ",contrib);
  }
  b[i]=i;
//  printf("%f\n",b[i]);
}
printf("Performing Householder reduction\n");
reduce(a,b);

printf("Extracting solutions:\n");
x[n]=b[n]/a[n][n];
for (i=n-1;i>=1;i--)
{
  sum=0.0;
  for (j=i+1;j<=n;j++)
    sum+=a[j][i]*x[j];
  x[i]=(b[i]-sum)/a[i][i];
}
/*
for (i=1;i<=n;i++)
  printf("%f\n",x[i]);
*/
// Scale so that the probabilities add to 1.  Also provides alpha.
sum=0.0;
for (i=1;i<=n;i++)
  sum+=x[i];
printf("Solution vector sums to %f\n",sum);
alpha=1.0/sum;
printf("alpha from Householder is %f\n",alpha);
for (i=1;i<=n;i++)
  printf("pi[%d] %f\n",i,x[i]*alpha);
}