// tree.c           PROCEDURES FOR BINARY TREE
// (GNU General Public License software, (C) Sivia and Skilling 2006)
// Example: m=5 rates (r0,r1,r2,r3,r4) are stored using
//          mplus=8 cells, and are preceded by partial sums.
//    +-------------------------------------------------------+
//  T |1                   r0+r1+r2+r3+r4                     |
//    |-------------------------------------------------------|
//  R |2       r0+r1+r2+r3        |3           r4             |
//    |---------------------------|---------------------------|
//  E |4   r0+r1    |5   r2+r3    |6    r4      |7     0      |
//    |-------------|-------------|-------------|-------------|
//  E |8  r0 |9  r1 |10 r2 |11 r3 |12 r4 |13  0 |14  0 |15  0 |
//    +------+------+------+------+------+------+------+------+
//          Tree[0] holds the value of mplus, here 8.
/*__________________________________________________________________*/
void PutRate(    // Update binary tree with new rate
double* Tree,    // Binary tree of rates
int     j,       // Input cell (not including mplus)
double  r)       // New value of rate
{
    j += (int)Tree[0];
    Tree[j] = r;
    for( ; j > 1; j >>= 1 )
        Tree[j>>1] = Tree[j] + Tree[j^1];
}
/*__________________________________________________________________*/
int GetRate(     // Select cell (not including mplus) with random rate
double* Tree)    // Binary tree of rates
{    // This procedure is crafted never to select a cell having rate=0
    double f;
    int    j;
    int    mplus = (int)Tree[0]; // Assuming Tree[1] > 0, then
    f = Tree[1] * UNIFORM;       // 0 < f <= Tree[1] for any rounding.
    for( j = 1; j < mplus; )
    {              // Enter with
        j <<= 1;   // f <= Tree[parent] = Tree[j]+Tree[j+1] >= Tree[j]
        if( f > Tree[j] ) // Only possible if Tree[j] < Tree[parent],
        {                        // so that Tree[j+1] > 0,
            f -= Tree[j++];      // hence j+1 is a safe destination.
            if( f > Tree[j] )    // Rounding error occasionally
                f = Tree[j];     // matters, so keep f <= Tree[j].
        }                        // Either way, f > 0 by construction.
    }                     // Else 0 < f <= Tree[j], so Tree[j] > 0,
    return j - mplus;     // hence the entry j was a safe destination.
}