#!/usr/bin/php
<?php

define('MAX_POWER', 10);
define('MAX_DISTANCE', 99);
ini_set('memory_limit', '2G');

# ---

function bcpow2($x) {
	return bcmul2($x,$x);
}

function bcpow3($x) {
	return bcmul3($x,$x,$x);
}

function bcmul3($a,$b,$c) {
	return bcmul2($a,bcmul2($b,$c));
}

function bcsub4($a,$b,$c,$d) {
	return bcsub2(bcsub2(bcsub2($a,$b),$c),$d);
}

function bcabs($x) {
	return str_replace('-', '', $x);
}

function bcadd2($a,$b) {
	return gmp_add($a,$b);
}

function bcsub2($a,$b) {
	return gmp_sub($a,$b);
}

function bcmul2($a,$b) {
	return gmp_mul($a,$b);
}

function quad_mul($a, $b) {
	// Powers explained here: https://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/functions/power/index.htm

	# $c0 = $a[0]*$b[0] - $a[1]*$b[1] - $a[2]*$b[2] - $a[3]*$b[3];
	# $c1 = $a[0]*$b[1] + $a[1]*$b[0] + $a[2]*$b[3] - $a[3]*$b[2];
	# $c2 = $a[0]*$b[2] - $a[1]*$b[3] + $a[2]*$b[0] + $a[3]*$b[1];
	# $c3 = $a[0]*$b[3] + $a[1]*$b[2] - $a[2]*$b[1] + $a[3]*$b[0];

	$c0 = bcsub2(bcsub2(bcsub2(bcmul2($a[0],$b[0]),bcmul2($a[1],$b[1])),bcmul2($a[2],$b[2])),bcmul2($a[3],$b[3]));
	$c1 = bcsub2(bcadd2(bcadd2(bcmul2($a[0],$b[1]),bcmul2($a[1],$b[0])),bcmul2($a[2],$b[3])),bcmul2($a[3],$b[2]));
	$c2 = bcadd2(bcadd2(bcsub2(bcmul2($a[0],$b[2]),bcmul2($a[1],$b[3])),bcmul2($a[2],$b[0])),bcmul2($a[3],$b[1]));
	$c3 = bcadd2(bcsub2(bcadd2(bcmul2($a[0],$b[3]),bcmul2($a[1],$b[2])),bcmul2($a[2],$b[1])),bcmul2($a[3],$b[0]));

	return [ $c0, $c1, $c2, $c3 ];
}

function quad_int_power($quad, $n) {
	$a = $quad[0];
	$b = $quad[1];
	$c = $quad[2];
	$d = $quad[3];
	if ($n == 1) {
		$a_ = $a;
		$b_ = $b;
		$c_ = $c;
		$d_ = $d;
		return array($a_,$b_,$c_,$d_);
	} else if ($n == 2) {
		$a_ = bcsub4(bcpow2($a),bcpow2($b),bcpow2($c),bcpow2($d));
		$b_ = bcmul3('2',$a,$b);
		$c_ = bcmul3('2',$a,$c);
		$d_ = bcmul3('2',$a,$d);
		return array($a_,$b_,$c_,$d_);
	} else if ($n == 3) {
		$a_ = bcsub4(bcpow3($a), bcmul3('3',bcpow2($b),$a),
		                         bcmul3('3',bcpow2($c),$a),
		                         bcmul3('3',bcpow2($d),$a));
		$b_ = bcsub4(bcmul3('3',bcpow2($a),$b), bcpow3($b),
		             bcmul3('3',bcpow2($c),$a),
		             bcmul3('3',bcpow2($d),$a));
		$c_ = bcsub4(bcmul3('3',bcpow2($a),$c), bcpow3($c),
		             bcmul3('3',bcpow2($b),$c),
		             bcmul3('3',bcpow2($d),$c));
		$d_ = bcsub4(bcmul3('3',bcpow2($a),$d), bcpow3($d),
		             bcmul3('3',bcpow2($b),$d),
		             bcmul3('3',bcpow2($c),$d));
		return array($a_,$b_,$c_,$d_);
	} else {
		assert($n >= 1);
		$res = $quad;
		for ($p=2; $p<=$n; $p++) $res = quad_mul($res, $quad);
		return $res;
	}
}

function quad_is_immortal($quad, $power) {
	$quadpow = quad_int_power($quad, $power);
	for ($i=0; $i<4; $i++) {
		$x1 = bcabs($quad[$i]);
		$x2 = bcabs($quadpow[$i]);
		if (substr($x2, -strlen($x1)) != $x1) return false;
	}
	return $quadpow;
}

function _go($a, $b, $c, $d, $power, $distance, &$progr, &$found) {
	$quad = [$a, $b, $c, $d];
	$quadpower = quad_is_immortal($quad, $power);
	if ($quadpower !== false) {
		$a_ = $quadpower[0];
		$b_ = $quadpower[1];
		$c_ = $quadpower[2];
		$d_ = $quadpower[3];
		if (($c!=0) || ($d!=0)) {
			$str = "QUAD ";
		} else if ($b!=0) {
			$str = "CPLX ";
		} else {
			$str = "REAL ";
		}
		$str .= "| Distance=$distance | ";
		$str .= "[$a + $b i + $c j + $d k] ^ $power = [$a_ + $b_ i + $c_ j + $d_ k]";
		if (!in_array($str,$found)) {
#			echo "$str\n";
			$found[] = $str;
		}
	}
	#if ($progr++ % 10000 == 0) {
	#	echo "(Dist=$distance) Processing ($a, $b, $c, $d) with power $power                 \r";
	#}
}

if (file_exists(__DIR__.'/savestate.txt')) {
	$start_distance = trim(file_get_contents(__DIR__.'/savestate.txt'));
} else {
	$start_distance = 0;
}
for ($distance=$start_distance; $distance<=MAX_DISTANCE; $distance++) {
	$str = "--------- Distance: $distance ---------";
	echo "                                                                              \r";
	echo "$str\n";
	$found = [ $str ];
	$progr = 0;

	$a = $distance;
	for ($b=0; $b<=$distance; $b++) {
	for ($c=0; $c<=$distance; $c++) {
	for ($d=0; $d<=$distance; $d++) {
	if ($d == 0) echo "\r".round(memory_get_usage()/1024/1024, 2)." MiB (".count($found)." found); ($a), $b, $c, $d                                ";
	for ($power=2; $power<=MAX_POWER; $power++) {
	_go($a, $b, $c, $d, $power, $distance, $progr, $found);
	_go($b, $a, $c, $d, $power, $distance, $progr, $found);
	_go($b, $c, $a, $d, $power, $distance, $progr, $found);
	_go($b, $c, $d, $a, $power, $distance, $progr, $found);
	}
	}
	}
	}

	file_put_contents(__DIR__.'/result.txt', implode("\n",$found)."\n", FILE_APPEND);
	file_put_contents(__DIR__.'/savestate.txt', $distance+1);
}