4

Is there an algorithm to find a representation of a given non-negative number $n$ as a sum of $k$, $p$ th powers where $k, p > 1$ ? That is a representation like,

$n = n_{1}^p + n_{2}^p + . . . + n_{k}^p$

Mathematica's, PowersRepresentations[n, k, p] solves this problem and I want to implement an algorithm myself to solve the problem.

I found that Waring's problem is closely related to this one but that's not exactly what I am looking for.

Prism
  • 11,162
thilinarmtb
  • 357
  • Mathematica is probably doing it the best way possible, but in the "Possible Issues" section from your link, they note that the computational time increases very quickly as the target number $n$ increases. So it seems there are no silver bullets to solve this problem quickly, unless perhaps if you make assumptions about $p$ and $k$ being small (and even then you might also need $n$ to be not too big). – user2566092 Sep 05 '13 at 20:30

3 Answers3

2

I am sending a dynamic programming implementation that solves this problem. This algorithm is pseudo-polynomial i.e. polynomial in the value of ${n}$ times $k$. This is exponential in the number of bits. With this kind of time complexity one would usually decide not to even bother, I decided to try it anyway as the Mathematica implementation also seems to be in this class or they would not have added the warning to the documentation that was pointed out above. My Perl implementation is faster than what Mathematica offers.

Here is some sample output:

$ ./powrep.pl 14000 3 2
12^2 + 20^2 + 116^2 = 14000
20^2 + 44^2 + 108^2 = 14000
20^2 + 60^2 + 100^2 = 14000
36^2 + 52^2 + 100^2 = 14000
44^2 + 60^2 + 92^2 = 14000
60^2 + 68^2 + 76^2 = 14000

$ ./powrep.pl 1729 2 3
1^3 + 12^3 = 1729
9^3 + 10^3 = 1729

$ ./powrep.pl 126  4 2
0^2 + 1^2 + 2^2 + 11^2 = 126
0^2 + 1^2 + 5^2 + 10^2 = 126
0^2 + 3^2 + 6^2 + 9^2 = 126
1^2 + 3^2 + 4^2 + 10^2 = 126
1^2 + 5^2 + 6^2 + 8^2 = 126
2^2 + 3^2 + 7^2 + 8^2 = 126
2^2 + 4^2 + 5^2 + 9^2 = 126
4^2 + 5^2 + 6^2 + 7^2 = 126

$ ./powrep.pl 28732 4 5
4^5 + 5^5 + 6^5 + 7^5 = 28732

And this is the code, written with brevity and simplicity in mind.

#! /usr/bin/perl -w

my %memo;
my @basepows;

sub powrep {
    my ($n, $k, $p) = @_;

    my $key = "$n-$k";
    return $memo{$key} if exists $memo{$key};

    return [] if $k == 1;

    my @res;
    foreach my $bpair (@basepows){
      my $bp = $bpair->[1];
      last if $bp >= $n;

      foreach my $r (@{powrep($n-$bp, $k-1, $p)}){
          push @res, [$bpair->[0], @$r]
            if $bpair->[0] <= $r->[0];
      }
    }

    $memo{$key} = \@res;
    return $memo{$key};
}


