/* rho2.c */
/* Rene Peralta: 6/1994 */
/* REVISION 7/2007 : Thanks to Andrew Sutherland for noticing errors in the estimation
   of rho2(u,v) when v is close to 1. The problem was that the coefficient c[0,1] was 
   not initialized. The computation of rho2 in the range of inputs of interest to 
   optimization of factoring methods IS NOT AFFECTED. */

/* the method used here is described in the paper
   "Asymptotic Semi-Smoothness Probabilities" - to appear
   in Math Comp. - by Eric Bach and Rene Peralta. Be aware
   that this is an iterative version of the (recursive)
   computation described in that paper */
	    
/* set ustart,uend,vstart,vend,delta ; compile cc rho2.c -lm */

/* you can also set parameters num_terms and rho_precision */
/* CAREFUL: if you do set these parameters then make sure num_terms <= rho_precision */

/* num_terms is the number of terms used in the summation for the
   two-dimensional rho function. rho_precision is the number of
   terms used in the summation for the rho function */

/* OUTPUT IS FOR TABULAR LATEX FORMATTING */

/* the range of u is [ustart,uend] */
/* the range of v is [vstart,vend] */
/* the step size for both u and v is delta */

/* no efforts have been made to optimize code (it is quite fast already) */

/* this will not handle large values of u and v */
/* e.g. you will cause a floating exception if you try to compute rho2[200,100] */
/* this shouldn't be hard to fix but hasn't been done as of 6/22/94 */


/* problems to peralta@cs.uwm.edu */


#include <stdio.h>
#include <math.h>
#define  delta 1.0
#define  ustart 1.0
#define  uend 20.0
#define range (int)(uend + 1) /* integer upper bound on value of u */
#define  vstart 1.0
#define  vend 10.0
#define num_terms 22
#define rho_precision 55

double c[rho_precision][range];

double power();
double rho();
double J();
double rho2();

double round();

main()
{
  double u,v;
  int i;
  
  coefficients(); /* this routine computes the c[i][j]'s */
  
  /*  This produces a table of rho values */
  u = ustart;
  v = vstart;
  printf("    &");
  while( v <= vend) 
  {
    if (v < vend ) printf("{\\bf %2.1lf }            &", v);
    else printf("{\\bf %2.1lf }           ", v);
    v = v + delta;
  };
  printf("\\\\ \\hline \n");

  while (u <= uend)
  {
    v = vstart;
    printf("{\\bf %3.1lf } & " , u);
    while ( v <= vend )
    {
      if (v < vend )
      {
	if ( (v <= vend) && (v <= u) ) printf("%le  & ", rho2(u,v));
        else printf("    -----     & ");
      }
      else
      {
        if ( (v <= vend) && (v <= u) ) printf("%le   ", rho2(u,v));
        else printf("    -----      ");
      }
      
      v = v + delta;
      v = round(v*10)/10.0;
    }
    printf("\\\\ \\hline \n");
    u = u + delta;
    u = round(u*10)/10.0;
  }
}

/* the 2-dimensional rho function */

double rho2(u,v)
double u,v;
{
  if (v > u) { printf("rho2: v should be <= u \n"); exit(0); }
  return(J(u,v) + rho(u));
}

/* calculates int_v^u t^(-1) rho(w - w/t) dt */

