Problem 382: Generating Polygons

View on Project Euler

Project Euler Problem 382 Solution

EulerSolve provides an optimized solution for Project Euler Problem 382, Generating Polygons, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The side lengths are generated by $$s_1=1,\quad s_2=2,\quad s_3=3,\quad s_n=s_{n-1}+s_{n-3}\quad (n\ge 4).$$ Because \(s_{n-3}\gt 0\), the sequence is strictly increasing, so every subset of \(\{s_1,\dots,s_n\}\) has a unique largest selected side. Let \(f(n)\) be the number of nonempty subsets that can be arranged into a nondegenerate polygon. The solver computes \(f(10^{18}) \bmod 10^9\). Mathematical Approach Step 1: Reduce Polygon Feasibility to One Inequality For any finite set of positive side lengths, a polygon exists if and only if the largest side is smaller than the sum of all remaining sides. If a subset has largest side \(s_n\), then it is bad exactly when the other chosen sides have total length at most \(s_n\). This immediately turns the geometric problem into a counting problem about subset sums. Step 2: Count Bad Subsets by Their Largest Index Define \(b_n\) as the number of subsets of \(\{s_1,\dots,s_{n-1}\}\) whose sum is at most \(s_n\): $$b_n=\#\left\{A\subseteq \{1,\dots,n-1\}:\sum_{i\in A}s_i\le s_n\right\}.$$ Every non-polygon subset of \(\{s_1,\dots,s_n\}\) has a unique largest element \(s_k\), and after removing that largest side the remaining subset contributes to \(b_k\)....

Detailed mathematical approach

Problem Summary

A polygon is a flat shape built from straight segments joined into a closed circuit; it has at least three sides and does not cross itself. A set \(S\) of positive numbers generates a polygon \(P\) when no two sides of \(P\) have the same length, every side length of \(P\) belongs to \(S\), and \(S\) contains nothing else — so the side lengths of \(P\) are exactly the elements of \(S\), each used once. [1]

The available lengths come from one sequence,

$$s_1=1,\quad s_2=2,\quad s_3=3,$$

$$s_n=s_{n-1}+s_{n-3}\quad(n\gt 3),$$

whose first terms are \(1,2,3,4,6,9,13,19,28,41\). With \(U_n=\{s_1,\dots,s_n\}\) — so \(U_{10}\) is that list — let \(f(n)\) count the subsets of \(U_n\) that generate at least one polygon. The statement supplies \(f(5)=7\), \(f(10)=501\) and \(f(25)=18635853\), and asks for the last nine digits of \(f(10^{18})\), that is \(f(10^{18})\bmod 10^9\).

The whole solution in one sentence: a set of distinct lengths closes into a polygon exactly when its largest element is smaller than the sum of the others, so we count the subsets that fail this test instead, show that those failure counts obey a short linear recurrence, and then leap \(10^{18}\) steps of that recurrence with a single matrix power modulo \(10^9\).

Mathematical Approach

1. When do given lengths close into a polygon?

Claim. Let \(a_1,\dots,a_k\) be positive numbers with sum \(\sigma\) and largest value \(m\). A polygon whose side lengths are exactly these numbers exists if and only if \(2m\lt\sigma\).

The condition already forces \(k\ge3\): for \(k=1\) it reads \(2m\lt m\), and for \(k=2\) the other length is at most \(m\), so \(\sigma\le2m\). The side count needs no separate test.

Necessity. Let the polygon have vertices \(V_0,\dots,V_{k-1}\) in order, the longest side joining \(V_{k-1}\) to \(V_0\). The other sides form a path of length \(\sigma-m\) between those two points, and no path is shorter than the segment joining its endpoints, so \(\sigma-m\ge m\). [2] Equality would force every intermediate vertex onto that segment, so the long side would overlap the path and the circuit would cross itself. Hence \(\sigma-m\gt m\).

Sufficiency. We build a convex polygon with every vertex on one circle — a cyclic polygon. On a circle of radius \(R\ge m/2\) a chord of length \(a_j\) subtends the central angle \(\theta_j(R)=2\arcsin(a_j/2R)\), as bisecting the isosceles triangle of that chord and two radii shows; each \(\theta_j\) is positive, continuous and decreasing in \(R\), and a longest side gives \(\theta_1(m/2)=\pi\). Chords laid end to end round the circle close up exactly when their arcs total \(2\pi\). If \(\theta_1+\dots+\theta_k\ge2\pi\) at \(R=m/2\), that sum tends to \(0\) for large \(R\), so some radius makes it exactly \(2\pi\). Otherwise \(\theta_2+\dots+\theta_k\lt\theta_1\) there, so \(D=\theta_2+\dots+\theta_k-\theta_1\) is negative at \(R=m/2\), while \(\arcsin x=x+O(x^3)\) gives

$$\lim_{R\to\infty}R\,D(R)=\sigma-2m\gt0,$$

so \(D\) vanishes at some radius. There give the longest chord the major arc \(2\pi-\theta_1\) and each other chord its arc \(\theta_j\); these again total \(2\pi\), and a chord has the same length on either of its arcs. All arcs are positive, so in both cases the \(k\) vertices are distinct and in circular order, and joining consecutive ones gives a convex — hence simple — polygon with the required sides.

