Fermat's polygonal number theorem

Question

Fermat's polygonal number theorem states that every positive integer can be expressed as the sum of at most \$n\$ \$n\$-gonal numbers. This means that every positive integer can be expressed as the sum of up to three triangle numbers, four square numbers, five pentagonal numbers etc. Your task is to take a positive integer \$x\$, and an integer \$s \ge 3\$, and to output the \$s\$-gonal integers which sum to \$x\$.

The \$n\$-th \$s\$-gonal integer, where \$n \ge 1\$ and \$s \ge 3\$, can be defined in a couple of ways. The non-math-y way is that the \$n\$th \$s\$-gonal number can be constructed as a regular polygon with \$s\$ sides, each of length \$n\$. For example, for \$s = 3\$ (triangular numbers):

See here for examples with a larger \$s\$.

The math-y definition is by using the formula for \$P(n, s)\$, which yields the \$n\$-th \$s\$-gonal number:

$$P(n, s) = \frac{n^2(s-2) - n(s-4)}{2}$$

which is given in the Wikipedia page here.

Input

Two positive integers, \$s\$ and \$x\$, with the condition \$s \ge 3\$. You may input these integers in the most natural representation in your language (decimal, unary, Church numerals, integer-valued floating point numbers etc.).

Output

A list of integers, \$L\$, with a maximum length of \$s\$, where the sum of \$L\$ is equal to \$x\$ and all integers in \$L\$ are \$s\$-gonal integers. Again, the integers may be outputted in the natural representation in your language, with any distinct, consistent separator (so non-decimal character(s) for decimal output, a character different from that used for unary output etc.)

Rules

• The inputs or outputs will never exceed the integer limit for your language
• \$L\$ does not have to be ordered
• In case of multiple possible outputs, any or all are acceptable
• This is code-golf so the shortest code in bytes wins

Test cases

   x,  s => L
   1,  s => 1
   2,  s => 1, 1
   5,  6 => 1, 1, 1, 1, 1
  17,  3 => 1, 6, 10
  17,  4 => 1, 16
  17,  5 => 5, 12
  36,  3 => 36
  43,  6 => 15, 28
 879, 17 => 17, 48, 155, 231, 428
4856, 23 => 130, 448, 955, 1398, 1925

Answer 1

Haskell, 55 bytes

n%s=[l|l<-mapM(\_->scanl(+)0[1,s-1..n])[1..s],sum l==n]

Try it online!

Outputs all possible solutions. Defines the s-gonal numbers as the cumulative sum of the arithmetic progression

1, s-2, 2*s-3, 3*s-4, ...

Answer 2

Haskell, 78 80 77 bytes

We compute the cartesian product of the first \$n\$ s-gonal numbers, and then find the first entry that sums to \$n\$.

s#n=[x|x<-mapM(map(\n->s*(n^2-n)`div`2+n*(2-n)))([0..n]<$[1..s]),sum x==n]!!0

Try it online!

Answer 3

JavaScript (ES6),  83  80 bytes

A fast recursive search which maximizes the smallest term of the output.

Takes input as (s)(x).

s=>g=(x,n=0,a=[],y=~n*(~-n-n*s/2))=>x<y?x|a[s]?0:a:g(x,n+1,a)||g(x-y,n,[...a,y])

Try it online!

Formula

It turns out to be shorter to use a 0-based formula to compute the \$s\$-gonal numbers in JS, i.e. to start with \$n=0\$ and to compute \$P(n+1,s)\$:

$$\begin{align}P(n+1,s)&=((n+1)^2(s-2)-(n+1)(s-4))/2\\ &=(n^2(s-2)+ns+2)/2\\ &=-(n+1)((n-1)-ns/2) \end{align}$$

which can be written in 14 bytes:

~n*(~-n-n*s/2)

Commented

s =>                         // main function taking s
  g = (                      // recursive function g
    x,                       // taking x
    n = 0,                   // start with n = 0
    a = [],                  // a[] = list of s-gonal numbers
    y =                      // y = P(n + 1, s)
      ~n * (~-n - n * s / 2) //   = -(n + 1) * ((n - 1) - n * s / 2)
  ) =>                       //
    x < y ?                  // if x is less than P(n + 1, s):
      x | a[s] ?             //   if x is not equal to 0 or a[] is too long:
        0                    //     failed: return 0
      :                      //   else:
        a                    //     success: return a[]
    :                        // else:
                             //   process recursive calls:
      g(x, n + 1, a) ||      //   preferred: try to increment n
      g(x - y, n, [...a, y]) //   fallback : try to use the current s-gonal number

Answer 4

05AB1E, 14 13 bytes

Port of Jonathan Allan's Jelly answer.

ÍIÝ*>.¥¹ã.ΔOQ

Try it online!

Answer 5

Ruby, 79 bytes

Computes the first \$n\$ \$s\$-gonal numbers and zero, gets the cartesian product of itself \$s\$ times, and finds the first entry that matches. Later test cases run out of memory on TIO but theoretically work.

\$\frac{n^2(s-2)-n(s-4)}{2}\$ is reorganized into \$\frac{n(ns-2n-s+4)}{2}\$ because it saves bytes in Ruby.

->n,s{a=(0..n).map{|i|i*(i*s-2*i-s+4)/2};a.product(*[a]*~-s).find{|a|a.sum==n}}

Try it online!

Answer 6

Retina, 111 bytes

