What is the correct algorithm to perform double-float division?

Reviewing the Pseudocode Algorithm as written in the book appears to support the C++ implementation of this algorithm, although my unfamiliarity with Cg means I can't prove that that implementation is correct for Cg.

So breaking down these steps into plain english:

1. The function takes two parameters, each of which are [pseudo-]double precision floating point values, and where the second parameter is not equal to 0
2. The variable x_n is assigned the arithmetic reciprocal of the higher order component of the [pseudo-]double divisor, calculated using single precision floating point math
3. The variable y_n is assigned the product of the higher order component of the [pseudo-]double dividend and x_n, calculated using single precision floating point math
4. The product of the [pseudo-]double Divisor and y_n is calculated
   • This is the first tricky part, because the paper doesn't describe an algorithm for [pseudo-]double x single multiplication. We can see in the Cg algorithm that the Cg algorithm clearly maps to this step 1-to-1, but the Cg rules for promoting a scalar value to a vector value are unknown.
   • What we can say, however, is that we do have a function for multiplying a double by a double, and a single can be promoted to a double by padding its lower order component with 0, so we can do that.
5. The difference between the Dividend and the product calculated in step 4 is calculated, and only the higher order component is kept as a single-precision floating point value
   • What makes this tricky is that the paper doesn't describe an algorithm for subtraction. However, it does describe an algorithm for converting a [IEEE754-]double into a [pseudo-]double, and an observation we can make is that negative [IEEE754-]doubles, when converted, have negative values for both its higher order and lower order components. So logically, a [pseudo-]double can be negated by simply negating both of its components. And a negated number added is mathematically equivalent to subtraction, so we can build a subtraction algorithm using this knowledge.
6. The product of x_n and step 5 is performed, preserving the extended precision that would otherwise be lost in a single x single multiplication.
   • The twoProd function exists for exactly this purpose.
7. The sum of step 6 and y_n is calculated
   • Again, we can use the [pseudo-]double addition algorithm if we simply promote y_n to a [pseudo-]double by padding the lower order component with 0
8. The result of step 7 is the returned value

So understanding this algorithm we can map each of these steps directly to the C++ algorithm I wrote:

//(1) Takes two [pseudo-]doubles, returns a [pseudo-]double
float2 df64_div(float2 B, float2 A) {
    //(2) single float divided by single float
    float xn = 1.0f / A.x;
    // (3) single float multiplied by single float
    float yn = B.x * xn;
    //                        (4) double x double multiplication
    //                                       (4a) yn promoted to [pseudo-]double
    //            (5) subtraction                           (5a) only higher order component kept
    float diff = (df64_diff(B, df64_mult(A, float2(yn, 0)))).x;
    // (6) single x single multiplication with extra precision preserved using twoProd
    float2 prod = twoProd(xn, diff);
    // (7) adding higher-order division to lower order division
    //              (7a) yn promoted to [pseudo-]double
    // (8) value is returned
    return df64_add(float2(yn, 0), prod);
}

float2 df64_diff(float2 a, float2 b) {
    //                 (5a) negating both components is a logical negation of the whole number
    return df64_add(a, float2(-b.x, -b.y));
}

From this, we can conclude that this is a correct implementation of the algorithm described in this paper, upheld by some testing I've done to validate that performing these operations in this manner yields results that appear to be correct.

Xirema's answer provides a faithful rendering of Thall's high-radix long-hand division algorithm into C++. Based on fairly extensive testing against a higher-precision reference, I found its maximum relative error to be on the order of 2^-45, provided there are no underflows in intermediate computation.

On platforms that provide a fused multiply-add operation (FMA), the following Newton-Raphson-based division algorithm due to Nagai et. al. is likely to be more efficient and achieves identical accuracy in my testing, that is, maximum relative error of 2^-45.

/*
  T. Nagai, H. Yoshida, H. Kuroda, Y. Kanada, "Fast Quadruple Precision 
  Arithmetic Library on Parallel Computer SR11000/J2." In: Proceedings 
  of the 8th International Conference on Computational Science, ICCS '08, 
  Part I, pp. 446-455.
*/
float2 div_df64 (float2 a, float2 b)
{
    float2 t, c;
    float r, s;
    r = 1.0f / b.x;
    t.x = a.x * r;
    s = fmaf (-b.x, t.x, a.x);
    t.x = fmaf (r, s, t.x);
    t.y = fmaf (-b.x, t.x, a.x);
    t.y = a.y + t.y;
    t.y = fmaf (-b.y, t.x, t.y);
    s = r * t.y;
    t.y = fmaf (-b.x, s, t.y);
    t.y = fmaf (r, t.y, s);
    c.x = t.x + t.y;
    c.y = (t.x - c.x) + t.y;
    return c;
}