The statement's four examples are the two cases of this claim: \(10\lt12\) for \(\{3,4,5\}\) and \(48\lt50\) for \(\{6,9,11,24\}\), against \(6=6\) for \(\{1,2,3\}\) and \(18=18\) for \(\{2,3,4,9\}\), both sitting exactly on the boundary — so the inequality must indeed be strict. More is true: for each prescribed cyclic order of the lengths the same condition gives a unique convex cyclic polygon, up to isometry. [3]

2. Counting the failures instead

Since \(s_{n-3}\gt0\) the sequence is strictly increasing, so its terms are distinct, every non-empty subset of \(U_n\) has a unique largest element, and the requirement that no two sides agree is automatic.

Let \(S\subseteq U_n\) be non-empty with largest element \(s_k\). By Section 1 it generates a polygon exactly when \(s_k\) is below the total of the other chosen elements, so it fails exactly when they total at most \(s_k\). Deleting \(s_k\) is a bijection between the failing subsets with largest element \(s_k\) and the sets counted by

$$b_k=\#\{A\subseteq U_{k-1}:\sigma(A)\le s_k\},$$

where \(\sigma(A)\) denotes the sum of \(A\) and \(U_0=\emptyset\). Adding the empty set, the number of subsets of \(U_n\) generating no polygon is \(1+b_1+\dots+b_n\), and therefore

$$\begin{aligned}f(n)&=2^n-1-B_n,\\B_n&=\sum_{k=1}^{n}b_k.\end{aligned}$$

The geometry is gone; what remains is a question about subset sums.

3. The sequence and its running total

The sequence is Narayana's cows sequence re-indexed, \(s_n=a(n+1)\) in OEIS A000930, and consecutive ratios tend to \(1.465571\ldots\), the supergolden ratio, the real root of \(x^3=x^2+1\). [4] [5] So the limits defining \(b_k\) grow exponentially and no table indexed by the sum can be built.

One identity makes everything work. Writing \(P_k=s_1+\dots+s_k\),

$$P_k=s_{k+3}-3 .$$

Read backwards the rule says \(s_{n-3}=s_n-s_{n-1}\) for \(n\gt3\); summing for \(n=4,\dots,k+3\) telescopes and leaves \(P_k=s_{k+3}-s_3\). For example \(P_1=1=s_4-3\) and \(P_{10}=126=s_{13}-3\).

The identity cuts both ways. The whole of \(U_{n-1}\) totals \(s_{n+2}-3\), far more than \(s_n\) — for \(n=10\), \(85\) against \(41\) — so the restriction \(\sigma(A)\le s_n\) really bites and \(b_n\) stays well below \(2^{n-1}\). Conversely a limit of \(s_{k+3}\) or more restricts subsets of \(U_k\) not at all, which is what stops the recursion of Section 5.

4. A worked example: \(n=5\)

Here \(U_5=\{1,2,3,4,6\}\); list the failing subsets by their largest element.

  • Largest \(1\): only the empty remainder, so \(b_1=1\).
  • Largest \(2\): remainders \(\emptyset,\{1\}\), so \(b_2=2\).
  • Largest \(3\): all four subsets of \(\{1,2\}\) have sum at most \(3\), so \(b_3=4\).
  • Largest \(4\): of the eight subsets of \(\{1,2,3\}\), the two of sum \(5\) and \(6\) are too heavy, so \(b_4=6\).
  • Largest \(6\): of the sixteen subsets of \(\{1,2,3,4\}\), the five with sums \(7,7,8,9,10\) are too heavy, so \(b_5=11\).

Hence \(B_5=1+2+4+6+11=24\) and \(f(5)=2^5-1-24=7\), matching the statement. The seven winners are \(\{2,3,4\}\), \(\{3,4,6\}\), \(\{1,2,3,4\}\), \(\{1,2,4,6\}\), \(\{1,3,4,6\}\), \(\{2,3,4,6\}\) and \(\{1,2,3,4,6\}\). Continuing the same count gives

$$\begin{aligned}(b_1,\dots,b_8)&=(1,2,4,6,\\&\quad 11,20,36,67),\end{aligned}$$

with \(B_8=147\) — exactly the constants the three solution files start from.

5. Peeling off the largest index

Generalise \(b_k\) slightly: for \(i\ge0\) and \(L\ge0\) put

$$C(i,L)=\#\{A\subseteq U_i:\sigma(A)\le L\},$$

so that \(b_n=C(n-1,s_n)\). Three rules evaluate it. First, \(C(0,L)=1\), the empty set alone. Second, if \(L\ge P_i\) then every subset qualifies and \(C(i,L)=2^i\); by Section 3 this happens as soon as \(L\ge s_{i+3}-3\). Third, otherwise split on whether the largest available term \(s_i\) is used:

$$\begin{aligned}C(i,L)&=C(i-1,L)\\&\quad+C(i-1,L-s_i)\end{aligned}$$

when \(s_i\le L\), and \(C(i,L)=C(i-1,L)\) when \(s_i\gt L\). This is the usual counting form of the subset-sum recursion, and alone it is useless here: the limits have size \(s_n\).

6. Only seven limit forms survive

