I thought the following code clean up might be worth illustrating. In the dcdflib method beta_inc, there are a number of constructs similar to the following:

    if (x == 0.0) goto S210;
    if (y == 0.0) goto S230;
    if (a == 0.0) goto S240;
    if (b == 0.0) goto S220;

S210:
    //
     //  TERMINATION OF THE PROCEDURE
     //
     if (a == 0.0)
    {
        // x == 0 && a == 0
         ierr = 6;
        return;
    }
S220:
    w = 0.0;
    w1 = 1.0;
    return;
S230:
    if (b == 0.0)
    {
        ierr = 7;
        return;
    }
S240:
    w = 1.0;
    w1 = 0.0;
    return;

When I approach modifying this sort of code, I find it helpful to trace code execution and add comments to the code to help me keep track of the conditions that apply when a code block is executed. Here is the code with my comments added:

    // If we get to this point, we know that x + y == 1
     // a and b are not negative
     // a and b are not both zero
     if (x == 0.0) goto S210;
    if (y == 0.0) goto S230;
    if (a == 0.0) goto S240;
    if (b == 0.0) goto S220;

S210:
    //
     //  TERMINATION OF THE PROCEDURE
     //
     // At this point, all we know is x == 0.0
     // The next test will fall through to label S220 if a != 0
     if (a == 0.0)
    {
        // x == 0 && a == 0
         ierr = 6;
        return;
    }
S220:
    // Two conditions can cause this code to execute
     // x == 0 and a != 0
     // or b == 0
     // There is a timing issue. The test for b == 0 occurs after
     // 'if' statements 2 & 3 have executed.
     w = 0.0;
    w1 = 1.0;
    return;
S230:
    // At this point, all we know is y == 0.0
     // The next test will fall through to label S240 if b != 0
     if (b == 0.0)
    {
        ierr = 7;
        return;
    }
S240:
    // Two conditions can cause this code to execute
     // y == 0 and b != 0
     // or a == 0
     w = 1.0;
    w1 = 0.0;
    return;

I used my notes from the above code block to layout a new set of if test blocks. The final code block is shown below. This code is not intended to be efficient, but to exactly imitate the existing code using modern constructs:

    if (x == 0.0)
    {
        if (a == 0.0)
        {
            ierr = 6;
            return;
        }
        else
        {
            w = 0.0;
            w1 = 1.0;
            return;
        }
    }

    if (y == 0.0)
    {
        if (b == 0.0)
        {
            ierr = 7;
            return;
        }
        else
        {
            w = 1.0;
            w1 = 0.0;
            return;
        }
    }

    if (a == 0.0)
    {
        w = 1.0;
        w1 = 0.0;
        return;
    }

    if (b == 0.0)
    {
        w = 0.0;
        w1 = 1.0;
        return;
    }

This code reduces to the following:

    if (x == 0.0)
    {
        if (a == 0.0)
        {
            ierr = 6;
            return;
        }
        else
        {
            w = 0.0;
            w1 = 1.0;
            return;
        }
    }

    if ((y == 0.0) || (a == 0.0))
    {
        if (b == 0.0)
        {
            ierr = 7;
            return;
        }
        else
        {
            w = 1.0;
            w1 = 0.0;
            return;
        }
    }

    if (b == 0.0)
    {
        w = 0.0;
        w1 = 1.0;
        return;
    }
 

How can we be assured that this code faithfully reproduces the original code paths? We have to return to my notes concerning what we know about the parameters a, b, x, and y. The notes are simply paraphrases of input parameter tests performed before any calculations take place. We know that x + y = 1 and we know that if a = 0, then b cannot be 0, similarly, if b = 0, then a cannot be zero.

    // If we get to this point, we know that x + y == 1
     // a and b are not negative
     // a and b are not both zero

The code does bring up an interesting issue. Should we continue to return an ierr for backward compatibility, or should the code raise execptions when input errors occur? As always, commments are welcome.

As this topic unfolds, I'll make my notes available here.