f03powersum.php
<?php
/************************
sum of a power e of the first x natural numbers
sum(y=0 to x: y**e) = p(x) = sum(f=e+1 downto 1: a[e][f] * x ** f)
p(x)-p(x-1) = x**e
p(x)-p(x-1) = sum(f=e+1 downto 1: a[e][f] * (x ** f - (x-1) ** f))
= - sum(f=e+1 downto 1: a[e][f] * sum(g=f-1 downto 0: binominal(f, g) x ** g * (-1) ** (f-g)))
= - sum(g=e downto 0: x ** g * sum(f=e+1 downto g+1: a[e][f] * binominal(f, g) * (-1) ** (f-g)))
yielding
for g=e: 1 = - a[e+1] * binominal(e+1, e) x ** g * (-1) ** 1 ==> = a[e][e+1] = 1 / (e+1)
for g<e: 0 = sum(f=e+1 downto g+1: a[e][f] * binominal(f, g) * (-1) ** (f-g))
= sum(f=e+1 downto g+1: a[e][f] * binominal(f, f-g) * (-1) ** (f-g))
= sum(f=e+1 downto g+1: a[e][f] * prod(h=f downto g+1: h) / (f-g)! * (-1) ** (f-g))
= sum(f=e+1 downto g+1: (cc[e, f] / (e+1) * prod(h=e+1 downto f+1: h)) * prod(h=f downto f-g+1: h) / (f-g)! * (-1) ** (f-g))
= prod(h=e+1 downto g+1) / (e+1) * sum(f=e+1 downto g+1: cc[e, f] / (f-g)! * (-1) ** (f-g))
0 = sum(f=e+1 downto g+1: cc[e, f] / (f-g)! * (-1) ** (f-g)) ==> equations depend on e-g, but not on e itself, except ccm, which thus can be simplified to c[e-f] = c [e-g-(f-g)]
0 = sum(f=e+1 downto g+1: c[e-f] / (f-g)! * (-1) ** (f-g)) ==> depends only on e-g
==> c[e-1] = sum(f=2 to e+1: c[e-f] / (f)! * (-1) ** (f)) by setting g=0 and solving for c[e-1] (f=1)
with c[-1] = 1 the assumption a[e][f] = c[e-f] / (e+1) * prod(h=e+1 downto f+1: h) is also correct for a[e][e+1]
************************/
require_once('env.php');
outBegin();
class Fract {
public static function gcd($a, $b) {
[$g, $s] = [abs($a), abs($b)];
if ($g < $s)
[$g, $s] = [$s, $g];
while (0 !== $m = $g % $s)
[$g, $s] = [$s, $m];
return $s;
}
readonly int $num, $den;
public function __construct ($n, $d=1) {
if ($d === 0) {
err("denominator 0 in ($n, $d)");
} elseif ($n === 0) {
[$this->num, $this->den] = [$n, 1];
} else {
$s = self::gcd($n, $d);
[$this->num, $this->den] = $d < 0 ? [- intdiv($n, $s), - intdiv($d, $s)] : [intdiv($n, $s), intdiv($d, $s)];
}
}
public function __tostring () {
return "$this->num" . ($this->den === 1 ? '' : "/$this->den");
}
public function reduce () {
if ($this->den === 0) {
err("denominator 0 in $this");
} elseif ($this->num === 0) {
$this->den = 1;
} else {
[$g, $s] = [abs($this->num), abs($this->den)];
if ($g < $s)
[$g, $s] = [$s, $g];
while (0 !== $m = $g % $s)
[$g, $s] = [$s, $m];
[$this->num, $this->den] = $this->den < 0 ? [- intdiv($this->num, $s), - intdiv($this->den, $s)] : [intdiv($this->num, $s), intdiv($this->den, $s)];
}
return $this;
}
public function add($s) {
try {
return is_object($s) ? (new Fract($this->num * $s->den + $this->den * $s->num, $this->den * $s->den)) : (new Fract($this->num + $this->den * $s, $this->den));
} catch (Throwable $e) {
dbq("catch add $e");
$g = self::gcd($this->den, $s->den);
$r = new Fract($this->num * ($s->den / $g) + ($d = $this->den / $g) * $s->num, $d * $s->den);
dbq("succes $this + $s = $r, gcd $g");
return $r;
}
}
public function sub($s) {
try {
return is_object($s) ? (new Fract($this->num * $s->den - $this->den * $s->num, $this->den * $s->den)) : (new Fract($this->num - $this->den * $s, $this->den));
} catch (Throwable $e) {
dbq("catch sub $e");
$g = self::gcd($this->den, $s->den);
$r = new Fract($this->num * ($s->den / $g) - ($d = $this->den / $g) * $s->num, $d * $s->den);
dbq("succes $this - $s = $r, gcd $g");
return $r;
}
}
public function mult($f) {
try {
return is_object($f) ? new Fract($this->num * $f->num, $this->den * $f->den) : new Fract($this->num * $f, $this->den);
} catch (Throwable $e) {
dbq("catch $this mult $f: $e");
is_object($f) or $f = new Fract($f);
$g = self::gcd($this->num, $f->den);
$h = self::gcd($this->den, $f->num);
$r = new Fract(($this->num/$g) * ($f->num/$h) , ($this->den/$h) * ($f->den/$g));
dbq("succes $this * $f = $r, gcd $g, $h");
return $r;
}
}
public function div($f) {
try {
return is_object($f) ? (new Fract($this->num * $f->den, $this->den * $f->num)) : (new Fract($this->num, $this->den * $f));
} catch (Throwable $e) {
dbq("catch $this div $f: $e");
is_object($f) or $f = new Fract($f);
$g = self::gcd($this->num, $f->num);
$h = self::gcd($this->den, $f->den);
$r = new Fract(($this->num/$g) * ($f->den/$h) , ($this->den/$h) * ($f->num/$g));
dbq("succes $this / $f = $r, gcd $g, $h");
return $r;
}
}
public function eq($o) {
return is_object($o) ? ($this->num === $o->num and $this->den === $o->den) : ($this->num === $o and $this->den === 1);
}
}
out("30/-45", $n = new Fract('30',-45), "$n");
out("30/45 + 3/36", $n = (new Fract(30,45))->add(new Fract(-3,-36)), "$n");
out("4/15 * 70", $n = (new Fract(4,15))->mult(70), "$n");
function powersum(int $e, int $x) {
static $c = [-1 => new Fract(1)], $a = [0 => [1=> new Fract(1)]];
if (!isset($a[$e]))
for ($t=count($a); $t <= $e; $t++) {
# ==> c[e-1] = sum(f=2 to e+1: c[e-f] / (f)! * (-1) ** (f))
$pF = 1;
$s = new Fract(0);
for ($f=2; $f <= $t+1; $f++) {
$pF *= $f;
$s = ($f % 2) ? $s->sub($c[$t-$f]->div($pF)) : $s->add($c[$t-$f]->div($pF));
}
$c[$t-1] = $s;
dbq("powersum t $t, s $s, c", array_map(fn ($v) => "$v", $c));
$pG = new Fract(1, $t+1);
for ($i=$t+1; $i>=1 ; $i--) {
$a[$t][$i] = $c[$t-$i]->mult($pG);
$pG = $pG->mult($i);
}
krsort($a[$t]);
dbq("powersum a[$t]", array_map(fn ($v) => "$v", $a[$t]));
}
$s = new Fract(0);
for ($f=$e+1;$f>0;$f--)
$s = $s->add($a[$e][$f])->mult($x);
return $s;
}
/*
$a[0][1] = $c[0] = new Fract(1);
for ($n=1; $n<16; $n++) { # compute coefficiehts for sum(i**n)
$pF = 1;
$s = new Fract(0);
for ($i=2; $i<= $n+1; $i++) {
$pF *= $i;
$s = ($i % 2) ? $s->sub($c[$n+1-$i]->div($pF)) : $s->add($c[$n+1-$i]->div($pF));
}
$c[$n] = $s;
$pG = new Fract(1, $n+1);
for ($i=$n+1; $i>=1 ; $i--) {
$a[$n][$i] = $c[$n+1-$i]->mult($pG);
$pG = $pG->mult($i);
}
krsort($a[$n]);
out("n $n, s $s, c", $c, ", a[$n]", array_map(fn ($v) => "$v", $a[$n]));
for ($x=0; $x<15; $x++) {
$s = new Fract(0);
for ($e=$n+1;$e>0;$e--)
$s = $s->add($a[$n][$e])->mult($x);
$q = 0;
for ($i=1; $i<=$x; $i++)
$q += $i**$n;
out("sum x**$n to $x = $s check $q", $s->eq($q) ? '' : 'differ ****************');
}
}
*/
try {
for ($e=0; $e<20; $e++)
for ($x=0; $x<8; $x++) {
$p = powersum($e, $x);
$q = 0;
for ($i=1; $i<=$x; $i++)
$q += $i**$e;
$p->eq($q) ? out("x**$e, 1..$x = $p") : err("x**$e, 1..$x = $p, check $q differs");
}
} catch (Throwable $t) {
out("catch $t");
}
outEnd(__file__);
?>