The saving observation: start from \(C(n-1,s_n)\) and apply the third rule repeatedly. A call whose limit reaches \(P_i\) is settled at once by the second rule and may carry any limit, but every call that survives both stopping tests at a level \(i\ge6\) has, against its own level, one of seven limits:

$$0,\ s_{i-1},\ s_i,\ s_{i+1},\ s_{i+2},$$

$$s_{i+1}+s_{i-4},\qquad s_{i+2}+s_{i-3}.$$

Give the six non-zero forms names, defined for \(i\ge5\):

  • \(p_i=C(i,s_{i-1})\), \(q_i=C(i,s_i)\), \(r_i=C(i,s_{i+1})\);
  • \(t_i=C(i,s_{i+1}+s_{i-4})\), \(u_i=C(i,s_{i+2})\), \(v_i=C(i,s_{i+2}+s_{i-3})\).

The quantity we want is \(b_{i+1}=r_i\). Four immediate consequences of the defining rule are needed: \(s_{i+1}-s_i=s_{i-2}\); \(s_{i-2}+s_{i-4}=s_{i-1}\); \(s_{i-1}+s_{i-2}=s_i+s_{i-5}\); and \(s_{i-1}+s_{i-2}+s_{i-3}=s_{i+1}\). With \(P_{i-1}=s_{i+2}-3\) they turn one application of the third rule into six exact equations:

  • \(p_i=q_{i-1}\), because \(s_i\gt s_{i-1}\) so \(s_i\) is never used;
  • \(q_i=r_{i-1}+1\), because the leftover limit is \(0\);
  • \(r_i=u_{i-1}+p_{i-1}\), because \(s_{i+1}-s_i=s_{i-2}\);
  • \(t_i=v_{i-1}+q_{i-1}\), because \(s_{i-2}+s_{i-4}=s_{i-1}\);
  • \(u_i=2^{i-1}+t_{i-1}\), because \(s_{i+2}\ge P_{i-1}\) and \(s_{i+2}-s_i=s_i+s_{i-5}\);
  • \(v_i=2^{i-1}+u_{i-1}\), because \(s_{i+2}+s_{i-3}\ge P_{i-1}\) and the leftover limit is \(s_{i+1}\).

All six hold for \(i\ge6\). Each is one use of the split, the second rule of Section 5 settling any unconstrained branch at once, and a zero limit admits only the empty subset. So a surviving call of one of the seven forms has only children that stop or take one of the same forms a level lower; descending from \(r_{n-1}\), induction closes the list down to the levels below \(6\), evaluated directly.

Numerically, \(r_6=C(6,13)=b_7\). Dropping \(s_6=9\) leaves \(C(5,13)=29\), since only three of the \(32\) subsets of \(U_5\) exceed \(13\); keeping \(9\) leaves the limit \(4\) and \(C(5,4)=7\). So \(r_6=29+7=36\), the third equation at \(i=6\).

7. One recurrence for the failure counts

Now eliminate. The second equation gives \(q_i=r_{i-1}+1\), and with the first, \(p_i=r_{i-2}+1\). Substituting into the third,

$$r_i=r_{i-3}+u_{i-1}+1 .$$

Chaining the fifth, fourth and sixth equations in turn,

$$u_i=u_{i-3}+r_{i-3}+2^{i-1}+2^{i-3}+1 .$$

The first of these two says \(u_{i-1}=r_i-r_{i-3}-1\). Putting that into the second, taken at index \(i-1\), removes \(u\) altogether and leaves

$$\begin{aligned}r_i&=2r_{i-3}+r_{i-4}-r_{i-6}\\&\quad+5\cdot2^{i-4}+1,\end{aligned}$$

using \(2^{i-2}+2^{i-4}=5\cdot2^{i-4}\). Since \(b_n=r_{n-1}\), this is

$$\begin{aligned}b_n&=2b_{n-3}+b_{n-4}-b_{n-6}\\&\quad+5\cdot2^{n-5}+1\end{aligned}$$

Used only in their stated range, the six equations give the two displayed identities for \(i\ge8\), and eliminating \(u\) also needs the first at \(i-3\), so this derivation reaches only \(i\ge11\), that is \(n\ge12\). The five remaining indices are checked directly: for \(n=7,8,9,10,11\) Section 5 counts \(36,67,128,247,480\), and the right-hand side returns the same — at \(n=9\), \(2\cdot20+11-4+5\cdot16+1=128\). So the relation holds for every \(n\ge7\).

8. From inhomogeneous to homogeneous

The three solution files use not that relation but a homogeneous one of order \(8\), which follows from ours by two cancellations. Let \(S\) be the shift operator, \((Sb)_n=b_{n+1}\), and set

$$P(S)=S^6-2S^3-S^2+1 .$$

Shifting that recurrence so its top index is \(n+6\) turns it into

$$(P(S)b)_n=10\cdot2^{n}+1\qquad(n\ge1).$$

Applying \(S-2\) kills the power of two, because \(10\cdot2^{n+1}+1-2(10\cdot2^{n}+1)=-1\); applying \(S-1\) then kills the constant, because \(-1-(-1)=0\). So \(b\) is annihilated by \((S-1)(S-2)P(S)\), whose expansion is

