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

x1 = 3291825.45;
y1 = -7473079.35;
z1 = 1;
 
x2 = 0;
y2 = -8000000+1;
z2 = 1;
 
x3 = -2127126.49;
y3 = -7519740.99+1;
z3 = 1;
 
R11 = 0;
R21 = -1847288.7;
R31 = -1055041.77;
L   = 6378137;
 
 
 
% syms x y R z ;
 
 
S = solve('x^2+y^2+z^2 = 6378137^2',...
          '(x-(3291825.45))^2+(y-(-7473079.35+1))^2+(z-(0+1))^2=(R)^2',...
          '(x-(0))^2+(y-(-8000000+1))^2+(z-(0+1))^2=(R+(-1847288.7))^2',...
          '(x-(-2127126.49))^2+(y-(-7519740.99+1))^2+(z-(0+1))^2=(R+(-1055041.77))^2');
 
 S.R
 
 
 for i=1:4
 S.x(i) =   S.x(i) - 1;
 S.y(i) =   S.y(i) - 1;
 S.z(i) =   S.z(i) - 1;
 end
 
 
A1 = 2*(x2-x1);
B1 = 2*(y2-y1);
C1 = 2*(z2-z1);
f1 = x2^2+y2^2+z2^2 - (x1^2+y1^2+z1^2) - R21^2;
g1 = -2*R21;
 
A2 = 2*(x3-x1);
B2 = 2*(y3-y1);
C2 = 2*(z3-z1);
f2 = x3^2+y3^2+z3^2 - (x1^2+y1^2+z1^2) - R31^2;
g2 = -2*R31;
 
A3 = 2*x1;
B3 = 2*y1;
C3 = 2*z1;
f3 = (x1^2+y1^2+z1^2) + L^2;
g3 = -1;
 
D1 = (A1*(B2*C3 - B3*C2) + A2*(B3*C1 - B1*C3) + A3*(B1*C2 - B2*C1));
D2 = D1^2;
 
m4_1 = (A1*B2*g3-A2*B1*g3)^2;
m4_2 = (A1*C2*g3-A2*C1*g3)^2;
m4_3 = (B1*C2*g3-B2*C1*g3)^2;
 
m3_1 = 2*(A1*B2*g3-A2*B1*g3)*(A1*B3*g2-A2*B3*g1-A3*B1*g2+A3*B2*g1);
m3_2 = 2*(A1*C2*g3-A2*C1*g3)*(A1*C3*g2-A2*C3*g1-A3*C1*g2+A3*C2*g1);
m3_3 = 2*(B1*C2*g3-B2*C1*g3)*(B1*C3*g2-B2*C3*g1-B3*C1*g2+B3*C2*g1);
 
m2_1 = (B1*C3*g2-B2*C3*g1-B3*C1*g2+B3*C2*g1)^2;
m2_2 = 2*(B1*C2*g3-B2*C1*g3)*(B1*C2*f3-B1*C3*f2-B2*C1*f3+B2*C3*f1+B3*C1*f2-...
        B3*C2*f1);
m2_3 = (A1*C3*g2-A2*C3*g1-A3*C1*g2+A3*C2*g1)^2;
m2_4 = 2*(A1*C2*g3-A2*C1*g3)*(A1*C2*f3-A1*C3*f2-A2*C1*f3+A2*C3*f1+A3*C1*f2-...
        A3*C2*f1);
m2_5 = (A1*B3*g2-A2*B3*g1-A3*B1*g2+A3*B2*g1)^2;
m2_6 = 2*(A1*B2*g3-A2*B1*g3)*(A1*B2*f3-A1*B3*f2-A2*B1*f3+A2*B3*f1+...
        A3*B1*f2-A3*B2*f1);
    
m1_1 = 2*(A1*C3*g2-A2*C3*g1-A3*C1*g2+A3*C2*g1)*(A1*C2*f3-A1*C3*f2-A2*C1*f3+...
        A2*C3*f1+A3*C1*f2-A3*C2*f1);
m1_2 = 2*(B1*C3*g2-B2*C3*g1-B3*C1*g2+B3*C2*g1)*(B1*C2*f3-B1*C3*f2-B2*C1*f3+...
        B2*C3*f1+B3*C1*f2-B3*C2*f1);
m1_3 = 2*(A1*B3*g2-A2*B3*g1-A3*B1*g2+A3*B2*g1)*(A1*B2*f3-A1*B3*f2-A2*B1*f3+...
        A2*B3*f1+A3*B1*f2-A3*B2*f1);
 