MAIN: {
    my $n = shift || 1729;
    my $k = shift || 2;
    my $p = shift || 3;

    die "all parameters (n, k, p) positive please"
      if $n<1 || $k<1 || $p<1;

    my ($q, $val) = (1);
    do {
      $val = $q ** $p;

      if($val<=$n){
          $memo{"$val-1"} = [[$q]];
          push @basepows, [$q, $val];
      }

      $q++;
    } while($val<$n);

    my @allsols;

    for(my $zterms = $k-1; $zterms >= 0; $zterms--){
      foreach my $sol (@{powrep($n, $k-$zterms, $p)}){
          push @allsols, 
          [(0) x $zterms,  @$sol];
      }
    }

    foreach my $sol (@allsols){
      my $check = 0;
      foreach my $t (@$sol){
          $check += $t ** $p;
      }

      my @terms = map { "$_^$p" } @$sol;
      print join(' + ', @terms);
      print " = $check\n";
    }

    exit 0;
}
Marko Riedel
  • 61,317
  • The memoization table for this problem has two dimensions that are mapped to an array of solutions, namely the integer $n$ that we are trying to represent and the number $k$ of the number of terms. Recursion reduces the number of terms by one in each step. – Marko Riedel Sep 05 '13 at 23:50
  • Interestingly enough the calculation with parameters as follows ./powrep.pl 405056 8 4 takes 11 seconds with the Perl program, whereas Mathematica takes 169 seconds (almost 3 minutes). – Marko Riedel Sep 06 '13 at 01:58
  • With parameters as follows ./powrep.pl 7065344 8 5 Perl takes 20 seconds and Mathematica 403 seconds (almost 7 minutes). The difference is surprising. – Marko Riedel Sep 06 '13 at 02:08
  • How large can $n$ be in your programme before a memory error is given? – thilinarmtb Sep 09 '13 at 04:55
  • I would guess that Perl integers will overflow before we run out of memory. – Marko Riedel Sep 09 '13 at 20:56
2

As there was some concern that the program would use too much memory, which was justified, I am sending a modified version that only stores and finds a single solution per n-k combination. This should make it possible to attack cases up to the limit imposed by the size of a Perl integer. 

#! /usr/bin/perl -w

my %memo;
my @basepows;

sub powrep {
    my ($n, $k, $p) = @_;

    my $key = "$n-$k";
    return $memo{$key} if exists $memo{$key};

    return [] if $k == 1;

    my @res = ();
    foreach my $bpair (@basepows){
      my $bp = $bpair->[1];
      last if $bp >= $n;

      my $resref = powrep($n-$bp, $k-1, $p);
      if(scalar(@$resref)){
          push @res, [$bpair->[0], @{$resref->[0]}];
        last;
      }
    }

    $memo{$key} = \@res;
    return $memo{$key};
}


MAIN: {
    my $n = shift || 1729;
    my $k = shift || 2;
    my $p = shift || 3;

    die "all parameters (n, k, p) positive please"
      if $n<1 || $k<1 || $p<1;

    my ($q, $val) = (1);
    do {
      $val = $q ** $p;

      if($val<=$n){
          $memo{"$val-1"} = [[$q]];
          push @basepows, [$q, $val];
      }

      $q++;
    } while($val<$n);

    my @allsols;

    for(my $zterms = $k-1; $zterms >= 0; $zterms--){
      my $resref = powrep($n, $k-$zterms, $p);
      if(scalar(@$resref)){
          push @allsols, [(0) x $zterms,  @{$resref->[0]}];
      }
    }

    foreach my $sol (@allsols){
      my $check = 0;
      foreach my $t (@$sol){
          $check += $t ** $p;
      }

      my @terms = map { "$_^$p" } @$sol;
      print join(' + ', @terms);
      print " = $check\n";
    }

    exit 0;
}

This modified program was able to compute the following table using half an hour of computation time and a very modest amout of memory.

0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 10000^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 2800^2 + 9600^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 192^2 + 1680^2 + 9856^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 8^2 + 40^2 + 1544^2 + 9880^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 1^2 + 1^2 + 1^2 + 1346^2 + 9909^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 1^2 + 1^2 + 1^2 + 5^2 + 3504^2 + 9366^2 = 100000000
0^2 + 0^2 + 0^2 + 1^2 + 1^2 + 1^2 + 2^2 + 2^2 + 5335^2 + 8458^2 = 100000000
0^2 + 0^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 4113^2 + 9115^2 = 100000000
0^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 9^2 + 3512^2 + 9363^2 = 100000000
1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 + 2^2 + 5335^2 + 8458^2 = 100000000
Marko Riedel
  • 61,317
1