$$\begin{aligned}&x^8-3x^7+2x^6-2x^5\\&\quad+5x^4-x^3-x^2-3x+2 .\end{aligned}$$

Read as a recurrence, that is

$$\begin{aligned}b_n&=3b_{n-1}-2b_{n-2}+2b_{n-3}-5b_{n-4}\\&\quad+b_{n-5}+b_{n-6}+3b_{n-7}-2b_{n-8}\end{aligned}$$

valid for all \(n\ge9\). The coefficient list \((3,-2,2,-5,1,1,3,-2)\) and the starting index \(9\) are exactly what the three programs hard-code, and the eight values of Section 4 are its initial data. This matters: the C++ program only tests the relation up to \(n=60\) and the other two test nothing, whereas the derivation above, with its boundary checks, proves it for every \(n\ge9\). The argument is the standard closure property that a sum of linearly recurrent sequences is again one, here applied to \(b\), \(2^n\) and the constants. [8]

9. What the characteristic roots say

The sixth-degree factor splits over the integers,

$$\begin{aligned}&x^6-2x^3-x^2+1\\&\quad=(x^3-x-1)(x^3+x-1),\end{aligned}$$

so the full characteristic polynomial is \((x-1)(x-2)\)\({}(x^3-x-1)(x^3+x-1)\). The first cubic is the minimal polynomial of the plastic ratio \(1.324717\ldots\) [6], and the second is the minimal polynomial of \(0.682327\ldots\), the reciprocal of the supergolden ratio that drives \(s\) itself. [5] Apart from \(2\) the moduli are \(1.324718\), \(1.210608\) twice, \(1\), \(0.868837\) twice and \(0.682328\).

The root \(2\) dominates. A particular solution \(b_n=A\cdot2^n+C\) of the order-six relation needs \(64A=16A+4A-A+10\) and \(C=2C+1\), so \(A=2/9\) and \(C=-1\); the difference obeys the homogeneous relation, whose roots are all below \(2\) in modulus, so \(b_n/2^n\to2/9\). Summing that leading term, \(B_n\) approaches \(\tfrac49 2^n\) and

$$\frac{f(n)}{2^n}\longrightarrow\frac59 .$$

Asymptotically five ninths of all subsets of \(U_n\) generate a polygon, and the convergence is fast: the ratio is \(0.5553917\) at \(n=25\) and \(0.5555556\) at \(n=100\).

10. Summing with the same matrix

We need \(B_n\), not \(b_n\), so the state carries one extra coordinate, taken as a column vector:

$$w_n=(b_n,b_{n-1},\dots,b_{n-7},B_n).$$

The \(9\times9\) transition \(T\) carries the eight coefficients in its first row, a shifted identity in rows two to eight, and in its last row those coefficients again plus a \(1\) in the last column — that row computes \(B_n=B_{n-1}+b_n\). Starting from

$$w_8=(67,36,20,11,6,4,2,1,147)$$

we get \(w_n=T^{\,n-8}w_8\) for \(n\ge8\). Binary exponentiation reads the exponent in binary: for \(13=8+4+1\), repeated squaring gives \(T^2,T^4,T^8\) and \(T^{13}=T^8T^4T\). That is at most two matrix products per bit, so \(O(\log n)\) in all. [7] All arithmetic is modulo \(10^9\), the negative coefficients \(-2\) and \(-5\) first becoming \(10^9-2\) and \(10^9-5\). Note \(10^9=2^9\cdot5^9\) is not prime, so most residues have no inverse; this never matters, as nothing divides.

11. The last nine digits

Finally,

$$f(n)\equiv 2^n-1-B_n\pmod{10^9},$$

with \(2^n\bmod 10^9\) from the same binary exponentiation and the difference shifted back into \(\{0,\dots,10^9-1\}\) if negative. For \(n\le8\) no matrix is needed: \(B_n\) is a partial sum of the eight base values. Recomputed from scratch the method reproduces \(f(5)=7\), \(f(10)=501\), \(f(25)=18635853\), and further out gives \(f(40)\equiv839417413\) and \(f(60)\equiv899724066\) modulo \(10^9\).

How the Code Works

The supplied C++ file holds both an exact validator and the fast solver. parse_arguments accepts --n= and --skip-checkpoints; by default \(n=10^{18}\) and the checks run. The modular helpers are mod_norm, which brings a signed 128-bit value into range and so absorbs the negative coefficients of Section 8, add_mod, which also detects 64-bit wraparound, and mul_mod and pow_mod, which multiply through __uint128_t.

ExactCounter implements Section 5. Its constructor fills s_ and the 128-bit prefix sums prefix_. count_subsets_leq(i, limit) returns \(1\) at \(i=0\), returns \(2^i\) once limit reaches prefix_[i] — the second rule of Section 5 — and otherwise splits on \(s_i\), caching in an unordered_map keyed by Key{i, limit}. Section 6 explains why that cache stays tiny.

run_checkpoints builds \(b_n\) exactly for \(n\le60\) through count_subsets_leq(n-1, s[n]), verifies the order-eight relation term by term from \(n=9\), rebuilds \(f(n)\) as \(2^n-1-B_n\), compares \(f(5)\), \(f(10)\) and \(f(25)\) with the published values, and checks the fast solver against exact arithmetic for every \(n\le60\). Any mismatch throws before anything is printed.

