# Fermat's polygonal number theorem

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}$$

## 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 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
53, 10 => 10, 15, 27
879, 17 => 17, 48, 155, 231, 428
4856, 23 => 130, 448, 955, 1398, 1925

• Sandbox post – caird coinheringaahing Oct 7 at 17:14
• Can the output have some zero padding? For instance if we consider x=17, s=5 could we output 5,12,0,0,0 instead of just 5,12? – flawr Oct 7 at 17:30
• @flawr So long as the length of the array doesn't exceed $s$, even with padding, that's fine – caird coinheringaahing Oct 7 at 17:42
• Are repeats allowed or should I add a Q to my submission? – Jonathan Allan Oct 7 at 18:46
• @JonathanAllan Repeated outputs are perfectly fine (if outputting multiple solutions) – caird coinheringaahing Oct 7 at 18:47

# 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)div2+n*(2-n)))([0..n]<[1..s]),sum x==n]!!0  Try it online! • @JonathanAllan I fixed it now - I misunderstood something. (It even saved me one byte:) – flawr Oct 7 at 18:52 # 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  • @AZTECCO I may try to fix it later. Removed for now. – Arnauld Oct 15 at 11:39 • Thanks. Waiting for it! – AZTECCO Oct 15 at 11:59 # 05AB1E, 14 13 bytes Port of Jonathan Allan's Jelly answer. ÍIÝ*>.¥¹ã.ΔOQ  Try it online! # Jelly, 17 bytes x’2;’ÄÄx⁸ŒPS⁼¥Ƈ⁹Ḣ  A (very very inefficient) dyadic Link accepting s on the left and x on the right which yields the shortest possible answer as a list of integers (sorted ascending). 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], ...] Ḣ - head [5,12]  • @AZTECCO That's perfectly fine, it times out on TIO there at 60 seconds (I'm pretty sure even far smaller input number than that will time out). As I point out in my answer, this is "very very inefficient" and that there's "not much point trying it for much higher values!". Remember, the code given for a code-golf solution only need work given infinite resources. – Jonathan Allan Oct 15 at 12:07 • ok I tested with s=3 and n=5 and it taken 12 seconds!! I like this inefficient solution and I'll trust you, even if it's almost impossible to test it :) thank you! – AZTECCO Oct 15 at 12:20 • @AZTECCO Indeed, it's hard to test - I originally did not have the ⁸ so was repeating the elements of the cumulative sum result $x$ times each rather than $s$ times each (which is incorrect) and it was hard to spot that it was actually doing anything wrong without running just parts of the code! – Jonathan Allan Oct 15 at 12:28 # 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! # 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, ...  # 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


# 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!

# 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