Once we stop requiring that the algorithm produce all solutions, and demand instead that it produce one solution, we can use a greedy technique to speed up the algorithm even more. Allocate terms starting with the largest one possible. This gives a very fast algorithm that nonetheless does not miss any solutions, in the sense that if there is at least one, one solution will be found.

Here is one example.

0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 10000^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 9600^2 + 2800^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 9984^2 + 512^2 + 240^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 9992^2 + 392^2 + 56^2 + 56^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 0^2 + 9999^2 + 141^2 + 10^2 + 3^2 + 3^2 = 100000000
0^2 + 0^2 + 0^2 + 0^2 + 9999^2 + 141^2 + 10^2 + 4^2 + 1^2 + 1^2 = 100000000
0^2 + 0^2 + 0^2 + 9999^2 + 141^2 + 10^2 + 3^2 + 2^2 + 2^2 + 1^2 = 100000000
0^2 + 0^2 + 9999^2 + 141^2 + 9^2 + 5^2 + 3^2 + 1^2 + 1^2 + 1^2 = 100000000
0^2 + 9999^2 + 141^2 + 10^2 + 2^2 + 2^2 + 2^2 + 2^2 + 1^2 + 1^2 = 100000000
9999^2 + 141^2 + 10^2 + 3^2 + 2^2 + 1^2 + 1^2 + 1^2 + 1^2 + 1^2 = 100000000

And here is another.

0^3 + 0^3 + 0^3 + 0^3 + 396^3 + 312^3 + 196^3 = 100000000
0^3 + 0^3 + 0^3 + 463^3 + 79^3 + 57^3 + 41^3 = 100000000
0^3 + 0^3 + 464^3 + 46^3 + 17^3 + 7^3 + 4^3 = 100000000
0^3 + 464^3 + 46^3 + 16^3 + 10^3 + 6^3 + 2^3 = 100000000
464^3 + 46^3 + 16^3 + 9^3 + 7^3 + 5^3 + 3^3 = 100000000

One last example.

0^5 + 39^5 + 22^5 + 21^5 + 14^5 + 3^5 + 1^5 = 100000000
37^5 + 27^5 + 27^5 + 18^5 + 8^5 + 8^5 + 5^5 = 100000000

This is the code.

#! /usr/bin/perl -w

my %memo;
my @basepows;

sub powrep {
    my ($n, $k, $p) = @_;

    my $key = "$n-$k";
    return $memo{$key} if exists $memo{$key};

    return [] if $k == 1;

    my @res = ();
    foreach my $bpair (@basepows){
      my $bp = $bpair->[1];
      next if $bp >= $n;

      my $resref = powrep($n-$bp, $k-1, $p);
      if(scalar(@$resref)){
          push @res, [$bpair->[0], @{$resref->[0]}];
        last;
      }
    }

    $memo{$key} = \@res;
    return $memo{$key};
}


MAIN: {
    my $n = shift || 1729;
    my $k = shift || 2;
    my $p = shift || 3;

    die "all parameters (n, k, p) positive please"
      if $n<1 || $k<1 || $p<1;

    my ($q, $val) = (1);
    do {
      $val = $q ** $p;

      if($val<=$n){
          $memo{"$val-1"} = [[$q]];
          push @basepows, [$q, $val];
      }

      $q++;
    } while($val<$n);
    @basepows = reverse(@basepows);

    my @allsols;

    for(my $zterms = $k-1; $zterms >= 0; $zterms--){
      my $resref = powrep($n, $k-$zterms, $p);
      if(scalar(@$resref)){
          push @allsols, [(0) x $zterms,  @{$resref->[0]}];
      }
    }

    foreach my $sol (@allsols){
      my $check = 0;
      foreach my $t (@$sol){
          $check += $t ** $p;
      }

      my @terms = map { "$_^$p" } @$sol;
      print join(' + ', @terms);
      print " = $check\n";
    }

    exit 0;
}
Marko Riedel
  • 61,317