The fast path is prefix_sum_b_mod. It answers \(n\le8\) by summing kBase, and otherwise assembles the state and matrix of Section 10: kRecurrence normalised by mod_norm in rows \(0\) and \(8\), ones below the diagonal in rows \(1\) to \(7\), and trans[8][8] = 1. apply_mat_pow builds the power with mat_mul, which skips zero factors — useful, as \(T\) is mostly zeros — then applies it with mat_vec_mul and returns the advanced state, whose component \(8\) is \(B_n\). solve_mod combines that with pow_mod(2, n, kMod) as in Section 11, and main prints the value, also zero-padded to nine digits.

The Python and Java files carry the fast path only, with the same constants, transition and congruence and \(n=10^{18}\) written in directly; neither accepts arguments nor validates anything. Integer handling differs. C++ uses __uint128_t for the exact counts and modular products, and signed 128-bit words in mod_norm — headroom rather than necessity, since for \(n\le60\) every \(b_n\) is below \(2^{59}\), every exact \(f(n)\) below \(2^{60}\) and every modular product below \(10^{18}\), all inside signed 64 bits. Java stays within long, its largest product being below \((10^9)^2\), comfortably under \(2^{63}-1\), and its modNorm adds the modulus back on a negative difference, as C++ does. Python needs neither precaution: its integers are unbounded and % returns a non-negative remainder, so its mod_norm is a bare modulus — which would be wrong in a language with truncating division.

Complexity and Verification

The production path is one binary exponentiation of a fixed \(9\times9\) matrix: \(O(\log n)\) products of at most \(9^3=729\) multiply-and-reduce steps, roughly \(10^5\) modular operations at the requested index, plus \(O(\log n)\) for the scalar power of two, in \(O(1)\) memory. The C++ validator runs only to \(n=60\), and Section 6 bounds its cached pairs by \(O(n)\), so it too finishes at once. A separate arbitrary-precision run with the same rules reaches \(n=200\), caching \(1368\) pairs \((i,L)\) and visiting \(1765\) distinct ones with terminal calls counted; it is not the C++ validator, whose 64-bit sequence words would overflow at \(s_{117}\).

The C++ checks were described above. Further checks were run for this article, none of them part of the solution files.

  • A brute force over all \(2^n\) subsets of \(U_n\), testing each directly against the criterion of Section 1, reproduces \(f(n)\) for every \(n\le20\) and \(b_n\) for every \(n\le18\).
  • Instrumenting the recursion confirms that every surviving call at a level \(i\ge6\) carries one of the seven limit forms, and that all six peel equations hold for \(6\le i\le150\). Both recurrences agree with exact subset-sum counts for every \(n\le200\).
  • A second \(9\times9\) matrix, built from the proved order-six relation by carrying \(2^{\,n-4}\) and \(1\) as extra coordinates, returns the same \(B_n\bmod 10^9\) as the programs' matrix for all \(n\le300\) and at \(n=10^6\), \(10^{12}\) and \(10^{18}\); plain iteration of that relation modulo \(10^9\) agrees too, to \(n=200000\).
  • The C++, Python and Java programs, plus an independent reimplementation, produce the same nine digits, and the C++ program exits cleanly, so its own checkpoints pass.
  • Numerically \(b_n/2^n\) is \(0.222222222\) to nine decimals by \(n=100\), matching the \(2/9\) of Section 9.

The finite checks test the implementation; the derivation of Sections 5 to 8, with its boundary cases, establishes the recurrences themselves.

Footnotes and References

The references below supply the background facts used above; the reduction to subset sums, the seven limit forms and both recurrences are derived explicitly in this article.

  1. Project Euler 382 — Generating Polygons. The official statement says when a set of positive numbers generates a polygon — a polygon has at least three sides and does not self-intersect — illustrates it with {3, 4, 5} and {6, 9, 11, 24} against {1, 2, 3} and {2, 3, 4, 9}, defines s by s₁ = 1, s₂ = 2, s₃ = 3, sₙ = sₙ₋₁ + sₙ₋₃, and asks for the last nine digits of f(10¹⁸), where f(n) counts the subsets of Uₙ = {s₁, …, sₙ} that generate a polygon, with f(5) = 7, f(10) = 501 and f(25) = 18635853.
  2. Wikipedia — Triangle inequality. The article states that three lengths bound a triangle of positive area exactly when each is less than the sum of the other two, that it degenerates to zero area at equality, and that the inequality extends to polygonal paths, which are never shorter than the segment joining their endpoints. Section 1 uses the extended form for necessity.
  3. I. Pinelis, Cyclic polygons with given edge lengths: existence and uniqueness, Journal of Geometry 82 (2005), no. 1–2, 156–171. For a prescribed cyclic order of positive edge lengths, the paper proves that a convex cyclic polygon with those lengths exists and is unique up to isometry when each length is less than the sum of the others — a condition it notes is necessary for closure — and that uniqueness can fail without convexity. Section 1 proves only existence.
  4. OEIS A000930 — Narayana's cows sequence. The entry defines a(0) = a(1) = a(2) = 1 and a(n) = a(n−1) + a(n−3), lists 1, 1, 1, 2, 3, 4, 6, 9, 13, 19, 28, 41, 60, 88, 129, and records that a(n+1)/a(n) tends to 1.465571…, the real root of x³ − x² − 1. The sequence here is that one shifted, sₙ = a(n+1).
  5. Wikipedia — Supergolden ratio. The supergolden ratio is the unique real solution of x³ = x² + 1, about 1.4655712; the page cites Narayana's cows sequence A000930 alongside it and gives x³ + x − 1 as the minimal polynomial of its reciprocal, one of the two factors found in Section 9.
  6. Wikipedia — Plastic ratio. The plastic ratio is the unique real solution of x³ = x + 1, about 1.32471795724, with minimal polynomial x³ − x − 1. This is the other factor in Section 9.
  7. cp-algorithms — Binary exponentiation. The page shows that repeated squaring forms a power with O(log n) multiplications, and its application section raises a transition matrix to the n-th power to obtain a term of a linear recurrence. Section 10 applies this to the 9 × 9 matrix.
  8. Wikipedia — Constant-recursive sequence. The article records that a sequence of order at most d can be written as uAⁿv for a d × d matrix A, that its closed form comes from the roots of the characteristic polynomial, that the powers of two and the all-ones sequence are constant-recursive, and that a term-wise sum of such sequences is again one. Section 8 carries out that closure with shift operators.

