/*********************************************************************************************************************************************************************/
float assign_value(ndad, nkid,i,j,l,p,q)
list ndad, nkid;
int i, j, l; 
float p[][8][8],q[][8][8]; {
  int k, n, m, biggest, rndmax(), round[8];
  float sum=0.0, tmp, cube=0.0, delta[8],
        dummy, dummy1, dummy2, pi[8];

  dummy2=0.0;
  for(m=0;m<8;m++) {
    dummy=0.0;
    for(n=0;n<8;n++)
      dummy+=(p[i][n][m]*(ndad->value[n]));
    dummy2+=(q[j][m][l]*dummy);
  }

  if(dummy2==0.0) {
    ndad->denom[i][j][l]=0.0;
    return (0.0);
  }
  else {
    ndad->denom[i][j][l]=dummy2;
    for(k=0;k<8;k++) {
      dummy1=0.0;
      for(n=0;n<8;n++)
        dummy1+=p[i][n][k]*ndad->value[n];
      tmp=(q[j][k][l]*dummy1/dummy2);                                                /*  Everything above this point in the procedure is calculation of the          */
      pi[k]=((float)(rint(tmp/rnd))*rnd);                                            /*  denom field and the pi vector.  Here, pi[k] is being 'rounded' to the       */
      delta[k]=fabs(pi[k]-tmp);                                                      /*  multiple of rnd. delta[k] holds the absolute value of the difference of tmp */
      if (pi[k]>tmp)                                                                 /*  and pi[k], for future calculation. Round[k] holds 1 if tmp was rounded 'up' */
        round[k]=1;                                                                  /*  to pi[k] and 0 if it was rounded 'down'.  Sum holds the sum of the pi       */
      else                                                                           /*  vector's elements.                                                          */
        round[k]=0;                                                                  /*                                                                              */
      sum+=pi[k];                                                                    /*                                                                              */
    }                                                                                /*                                                                              */
    if (sum>(1.0+rnd-0.001)) {                                                       /*  If sum is not equal to one, then it is not a good pi vector, and must be    */
      while (sum>1.0) {                                                              /*  changed to the closest pi vector.  Closest here means the pi vector that    */
        biggest=rndmax(delta,round,1);                                               /*  minimizes the absolute value of the difference between the values of tmp    */
        pi[biggest]-=rnd;                                                            /*  and a pi vector of multiples of rnd.  This is done by finding the largest   */
        delta[biggest]=0.0;                                                          /*  rounding error that is in the same  direction of sum's error, i.e. if       */
        sum-=rnd;                                                                    /*  sum is equal to 1.05 it is optimal to change the pi value that was changed  */
      }                                                                              /*  the most to be rounded up (like 0.426 and not .449 when rnd =0.05).         */
    }                                                                                /*  Rndmax() does this, taking round and delta as arguments.  The while loops   */
    if (sum<(1.0-rnd+0.001)) {                                                       /*  under the if statements assure that if sum is equal to 1.10 or 1.15, this   */
      while (sum<1.0) {                                                              /*  large of an error will will be compensated for.                             */
        biggest=rndmax(delta,round,0);
        pi[biggest]+=rnd;
        delta[biggest]=0.0;
        sum+=rnd;
      }
    }
    j=1;
    for(i=0;i<8;i++) {                                                               /*  This is the calculation of cube, the signature value.  It is fairly         */
      cube+=(j*j*j)*pi[i]*pi[i]*pi[i];                                               /*  straight-forward.                                                           */
      nkid->value[i]=pi[i];
      j++;
    }
    return(cube);
  }
}
/*********************************************************************************************************************************************************************/
float round(pi,rnd)                                                                  /*  This is the procedure that rounds a vector to the closet one whose elements */
float pi[8], rnd; {                                                                  /*  are multiples of rnd.  It is not used except in init_tree and it's          */
  int i,biggest, rndmax();                                                           /*  explanation is quite similar to that of the one above.                      */
  register int j=1;
  static int round[8];
  float sum=0.0, tmp,cube=0.0;
  static float delta[8];
  
  for(i=0;i<8;i++) {
    tmp=(float)(rint(pi[i]/rnd)*rnd);
    delta[i]=fabs(pi[i]-tmp);
    if (pi[i]<tmp)
      round[i]=1;
    else 
      round[i]=0;
    pi[i]=tmp;
    sum+=pi[i];

  }
  if (sum>(1.0+rnd-0.001)) {
    while (sum>1.0) {
      biggest=rndmax(delta,round,1);
      pi[biggest]-=rnd;
      delta[biggest]=0.0;
      sum-=rnd;
    }
  }
  if (sum<(1.0-rnd+0.001)) {
    while (sum<1.0) {
      biggest=rndmax(delta,round,0);
      pi[biggest]+=rnd;
      delta[biggest]=0.0;
      sum+=rnd;
    }
  }
  for(i=0;i<8;i++) {
    cube+=(j*j*j)*pi[i]*pi[i]*pi[i];
    j++;
  }
  return(cube);
}
/*****************************************************************************/
int rndmax(delta, round, val)
float delta[];
int round[],val; {

  int i, choice=0;
  float diff=0.0;

  for(i=0;i<8;i++) {
    if(round[i]==val) {
      if(delta[i]>diff) {
        diff=delta[i];
        choice=i;
      }
    }
  }
  return(choice);
}