\d+
*
~(`$
$"
0%["^_+ "|""]'$L$`\G_(?<=(?=___(_*))_+)
((_(?($.(2*$>`))$1\$.(2*$>`)))$*)
1%|' L$`\G_
$$.$.($`$>`

Try it online! Link includes test cases. Takes input in the order s n. Explanation:

\d+
*

Convert to unary.

~(`

After processing the remaining stages, treat them as a Retina program and execute them on the same input.

$
$"

Duplicate the line.

0%["^_+ "|""]'$L$`\G_(?<=(?=___(_*))_+)
((_(?($.(2*$>`))$1\$.(2*$>`)))$*)

Replace the first copy with a regular expression that skips over the first number and then matches s s-gonal numbers. The numbers themselves are captured in odd capture groups and the even capture groups are used to ensure that all of the numbers are s-gonal.

1%|' L$`\G_
$$.$.($`$>`

Replace the second copy with a space-separated list of the odd capture groups.

As an example, the generated code for an input of 5 17 is as follows:

^_+ ((_(?(2)__\2))*)((_(?(4)__\4))*)((_(?(6)__\6))*)((_(?(8)__\8))*)((_(?(10)__\10))*)$
$.1 $.3 $.5 $.7 $.9

Answer 7

Python 3, 105 bytes

def f(s,x,n=0,a=[]):y=~n*(~-n-n*s/2);return(x<1)*a*(len(a)<=s)if x<y else f(s,x,n+1,a)or f(s,x-y,n,a+[y])

Try it online!

Port of Arnauld's JavaScript answer, so make sure to upvote him!

Answer 8

Jelly,  17  15 bytes

-2 thanks to caird coinheringaahing

x’2;’ÄÄx⁸ŒPS⁼¥Ƈ

A (very very inefficient) dyadic Link accepting s on the left and x on the right which yields a list of many possible L lists, with repeats.

Try it online! - not much point trying it for much higher values!

How?

x’2;’ÄÄx⁸ŒPS⁼¥Ƈ - Link: s, x                    e.g.  5, 17
x               - repeat (s) (x) times                [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]
 ’              - decrement (vectorises)              [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]
  2;            - prepend a two                       [2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]
    ’           - decrement (vectorises)              [1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
     Ä          - cumulative sums                     [1, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31, 34, 37, 40, 43, 46, 49, 52]
      Ä         - cumulative sums                     [1, 5, 12, 22, 35, 51, 70, 92, 117, 145, 176, 210, 247, 287, 330, 376, 425, 477]
       x⁸       - repeat (each of those) (s) times    [1, 1, 1, 5, ..., 425, 477, 477, 477]
         ŒP     - power-set                           [[], [1], [1], ..., [1, 1], ..., [5, 22, 70], ... etc]
                    (this has 2^(x(s+1)) entries ...this example would have 2^(17(5+1)) = 2^102 = 5070602400912917605986812821504 entries!)
                    (Note: the lengths increase left to right)
              Ƈ - filter keep if:
             ¥  -   last two links as a dyad:
           S    -     sum
            ⁼   -     equals (x)?                     [[5,12], ... , [5,12], [1, 1, 5, 5, 5], ... , [1, 1, 5, 5, 5], [1, 1, 1, 1, 1, 12], ...]

Answer 9

C++ (clang), 198 bytes

#import<vector>
using V=std::vector<int>;V f(int n,int s){V _{0};int z=1,a=0,b=1,i,t;for(;a<n;b+=s-2)_.push_back(a+=b),++z;V o;for(t=a=0;t-n;b=++a)for(o=V(s),t=i=0;b;b/=z)t+=o[i++]=_[b%z];return o;}

Try it online!

V=vector<int> 
V _{0}; // initialized with one element =0 
int z=1, // vector size 
a=0,b=1,i,t;for(;a<n;b+=s-2)_.push_back(a+=b),++z;
// pushes polygons in V
V o; // vector to be returned 
for(t=a=0;t-n;b=++a) // ends when t=n
// loop to generate multi-dimension indexes
// for example a=1234 z=10
// a%z->4 , a/=z , a%z-> 3 , ... 2 , 1
for(o=V(s),t=i=0;b;b/=z)// loop to extract indexes
t+=o[i++]=_[b%z]; // put the sum in t and values in o
return o

Answer 10

APL(NARS), 149 chars, 298 bytes

r←f w;n;s;i;k
(n s)←w⋄r←⍬⋄→0×⍳s<3⋄i←1
→0×⍳n<k←2÷⍨(i×i×s-2)-i×s-4⋄r←r,k⋄i+←1⋄→2

h←{0=≢b←((v←↑⍵)=+/¨a)/a←{0=≢⍵:⊂⍬⋄m,(⊂1⌷⍵),¨m←∇1↓⍵}f⍵:v⍴1⋄k←↑⍋≢¨b⋄k⊃b}

if not find solutions "0=≢b" than return for (n s) input, n times 1; else it would return the sum of s numbers that has less addend...

test:

  h 1 3
1 
  h 2 8
1 1 
  h 5 6
1 1 1 1 1 
  h 17 3
1 6 10 
  h 17 4
1 16 
  h 17 5
5 12 
  h 36 3
36 
  h 43 6
15 28 
  h 879 17
17 48 155 231 428 
  h 4856 23
321 448 596 955 2536 
  +/h 4856 23
4856

The problem of this: It not find some solution has some number repeat in the sum...