Problem 382 source code

C++

#include <array>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = __uint128_t;
using i128 = __int128_t;

constexpr u64 kDefaultN = 1'000'000'000'000'000'000ULL;
constexpr u64 kMod = 1'000'000'000ULL;
constexpr int kOrder = 8;
constexpr std::array<i64, kOrder> kRecurrence = {3, -2, 2, -5, 1, 1, 3, -2};
constexpr int kValidationN = 60;

struct Options {
    u64 n = kDefaultN;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

u64 mod_norm(const i128 value, const u64 mod) {
    i128 x = value % static_cast<i128>(mod);
    if (x < 0) {
        x += static_cast<i128>(mod);
    }
    return static_cast<u64>(x);
}

u64 add_mod(const u64 a, const u64 b, const u64 mod) {
    const u64 sum = a + b;
    if (sum >= mod || sum < a) {
        return sum % mod;
    }
    return sum;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) %
                            static_cast<u128>(mod));
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }

    u64 result = 1ULL % mod;
    base %= mod;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1U;
    }
    return result;
}

struct Key {
    int i = 0;
    u64 limit = 0ULL;

    bool operator==(const Key& other) const {
        return i == other.i && limit == other.limit;
    }
};

struct KeyHash {
    std::size_t operator()(const Key& key) const {
        const std::size_t h1 = std::hash<int>{}(key.i);
        const std::size_t h2 = std::hash<u64>{}(key.limit);
        return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6U) + (h1 >> 2U));
    }
};

class ExactCounter {
public:
    explicit ExactCounter(const int max_n) : s_(static_cast<std::size_t>(max_n + 1), 0ULL),
                                             prefix_(static_cast<std::size_t>(max_n + 1), 0U) {
        if (max_n < 3) {
            throw std::runtime_error("ExactCounter requires max_n >= 3.");
        }

        s_[1] = 1ULL;
        s_[2] = 2ULL;
        s_[3] = 3ULL;
        for (int n = 4; n <= max_n; ++n) {
            s_[n] = s_[n - 1] + s_[n - 3];
        }
        for (int i = 1; i <= max_n; ++i) {
            prefix_[i] = prefix_[i - 1] + static_cast<u128>(s_[i]);
        }
    }

    const std::vector<u64>& sequence() const {
        return s_;
    }

    u128 count_subsets_leq(const int i, const u64 limit) {
        if (i <= 0) {
            return 1U;
        }

        if (static_cast<u128>(limit) >= prefix_[i]) {
            return static_cast<u128>(1U) << i;
        }

        const Key key{i, limit};
        const auto it = memo_.find(key);
        if (it != memo_.end()) {
            return it->second;
        }

        u128 result = 0U;
        if (s_[i] > limit) {
            result = count_subsets_leq(i - 1, limit);
        } else {
            result = count_subsets_leq(i - 1, limit) + count_subsets_leq(i - 1, limit - s_[i]);
        }

        memo_.emplace(key, result);
        return result;
    }

private:
    std::vector<u64> s_;
    std::vector<u128> prefix_;
    std::unordered_map<Key, u128, KeyHash> memo_;
};

struct Mat9 {
    std::array<std::array<u64, 9>, 9> a{};
};

Mat9 identity9() {
    Mat9 id;
    for (int i = 0; i < 9; ++i) {
        id.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(i)] = 1ULL;
    }
    return id;
}

Mat9 mat_mul(const Mat9& lhs, const Mat9& rhs, const u64 mod) {
    Mat9 out;
    for (int i = 0; i < 9; ++i) {
        for (int k = 0; k < 9; ++k) {
            const u64 lv = lhs.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)];
            if (lv == 0ULL) {
                continue;
            }
            for (int j = 0; j < 9; ++j) {
                const u64 rv = rhs.a[static_cast<std::size_t>(k)][static_cast<std::size_t>(j)];
                if (rv == 0ULL) {
                    continue;
                }
                const u64 add = mul_mod(lv, rv, mod);
                u64& cell = out.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
                cell = add_mod(cell, add, mod);
            }
        }
    }
    return out;
}