m0_1 = (A1*B2*f3-A1*B3*f2-A2*B1*f3+A2*B3*f1+A3*B1*f2-A3*B2*f1)^2;
m0_2 = (A1*C2*f3-A1*C3*f2-A2*C1*f3+A2*C3*f1+A3*C1*f2-A3*C2*f1)^2;
m0_3 = (B1*C2*f3-B1*C3*f2-B2*C1*f3+B2*C3*f1+B3*C1*f2-B3*C2*f1)^2;
m0_4 = L^2;
 
a = (m4_1 + m4_2 + m4_3)/D2;
b = (-m3_1 - m3_2 - m3_3)/D2;
c = (m2_6 + m2_5 + m2_4 + m2_3 + m2_2 + m2_1)/D2;
d = (-m1_1 - m1_2 - m1_3)/D2;
e = (m0_1 + m0_2 + m0_3)/D2 - m0_4;
 
p = (8*a*c - 3*b^2)/(8*a^2);
q = (8*a^2*d+b^3-4*a*b*c)/(8*a^3);
r = (16*a*b^2*c-64*a^2*b*d-3*b^4+256*a^3*e)/(256*a^4);
 
 
if q ~= 0
    
       
    a33 = 1
    b33 = p
    c33 = (p^2-4*r)/4
    d33 = - q^2/8
    
    T = solve('t^3 -880256255399764*t^2 -1.79792141373932e+016*t -9.24139523145728e+016 = 0')
    ff = double(T(1))
    p33 = (3*a33*c33-b33^2)/(3*a33^2);
    q33 = (2*b33^3-9*a33*b33*c33+27*a33^2*d33)/(27*a33^3)
        
    delta = (q33/2)^2 + (p33/3)^3
     
    if delta > 0                   
       y0 = (-q33/2 + sqrt(delta))^(1/3) + (-q33/2 - sqrt(delta))^(1/3);
    
    else        
        
        if delta < 0
            
            if q33 < 0
                fi = atan2(-sqrt(delta),-q33/2);
            else
                if q33 > 0
                    fi = atan2(-sqrt(delta),-q33/2) + pi;
                else
                    fi = pi / 2;
                end
            y0 = 2*sqrt(-p33/2)*cos(fi/3);
            end
             
        else
            y0 = 2*(-q33/2)^(1/3);
        end
 
    end
   
    z0 =  y0 - b33 / (3.0*a33)
    z0 = 880256255399784.4;
    
    
    y_1 = (sqrt(2*z0) - sqrt(2*z0-4*(p/2+z0 + q / (2*sqrt(2*z0)) ) ) )/2;    
    
    y_2 = (sqrt(2*z0) + sqrt(2*z0-4*(p/2+z0 + q / (2*sqrt(2*z0)) ) ) )/2;    
    
    y_3 = (-sqrt(2*z0) - sqrt(2*z0-4*(p/2+z0 - q / (2*sqrt(2*z0)) ) ) )/2;    
    
    y_4 = (-sqrt(2*z0) + sqrt(2*z0-4*(p/2+z0 - q / (2*sqrt(2*z0)) ) ) )/2;
            
else
    
    y_1 =  sqrt((-p-sqrt(p^2-4*r))/2);
    y_2 = -sqrt((-p-sqrt(p^2-4*r))/2);
    y_3 =  sqrt((-p+sqrt(p^2-4*r))/2);
    y_4 = -sqrt((-p+sqrt(p^2-4*r))/2);
    
end
 
R1 = y_1 - b/(4*a);
R2 = y_2 - b/(4*a);
R3 = y_3 - b/(4*a);
R4 = y_4 - b/(4*a);
R = [R1 R2 R3 R4];
 
 
for i = 1:4
    
    if R(i)  > 0
        
        R(i) 
    
        x = ((B1*C2*g3-B2*C1*g3)*R(i)^2 - (B1*C3*g2-B2*C3*g1-B3*C1*g2+B3*C2*g1)*R(i) +...
             (B1*C2*f3-B1*C3*f2-B2*C1*f3+B2*C3*f1+B3*C1*f2-B3*C2*f1))/D1
 
        y = (-(A1*C2*g3-A2*C1*g3)*R(i)^2 + (A1*C3*g2-A2*C3*g1-A3*C1*g2+A3*C2*g1)*R(i) -...
            (A1*C2*f3-A1*C3*f2-A2*C1*f3+A2*C3*f1+A3*C1*f2-A3*C2*f1))/D1
 
        z = ((A1*B2*g3-A2*B1*g3)*R(i)^2+(-A1*B3*g2+A2*B3*g1+A3*B1*g2-A3*B2*g1)*R(i)+...
            (A1*B2*f3-A1*B3*f2-A2*B1*f3+A2*B3*f1+A3*B1*f2-A3*B2*f1))/D1
    end
    
end

T = solve('t^3 -880256255399764*t^2 -1.79792141373932e+016*t -9.24139523145728e+016 = 0')