double J(u,v)
double u,v;
{
  double A,B,C,sum,res,luv,acc;
  double A_C,B_C;
  long k;
  double H[num_terms];
  int i,j;
  double u_i, v_i;  /* values of u,v at ith partition */

  k = ceil( u - 1 );
  C = k - u;
  A = u/v + C;
  B = 1 + C;
  A_C = A/C;
  B_C = B/C;

  
  /* If v >= u / (1 - C)  there is no partitioning */
  if ( v >=  u / (1 - C) )
  {
      luv = log(u/v);
      /* calculate H[0] */
      H[0] = luv;
      
      /* calculate H[i] */
      if ( C > -0.2 )  /* danger of overflow , so bring C^i into summations */
      for ( i = 1; i < num_terms; i++ )
      {
	sum = power(C,i) * luv;
	for ( j = 1; j <= i; j++)
	  sum = sum + power(A,j)*power(C,i-j)/ j - power(B,j)*power(C,i-j) / j ;
        
	H[i] = sum;
      }
      else /* C <= -0.2 */
      for ( i = 1; i < num_terms; i++ )
      {
	sum = luv;
	for ( j = 1; j <= i; j++)
	  sum = sum + power(A_C,j)/j - power(B_C,j) / j ;
        
	H[i] = sum * power(C,i);;
      };

      res = 0;
      for (i = num_terms - 1; i > -1 ; i--) res = res + c[i][k] * H[i];
      return(res);
   }  /* no partitioning case */
  
   
   /* case where there is partitioning */
   
   acc = 0; /* I will sum to acc the contribution of each piece */
   
      /* first piece */
      
      /* set variables for first piece */
      k = ceil(u-1);
      u_i = u;
      v_i = u/(u - k + 1);
      C = k - u;
      A = 1;
      B = 1 + k - u;
      A_C = 1/C;
      B_C = 1 + 1/C;
      luv = log(u/v_i);
      
      /* calculate H[0] */
      H[0] = luv;
      
      /* calculate H[i] */
      if ( C > -0.2 )  /* danger of overflow , so bring C^i into summations */
      for ( i = 1; i < num_terms; i++ )
      {
	sum = power(C,i) * luv;
	for ( j = 1; j <= i; j++)
	  sum = sum + power(A,j)*power(C,i-j) / j - power(B,j)*power(C,i-j) / j ;
        
	H[i] = sum;
      }
      else /* C <= -0.2 */
      for ( i = 1; i < num_terms; i++ )
      {
	sum = luv;
	for ( j = 1; j <= i; j++)
	  sum = sum + power(A_C,j)/j - power(B_C,j) / j ;
        
	H[i] = sum * power(C,i);;
      };

      res = 0;
      for (i = num_terms - 1; i > -1 ; i--) res = res + c[i][k] * H[i];
      acc = res;
      /* intermediate pieces */
      /* set variables for intermediate piece */
      while (1)
      {
	u_i = v_i;
        k = k - 1;
	C = k - u;
	if ( v >= u/(1-C) ) break; /* go to last level */
	v_i = u/(1 - C);
	luv = log(u_i/v_i); 
        /* A = 1; B = 0 */
        
	/* the previous test should be equivalent to v_i < v */

        /* calculate H[0] */
        H[0] = luv;
        
	/* calculate H[i] */
        for ( i = 1; i < num_terms; i++ )
        {
	  sum = 0;
	  for ( j = num_terms ; j > i ; j--) sum = sum - ((double) 1.0)/(power(C,j-i)*j);
	  H[i] = sum;
	}
        res = 0;
        for (i = num_terms - 1; i > -1 ; i--) res = res + c[i][k] * H[i];
        acc = acc + res;
      };

      /* last piece*/
      /* set variables for last piece */
      /* u_i and k and C come from above already */
      v_i = v;
      A = u/v + C;
      luv = log(u_i/v_i); 
      /*  B = 0 */
      
      /* calculate H[0] */
      H[0] = luv;
        
      /* calculate H[i] */
      for ( i = 1; i < num_terms; i++ )
      {
	  sum = 0;
	  for ( j = i+1; j <= num_terms; j++) sum = sum - power(A,j)/(power(C,j-i)*j);
	  H[i] = sum;
       }
       res = 0;
       for (i = num_terms - 1; i > -1 ; i--) res = res + c[i][k] * H[i];
       acc = acc + res;
       return(acc);
}


double round(r)
double r;
{
  return((double) (int) (r + 0.5) );
}

double rho(x)
double x;
{ int i,k;
  double u, sum;
  
  if (x <= 1.0) return(1.0);

  k = ceil(x);
  u = k - x;

  sum = 0;
  for (i = rho_precision - 1 ; i > 0 ; i--) 
    { sum = sum + c[i][k]*power(u,i);
    } 
  sum = sum + c[0][k];
  
  return(sum);
}

coefficients()
{
  int i,j,k;
  
  c[0][1] = 1.0;  /* FIX 7/2007 */
  c[0][2] = 1 - log(2.0);

  for (i = 1; i < rho_precision; i++) c[i][2] = 1.0/ (i * power(2.0,i));
  
  for (k = 3; k < range; k++)
  {
    for (i = 1; i < rho_precision; i++)
    { 
       c[i][k] = 0;
       
       for (j = 0; j <= (i - 1) ; j++) 
          c[i][k] = c[i][k] + c[j][k-1]* power((double)k,j);
       
       c[i][k] = c[i][k]/(i*power((double)k,i));
    }

    c[0][k] = 0;
    for (j = 1; j< rho_precision ; j++) c[0][k] = c[0][k] + c[j][k]/(j+1) ;
    c[0][k] = c[0][k]/ (k-1);
  }
}

double power(x,i)
double x;
int i;
{ 
  /* safety */
  if (i < 0) 
  { printf("negative exponent in power function\n");
    exit(0);
  }
  
  if (i == 0) 
  { 
    /* safety */ 
    if (x == 0) { printf("0^0 \n"); exit(0);}
    else return(1.0);
  }
  
  if (x == 0) return((double) 0);
  if (i == 1) return(x);
  if (x > 0) return(exp(i*log(x)) );
  if ( (x < 0) && ((i % 2) == 0) ) return( exp(i*log(-x)) );
  if ( (x < 0) && ((i % 2) == 1) ) return( - exp(i*log(-x)) );
  
} /* power */