Numerical Instability in Borwein's Quartic Algorithm for Pi Calculation

[DeltaResero]

Description:

I've been trying to implement the Borwein quartic algorithm for calculating Pi as another possible method of calculating Pi for WPCPP, but have encountered significant issues with precision loss and numerical instability. Despite multiple refinements, the algorithm is consistently returning incorrect results (around 3.092...), even after increasing precision and iterations dramatically. y was even modified to using a more stable form to attempt to avoid the catastrophic cancellation (sqrt(1 - y^4) issue outline here on Stack Exchange: https://math.stackexchange.com/questions/4325439/cancellation-problem-in-borweins-quartic-algorithm-for-computing-pi

 

Code Attempt:

mpf_class calculate_pi_borwein() {
    mpf_set_default_prec(1024);  // Manually set here for high precision for more accuracy while testing
    mpf_class a = mpf_class(6) - 4 * sqrt(mpf_class(2));
    mpf_class y = sqrt(mpf_class(2)) - 1;
    int iterations = 5; // Tried as high as 100 and seen no difference

    for (int i = 0; i < iterations; ++i) {
        mpf_class y2 = y * y;
        mpf_class y4 = y2 * y2;
        mpf_class sqrt1_y4 = sqrt(1 - y4);
        mpf_class y_new = (1 - sqrt1_y4) / (1 + sqrt1_y4);
        mpf_class one_plus_y_new = 1 + y_new;
        mpf_class one_plus_y_new_pow_4;
        mpf_pow_ui(one_plus_y_new_pow_4.get_mpf_t(), one_plus_y_new.get_mpf_t(), 4);
        mpf_class two_pow;
        mpf_pow_ui(two_pow.get_mpf_t(), mpf_class(2).get_mpf_t(), 2 * (i + 1));
        mpf_class a_new = a * one_plus_y_new_pow_4 - two_pow * y_new * (1 + y_new + y_new * y_new);
        a = a_new;
        y = y_new;
    }
    return 1 / a;
}

 

Potential Root Cause:

• Despite attempts to stabilize the recursive terms and increase precision, the algorithm still suffers from rapid precision loss, particularly in the recursive update for y and the a term

• Numerical instability in the square root term (sqrt(1 - y^4)) and other recursive terms may be causing the collapse of y prematurely

Any insights on why this might be happening and what adjustments can be made to fix the numerical stability issues would be much appreciated. Until this is fixed properly, the Borwein quartic algorithm will be left out of WPCPP.