[PHP-CLI] Riemann zeta function zeroes "test" in ℍ

<?php
//not really efficient but it does the job relatively well enough lol
function conj($q) { return [$q[0],-$q[1],-$q[2],-$q[3]]; }
function mul($q1,$q2) { return [$q1[0]*$q2[0]-$q1[1]*$q2[1]-$q1[2]*$q2[2]-$q1[3]*$q2[3],
                                $q1[1]*$q2[0]+$q1[0]*$q2[1]-$q1[3]*$q2[2]+$q1[2]*$q2[3],
                                $q1[2]*$q2[0]+$q1[3]*$q2[1]+$q1[0]*$q2[2]-$q1[1]*$q2[3],
                                $q1[3]*$q2[0]-$q1[2]*$q2[1]+$q1[1]*$q2[2]+$q1[0]*$q2[3]]; }
function add($q1,$q2) { return [$q1[0]+$q2[0],$q1[1]+$q2[1],$q1[2]+$q2[2],$q1[3]+$q2[3]]; }
function opp($q) { return [-$q[0],-$q[1],-$q[2],-$q[3]]; }
function norm($q) { return sqrt(pow($q[0],2)+pow($q[1],2)+pow($q[2],2)+pow($q[3],2)); }
function a($q) { return mul([1/(norm($q)**2),0,0,0],conj($q)); }
function b($q) { return [0,$q[1],$q[2],$q[3]]; }
function c($q) { if (b($q) != [0,0,0,0]) {
                        return mul(
                            [exp($q[0]),0,0,0],
                            add (
                                [cos(norm(b($q))),0,0,0],
                                opp (
                                        mul([sin(norm(b($q)))/norm(b($q)),0,0,0],b($q))
                                    )
                                )
                        ); } else {
                            return [exp($q[0]),0,0,0];
                        }
                }
function zeta($q,$up) {
    $eta = 0;
    for ($i = 1; $i <= $up; $i++) {
        if ($i % 2 == 1) { $eta = add($eta,(c(mul($q,[-log($i),0,0,0])))); }
        if ($i % 2 == 0) { $eta = add($eta,opp((c(mul($q,[-log($i),0,0,0]))))); }
    }
    return  mul(
                a(add(
                      [1,0,0,0],
                      opp(c(add(
                                [log(2),0,0,0],
                                opp(mul(
                                        $q,
                                        [log(2),0,0,0]
                                ))
                      )))
                 )),
                $eta
            );
}
function show($q) { return $q[0]."+".$q[1]."i+".$q[2]."j+".$q[3]."k"; }
print "This program computes random instances of Z(a+bi+cj+dk) and gives potential zeroes approximations\n";
print "Please use a precision sufficiently larger than 0.0000000001 for best accuracy\n";
$scal = readline("Precision (smallest decimal k s.t. a/k, b/k, c/k, d/k are integers) : ");
$up = readline("Zeta sum upper bound (the larger the better) : ");
$maxa = readline("0  < a <= ");
$maxb = readline("0 <= b <= ");
$maxc = readline("0 <= c <= ");
$maxd = readline("0 <= d <= ");
$bruh = [];
while (0 < 1) {
    $p = [rand(1,$maxa/$scal)*$scal,rand(0,$maxb/$scal)*$scal,rand(0,$maxc/$scal)*$scal,rand(0,$maxd/$scal)*$scal];
    $q = norm(zeta($p,$up));
    if ($q < $scal*ceil(0.0000000001/$scal) & !in_array($p,$bruh)) { print "Z(".show($p).") = $q                \n"; $bruh[count($bruh)] = $p; }
    else { print "Z(".show($p).") = $q                \r"; }
}
?>