四次方程求解器不适用于大系数

Question

方程模式：

ax^4+bx^3+cx^2+dx+e=0

求解器代码：

<?php

   function ferrariQuartic($a, $b, $c, $d, $e) {
    // Normalize the coefficients
    if ($a == 0) {
        throw new Exception("Coefficient 'a' cannot be zero for a quartic equation.");
    }

    // Step 1: Normalize the equation
    $b /= $a;
    $c /= $a;
    $d /= $a;
    $e /= $a;

    // Step 2: Substitute x = y - (b / 4)
    $p = $c - (3 * $b**2) / 8;
    $q = ($b * $d) / 2 - ($b**3) / 8 + $e;
    
    // Step 3: Solve y^4 + p*y^2 + q = 0 as a quadratic in y^2
    $D = $p**2 / 4 - $q; // Discriminant
    if ($D < 0) {
       // throw new Exception("The quartic equation has complex roots.");
    }

    // Roots for y^2
    $y2_1 = -$p / 2 + sqrt($D);
    $y2_2 = -$p / 2 - sqrt($D);
    
    // Step 4: Solve for y
    $y1 = sqrt($y2_1);
    $y2 = sqrt($y2_2);
    
    // Step 5: Back substitute to find x
    $roots = [];
    $roots[] = $y1 - ($b / 4);
    $roots[] = -$y1 - ($b / 4);
    $roots[] = $y2 - ($b / 4);
    $roots[] = -$y2 - ($b / 4);
    
    // Return the roots
    return array_filter($roots, function($root) {
        return is_finite($root); // Filter out any infinite or NaN values
    });
}



// Working with small coefficients
$a = 1; // Coefficient of x^4
$b = 0; // Coefficient of x^3
$c = -5; // Coefficient of x^2
$d = 0; // Coefficient of x
$e = 4; // Constant term


$roots = ferrariQuartic($a, $b, $c, $d, $e);

$lefHandSide = $a*($roots[1]**4) +  $b*($roots[1]**3)+ $c*($roots[1]**2)+$d*($roots[1]**1) + $e;
echo $lefHandSide;

?>

这里的系数是

// Working with small Coefficients
$a = 1; // Coefficient of x^4
$b = 0; // Coefficient of x^3
$c = -5; // Coefficient of x^2
$d = 0; // Coefficient of x
$e = 4; // Constant term

代码运行正常。 $lefHandSide = 0 = RightHandSide

但是当我输入大系数时，$lefHandSide != 0.

示例：

<?php

function ferrariQuartic($a, $b, $c, $d, $e) {
    // Normalize the coefficients
    if ($a == 0) {
        throw new Exception("Coefficient 'a' cannot be zero for a quartic equation.");
    }

    // Step 1: Normalize the equation
    $b /= $a;
    $c /= $a;
    $d /= $a;
    $e /= $a;

    // Step 2: Substitute x = y - (b / 4)
    $p = $c - (3 * $b**2) / 8;
    $q = ($b * $d) / 2 - ($b**3) / 8 + $e;
    
    // Step 3: Solve y^4 + p*y^2 + q = 0 as a quadratic in y^2
    $D = $p**2 / 4 - $q; // Discriminant
    if ($D < 0) {
       // throw new Exception("The quartic equation has complex roots.");
    }

    // Roots for y^2
    $y2_1 = -$p / 2 + sqrt($D);
    $y2_2 = -$p / 2 - sqrt($D);
    
    // Step 4: Solve for y
    $y1 = sqrt($y2_1);
    $y2 = sqrt($y2_2);
    
    // Step 5: Back substitute to find x
    $roots = [];
    $roots[] = $y1 - ($b / 4);
    $roots[] = -$y1 - ($b / 4);
    $roots[] = $y2 - ($b / 4);
    $roots[] = -$y2 - ($b / 4);
    
    // Return the roots
    return array_filter($roots, function($root) {
        return is_finite($root); // Filter out any infinite or NaN values
    });
}

$x =89565891926547004231252920425935692360644145829622209833684329913297188986597.0;
$y =12158399299693830322967808612713398636155367887041628176798871954788371653930.0; 

// Does not work with big coefficients
$a = 1; // Coefficient of x^4
$b = 0; // Coefficient of x^3
$c = -6*$x; // Coefficient of x^2
$d = -8*$y; // Coefficient of x
$e = -3*($x**2); // Constant term

$roots = ferrariQuartic($a, $b, $c, $d, $e);

$lefHandSide = $a*($roots[1]**4) +  $b*($roots[1]**3)+ $c*($roots[1]**2)+$d*($roots[1]**1) + $e;
echo $lefHandSide;


?>

这里的系数是

$x =89565891926547004231252920425935692360644145829622209833684329913297188986597.0;
$y =12158399299693830322967808612713398636155367887041628176798871954788371653930.0; 

// Does not work with big coefficients
$a = 1; // Coefficient of x^4
$b = 0; // Coefficient of x^3
$c = -6*$x; // Coefficient of x^2
$d = -8*$y; // Coefficient of x
$e = -3*($x**2); // Constant term

$lefHandSide = -7.7405416782585E+139 != RightHandSide。为什么会发生这种情况以及如何克服？

Answer 1

假设算法的实现是正确的。您可能会遇到 IEEE754 浮点精度限制。

在这种情况下，您应该使用 BCMath 扩展，它可以处理任意小数。在 PHP 8.4 中使用 BCMath 会更容易、更快，因为它有一个支持运算符重载的新类。