std::array<u64, 9> mat_vec_mul(const Mat9& m, const std::array<u64, 9>& v, const u64 mod) {
    std::array<u64, 9> out{};
    for (int i = 0; i < 9; ++i) {
        u64 acc = 0ULL;
        for (int j = 0; j < 9; ++j) {
            const u64 prod = mul_mod(m.a[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)],
                                     v[static_cast<std::size_t>(j)],
                                     mod);
            acc = add_mod(acc, prod, mod);
        }
        out[static_cast<std::size_t>(i)] = acc;
    }
    return out;
}

std::array<u64, 9> apply_mat_pow(Mat9 base,
                                 u64 exp,
                                 std::array<u64, 9> vec,
                                 const u64 mod) {
    Mat9 acc = identity9();
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            acc = mat_mul(base, acc, mod);
        }
        base = mat_mul(base, base, mod);
        exp >>= 1U;
    }
    return mat_vec_mul(acc, vec, mod);
}

u64 prefix_sum_b_mod(const u64 n, const u64 mod) {
    constexpr std::array<u64, kOrder> kBase = {1ULL, 2ULL, 4ULL, 6ULL, 11ULL, 20ULL, 36ULL,
                                                67ULL};

    if (n == 0ULL) {
        return 0ULL;
    }

    if (n <= static_cast<u64>(kOrder)) {
        u64 sum = 0ULL;
        for (u64 i = 0ULL; i < n; ++i) {
            sum = add_mod(sum, kBase[static_cast<std::size_t>(i)] % mod, mod);
        }
        return sum;
    }

    u64 s8 = 0ULL;
    for (const u64 v : kBase) {
        s8 = add_mod(s8, v % mod, mod);
    }

    std::array<u64, 9> state{
        kBase[7] % mod, kBase[6] % mod, kBase[5] % mod, kBase[4] % mod, kBase[3] % mod,
        kBase[2] % mod, kBase[1] % mod, kBase[0] % mod, s8};

    Mat9 trans;
    for (int j = 0; j < kOrder; ++j) {
        trans.a[0][static_cast<std::size_t>(j)] = mod_norm(static_cast<i128>(kRecurrence[j]), mod);
    }
    for (int r = 1; r < kOrder; ++r) {
        trans.a[static_cast<std::size_t>(r)][static_cast<std::size_t>(r - 1)] = 1ULL;
    }
    for (int j = 0; j < kOrder; ++j) {
        trans.a[8][static_cast<std::size_t>(j)] = mod_norm(static_cast<i128>(kRecurrence[j]), mod);
    }
    trans.a[8][8] = 1ULL;

    const std::array<u64, 9> advanced = apply_mat_pow(trans, n - 8ULL, state, mod);
    return advanced[8];
}

u64 solve_mod(const u64 n, const u64 mod) {
    const u64 sum_b = prefix_sum_b_mod(n, mod);
    const u64 pow2 = pow_mod(2ULL, n, mod);
    i128 ans = static_cast<i128>(pow2) - 1 - static_cast<i128>(sum_b);
    return mod_norm(ans, mod);
}

u64 to_u64(const u128 x) {
    return static_cast<u64>(x);
}

void run_checkpoints() {
    ExactCounter exact(kValidationN);
    const auto& s = exact.sequence();

    std::vector<u128> b_exact(static_cast<std::size_t>(kValidationN + 1), 0U);
    for (int n = 1; n <= kValidationN; ++n) {
        b_exact[static_cast<std::size_t>(n)] = exact.count_subsets_leq(n - 1, s[n]);
    }

    for (int n = 9; n <= kValidationN; ++n) {
        i128 rhs = 0;
        for (int j = 0; j < kOrder; ++j) {
            rhs += static_cast<i128>(kRecurrence[static_cast<std::size_t>(j)]) *
                   static_cast<i128>(b_exact[static_cast<std::size_t>(n - 1 - j)]);
        }
        const u128 lhs = b_exact[static_cast<std::size_t>(n)];
        if (rhs < 0 || static_cast<u128>(rhs) != lhs) {
            throw std::runtime_error("Derived linear recurrence failed validation.");
        }
    }

    const auto f_exact = [&](const int n) -> u128 {
        u128 sum_b = 0U;
        for (int i = 1; i <= n; ++i) {
            sum_b += b_exact[static_cast<std::size_t>(i)];
        }
        return (static_cast<u128>(1U) << n) - 1U - sum_b;
    };

    const u64 f5 = to_u64(f_exact(5));
    const u64 f10 = to_u64(f_exact(10));
    const u64 f25 = to_u64(f_exact(25));
    if (f5 != 7ULL || f10 != 501ULL || f25 != 18'635'853ULL) {
        throw std::runtime_error("Problem statement checkpoints failed.");
    }

    for (int n = 1; n <= kValidationN; ++n) {
        const u64 fast = solve_mod(static_cast<u64>(n), kMod);
        const u64 slow = static_cast<u64>(f_exact(n) % static_cast<u128>(kMod));
        if (fast != slow) {
            throw std::runtime_error("Fast modular solver mismatch on validation range.");
        }
    }
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    try {
        if (options.run_checkpoints) {
            run_checkpoints();
        }
    } catch (const std::exception& ex) {
        std::cerr << "Checkpoint failure: " << ex.what() << '\n';
        return 1;
    }

    const u64 answer = solve_mod(options.n, kMod);
    std::cout << "f(" << options.n << ") mod 10^9 = " << answer << '\n';
    std::cout << "Last 9 digits: " << std::setw(9) << std::setfill('0') << answer << '\n';
    return 0;
}

