Теряется точность при вычислении с double

#include <iostream>
#include <iomanip>
#include "mpreal.h"
 
int main() {
    mpfr::mpreal::set_default_prec( mpfr::digits2bits( 100 ) );
 
    mpfr::mpreal pi = 3.1415926535897932384626433832795;
    mpfr::mpreal a3 = 1;
    mpfr::mpreal b3 = -8.800775976179623e+014;
    mpfr::mpreal c3 = -1.794402976530432e+016;
    mpfr::mpreal d3 = -9.331775087260467e+016;
 
    mpfr::mpreal p3 = (3 * a3 * c3 - b3 * b3) / (3 * a3 * a3);
    mpfr::mpreal q3 = (2 * b3 * b3 * b3 - 9 * a3 * b3 * c3 + 27 * a3 * a3 * d3) / (27 * a3 * a3 * a3);
 
    mpfr::mpreal q32 = q3 / 2;
    mpfr::mpreal p33 = p3 / 3;
    mpfr::mpreal delta = q32 * q32 + p33 * p33 * p33;
 
    mpfr::mpreal fi, y0;
    if ( delta > 0 ) {
        y0 = mpfr::root( -q3 / 2 + mpfr::sqrt( delta ), 3 ) +
             mpfr::root( -q3 / 2 - mpfr::sqrt( delta ), 3 );
    } else {
        if ( delta < 0 ) {
            if ( q3 < 0 ) {
                fi = mpfr::atan2( -mpfr::sqrt( delta ), -q3 / 2 );
            } else {
                if ( q3 > 0 ) {
                    fi = mpfr::atan2( -mpfr::sqrt( delta ), -q3 / 2 ) + pi;
                } else {
                    fi = pi / 2;
                }
            }
            y0 = 2 * mpfr::sqrt( -p3 / 2 ) * mpfr::cos( fi / 3 );
        } else {
            y0 = mpfr::root( -q3 / 2, 3 ) * 2;
        }
    }
 
    std::cout << std::setprecision( 32 );
    std::cout << y0 - b3 / 3 << std::endl;
}