polyfit equivalent in PHP

polyfit.php

<?php
function polyfit($dependentValues, $independentValues, $countOfElements, $order, &$coefficients)
{
    // Declarations...
    // ----------------------------------
    $maxOrder = 5;
    
    $B = array_fill(0, $maxOrder + 1, 0.0);
    $P = array_fill(0, (2 * ($maxOrder + 1)) + 1, 0.0);
    $A = array_fill(0, (2 * ($maxOrder + 1)) * ($maxOrder + 1), 0.0);

    $x = 0.0;
    $y = 0.0;
    $powx = 0.0;

    // Verify initial conditions....
    // ----------------------------------

    // This method requires that the countOfElements > 
    // (order+1) 
    if ($countOfElements <= $order)
        return -1;

    // This method has imposed an arbitrary bound of
    // order <= maxOrder.  Increase maxOrder if necessary.
    if ($order > $maxOrder)
        return -1;

    // Begin Code...
    // ----------------------------------

    // Identify the column vector
    for ($ii = 0; $ii < $countOfElements; $ii++)
    {
        $x = $dependentValues[$ii];
        $y = $independentValues[$ii];
        $powx = 1;

        for ($jj = 0; $jj < ($order + 1); $jj++)
        {
            $B[$jj] = $B[$jj] + ($y * $powx);
            $powx = $powx * $x;
        }
    }

    // Initialize the PowX array
    $P[0] = $countOfElements;

    // Compute the sum of the Powers of X
    for ($ii = 0; $ii < $countOfElements; $ii++)
    {
        $x = $dependentValues[$ii];
        $powx = $dependentValues[$ii];

        for ($jj = 1; $jj < ((2 * ($order + 1)) + 1); $jj++)
        {
            $P[$jj] = $P[$jj] + $powx;
            $powx = $powx * $x;
        }
    }

    // Initialize the reduction matrix
    for ($ii = 0; $ii < ($order + 1); $ii++)
    {
        for ($jj = 0; $jj < ($order + 1); $jj++)
        {
            $A[($ii * (2 * ($order + 1))) + $jj] = $P[$ii + $jj];
        }

        $A[($ii * (2 * ($order + 1))) + ($ii + ($order + 1))] = 1;
    }

    // Move the Identity matrix portion of the redux matrix
    // to the left side (find the inverse of the left side
    // of the redux matrix
    for ($ii = 0; $ii < ($order + 1); $ii++)
    {
        $x = $A[($ii * (2 * ($order + 1))) + $ii];
        if ($x != 0)
        {
            for ($kk = 0; $kk < (2 * ($order + 1)); $kk++)
            {
                $A[($ii * (2 * ($order + 1))) + $kk] = 
                    $A[($ii * (2 * ($order + 1))) + $kk] / $x;
            }

            for ($jj = 0; $jj < ($order + 1); $jj++)
            {
                if (($jj - $ii) != 0)
                {
                    $y = $A[($jj * (2 * ($order + 1))) + $ii];
                    for ($kk = 0; $kk < (2 * ($order + 1)); $kk++)
                    {
                        $A[($jj * (2 * ($order + 1))) + $kk] = 
                            $A[($jj * (2 * ($order + 1))) + $kk] -
                            $y * $A[($ii * (2 * ($order + 1))) + $kk];
                    }
                }
            }
        }
        else
        {
            // Cannot work with singular matrices
            return -1;
        }
    }

    // Calculate and Identify the coefficients
    for ($ii = 0; $ii < ($order + 1); $ii++)
    {
        for ($jj = 0; $jj < ($order + 1); $jj++)
        {
            $x = 0;
            for ($kk = 0; $kk < ($order + 1); $kk++)
            {
                $x = $x + ($A[($ii * (2 * ($order + 1))) + ($kk + ($order + 1))] *
                    $B[$kk]);
            }
            $coefficients[$ii] = $x;
        }
    }

    return 0;
}

$coefficients=[];
$x=[1,2,3,4,5];
$y=[3,8,12,20,30];
$r = polyfit($x,$y, 5,2, $coefficients);
print_r($coefficients);