f03powersum.php

End f03powersum.php

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__);
?>