Python

def solve():
    kMod = 1000000000
    kOrder = 8
    kRecurrence = [3, -2, 2, -5, 1, 1, 3, -2]
    
    def zero_mat():
        return [[0]*9 for _ in range(9)]
        
    def mat_mul(a, b):
        c = zero_mat()
        for i in range(9):
            for k in range(9):
                if a[i][k] == 0: continue
                for j in range(9):
                    if b[k][j] == 0: continue
                    c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % kMod
        return c

    def apply_mat_pow(base, exp, vec):
        acc = zero_mat()
        for i in range(9): acc[i][i] = 1
        
        while exp > 0:
            if exp & 1:
                acc = mat_mul(base, acc)
            base = mat_mul(base, base)
            exp >>= 1
            
        out = [0]*9
        for i in range(9):
            s = 0
            for j in range(9):
                s = (s + acc[i][j] * vec[j]) % kMod
            out[i] = s
        return out

    def mod_norm(val):
        return val % kMod

    kBase = [1, 2, 4, 6, 11, 20, 36, 67]
    n = 1000000000000000000
    
    s8 = sum(kBase) % kMod
    state = [
        kBase[7]%kMod, kBase[6]%kMod, kBase[5]%kMod, kBase[4]%kMod,
        kBase[3]%kMod, kBase[2]%kMod, kBase[1]%kMod, kBase[0]%kMod, s8
    ]
    
    trans = zero_mat()
    for j in range(kOrder):
        trans[0][j] = mod_norm(kRecurrence[j])
    for r in range(1, kOrder):
        trans[r][r-1] = 1
    for j in range(kOrder):
        trans[8][j] = mod_norm(kRecurrence[j])
    trans[8][8] = 1
    
    advanced = apply_mat_pow(trans, n - 8, state)
    sum_b = advanced[8]
    
    pow2 = pow(2, n, kMod)
    ans = (pow2 - 1 - sum_b) % kMod
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

public class Euler382 {
    private static final long kMod = 1000000000L;
    private static final int kOrder = 8;
    private static final long[] kRecurrence = { 3, -2, 2, -5, 1, 1, 3, -2 };
    private static final long[] kBase = { 1, 2, 4, 6, 11, 20, 36, 67 };

    private static long[][] matMul(long[][] a, long[][] b) {
        long[][] c = new long[9][9];
        for (int i = 0; i < 9; ++i) {
            for (int k = 0; k < 9; ++k) {
                if (a[i][k] == 0)
                    continue;
                for (int j = 0; j < 9; ++j) {
                    if (b[k][j] == 0)
                        continue;
                    c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % kMod;
                }
            }
        }
        return c;
    }

    private static long modNorm(long val) {
        long x = val % kMod;
        if (x < 0)
            x += kMod;
        return x;
    }

    private static long[] applyMatPow(long[][] base, long exp, long[] vec) {
        long[][] acc = new long[9][9];
        for (int i = 0; i < 9; ++i)
            acc[i][i] = 1;

        while (exp > 0) {
            if ((exp & 1) != 0) {
                acc = matMul(base, acc);
            }
            base = matMul(base, base);
            exp >>= 1;
        }

        long[] out = new long[9];
        for (int i = 0; i < 9; ++i) {
            long s = 0;
            for (int j = 0; j < 9; ++j) {
                s = (s + acc[i][j] * vec[j]) % kMod;
            }
            out[i] = s;
        }
        return out;
    }

    private static long powMod(long base, long exp) {
        long res = 1;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * base) % kMod;
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return res;
    }

    public static String solve() {
        long n = 1000000000000000000L;
        long s8 = 0;
        for (long v : kBase) {
            s8 = (s8 + v) % kMod;
        }

        long[] state = {
                kBase[7] % kMod, kBase[6] % kMod, kBase[5] % kMod, kBase[4] % kMod,
                kBase[3] % kMod, kBase[2] % kMod, kBase[1] % kMod, kBase[0] % kMod, s8
        };

        long[][] trans = new long[9][9];
        for (int j = 0; j < kOrder; ++j) {
            trans[0][j] = modNorm(kRecurrence[j]);
        }
        for (int r = 1; r < kOrder; ++r) {
            trans[r][r - 1] = 1;
        }
        for (int j = 0; j < kOrder; ++j) {
            trans[8][j] = modNorm(kRecurrence[j]);
        }
        trans[8][8] = 1;

        long[] advanced = applyMatPow(trans, n - 8, state);
        long sumB = advanced[8];

        long pow2 = powMod(2, n);
        long ans = modNorm(pow2 - 1 - sumB);
        return String.valueOf(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}