Problem 841: Regular Star Polygons

View on Project Euler

Project Euler Problem 841 Solution

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

Problem Summary For each relevant pair \((p,q)\), the task is to evaluate the oriented area of a regular star polygon with \(p\)-fold symmetry and step parameter \(q\), then add those area values over a Fibonacci-driven family of pairs. The C++, Python, and Java implementations do not build the polygon vertex by vertex and they do not solve line intersections explicitly. Instead, they use a trigonometric formula for one repeated sector contribution and then scale it by symmetry. With Fibonacci numbers \(F_1=F_2=1\), the implementations ultimately sum the area quantity for the pairs \((F_m,F_{m-2})\) with \(m=4,5,\dots,35\). The mathematical problem is therefore to turn each star polygon area into a stable one-dimensional sum and then accumulate those sums accurately. Mathematical Approach Let \(\mathcal{A}(p,q)\) denote the oriented area quantity used by the implementations. The core idea is that regular symmetry reduces the geometry to a single angular increment, and that increment can be handled with tangent identities instead of coordinate-heavy geometry....

Detailed mathematical approach

Problem Summary

For each relevant pair \((p,q)\), the task is to evaluate the oriented area of a regular star polygon with \(p\)-fold symmetry and step parameter \(q\), then add those area values over a Fibonacci-driven family of pairs. The C++, Python, and Java implementations do not build the polygon vertex by vertex and they do not solve line intersections explicitly. Instead, they use a trigonometric formula for one repeated sector contribution and then scale it by symmetry.

With Fibonacci numbers \(F_1=F_2=1\), the implementations ultimately sum the area quantity for the pairs \((F_m,F_{m-2})\) with \(m=4,5,\dots,35\). The mathematical problem is therefore to turn each star polygon area into a stable one-dimensional sum and then accumulate those sums accurately.

Mathematical Approach

Let \(\mathcal{A}(p,q)\) denote the oriented area quantity used by the implementations. The core idea is that regular symmetry reduces the geometry to a single angular increment, and that increment can be handled with tangent identities instead of coordinate-heavy geometry.

Step 1: Use the Fundamental Angle

For a \(p\)-fold regular figure, the natural angular unit is

$$\delta=\frac{\pi}{p}.$$

The implementations work with the sequence of angles

$$\theta_k=k\delta \qquad (k\ge 1).$$

Rather than describing the star polygon through Cartesian coordinates, the area computation is expressed directly in terms of \(\sin(\theta_k)\), \(\cos(\theta_k)\), and \(\tan(\theta_k)\). This is enough because every contribution depends only on consecutive angles \((k-1)\delta\) and \(k\delta\).

Step 2: Convert Each Local Contribution into a Tangent Difference

The local building block appearing in the implementations is

$$\frac{\sin\delta}{\cos((k-1)\delta)\cos(k\delta)}.$$

Now apply the standard identity

$$\tan u-\tan v=\frac{\sin(u-v)}{\cos u\cos v}.$$

With \(u=k\delta\) and \(v=(k-1)\delta\), we get \(u-v=\delta\), so

$$\frac{\sin\delta}{\cos((k-1)\delta)\cos(k\delta)}=\tan(k\delta)-\tan((k-1)\delta).$$

This matters because the geometry becomes a sum of simple first differences of tangent values. No explicit polygon clipping or line-intersection algebra is required.

Step 3: Assemble the Signed Sector Sum

The implementations form the alternating sum

$$T(p,q)=\sum_{k=2}^{q}(-1)^k\frac{\sin\delta}{\cos((k-1)\delta)\cos(k\delta)}.$$

Using the identity above, this is equivalently

$$T(p,q)=\sum_{k=2}^{q}(-1)^k\left(\tan(k\delta)-\tan((k-1)\delta)\right).$$

The final area quantity is then

$$\mathcal{A}(p,q)=p(-1)^q\left(T(p,q)-\tan\delta\right).$$

The factor \(p\) reflects the \(p\)-fold rotational symmetry: once one representative angular contribution is known, the whole figure is obtained by repeating it \(p\) times. The alternating sign inside the sum and the outer factor \((-1)^q\) encode the orientation pattern used by the formula.

Step 4: Advance the Angles Recursively

Computing \(\sin(k\delta)\) and \(\cos(k\delta)\) from scratch for every \(k\) would be unnecessary. The implementations instead use the addition formulas

$$\sin((k+1)\delta)=\sin(k\delta)\cos\delta+\cos(k\delta)\sin\delta,$$

$$\cos((k+1)\delta)=\cos(k\delta)\cos\delta-\sin(k\delta)\sin\delta.$$

So after one initial evaluation of \(\sin\delta\) and \(\cos\delta\), the next angle is generated by recurrence. This keeps the inner loop lightweight and matches exactly what the code does.

Step 5: Sum over Fibonacci-Indexed Pairs

The generated Fibonacci list begins

$$1,1,2,3,5,8,\dots$$

and the implementations take pairs separated by two positions. Rewritten in standard Fibonacci notation, the total is

$$\sum_{m=4}^{35}\mathcal{A}(F_m,F_{m-2}).$$

This indexing is important: it is more precise than saying merely that the inputs are “Fibonacci-based”. Each area evaluation is tied to one specific pair \((F_m,F_{m-2})\).

Worked Example: \(\mathcal{A}(8,3)\)

For \(p=8\) and \(q=3\),

$$\delta=\frac{\pi}{8},\qquad \tan\delta=\sqrt{2}-1,\qquad \tan(2\delta)=1,\qquad \tan(3\delta)=\sqrt{2}+1.$$

The alternating inner sum has two terms:

$$\begin{aligned} T(8,3) &=\left(\tan\frac{\pi}{4}-\tan\frac{\pi}{8}\right)-\left(\tan\frac{3\pi}{8}-\tan\frac{\pi}{4}\right)\\ &=\left(1-(\sqrt{2}-1)\right)-\left((\sqrt{2}+1)-1\right)\\ &=2-2\sqrt{2}. \end{aligned}$$

Subtract the extra base term:

$$T(8,3)-\tan\delta=(2-2\sqrt{2})-(\sqrt{2}-1)=3-3\sqrt{2}.$$

Because \(q=3\) is odd, \((-1)^q=-1\). Therefore

$$\mathcal{A}(8,3)=8(-1)(3-3\sqrt{2})=24(\sqrt{2}-1)\approx 9.9411254970.$$

This is exactly the small numerical checkpoint reproduced by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they generate the needed Fibonacci values. For each pair \((p,q)\), they compute \(\delta\), \(\sin\delta\), \(\cos\delta\), and \(\tan\delta\). Then an inner loop runs from \(k=2\) to \(q\), advances \(\sin(k\delta)\) and \(\cos(k\delta)\) with the recurrence above, forms the quotient \(\sin\delta/(\cos((k-1)\delta)\cos(k\delta))\), applies the alternating sign, and adds the result to a compensated running total.

After the inner loop ends, the implementation subtracts \(\tan\delta\), multiplies by \(p(-1)^q\), and obtains one area value. The outer Fibonacci sum is also accumulated with compensated summation, because many contributions have opposite signs and similar magnitudes. Finally the result is formatted to 10 decimal places.

Complexity Analysis

One evaluation of \(\mathcal{A}(p,q)\) uses \(O(q)\) arithmetic steps and \(O(1)\) extra memory. Over the full problem range, the total cost is

$$O\left(\sum_{m=4}^{35}F_{m-2}\right)=O(F_{35}),$$

which is a fixed and very manageable amount of work. In exact terms, the inner loop executes

$$\sum_{j=2}^{33}F_j=F_{35}-2=9{,}227{,}463$$

iterations in total. Memory usage remains \(O(1)\) apart from the short Fibonacci table.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=841
  2. Regular star polygons: Wikipedia - Star polygon
  3. Tangent identities: Wikipedia - List of trigonometric identities
  4. Fibonacci numbers: Wikipedia - Fibonacci number
  5. Kahan summation algorithm: Wikipedia - Kahan summation algorithm

Problem 841 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(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;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

double area_A(const int p, const int q) {
    const double dp = static_cast<double>(p);
    const double delta = std::acos(-1.0) / dp;
    const double sd = std::sin(delta);
    const double cd = std::cos(delta);
    const double tan_delta = sd / cd;

    double sum_diff = 0.0;
    if (q >= 2) {
        int k = 2;
        double sin_k = sd;
        double cos_k = cd;
        double sign = 1.0;
        double corr = 0.0;  // Kahan compensation.
        while (k <= q) {
            const double sin_next = sin_k * cd + cos_k * sd;
            const double cos_next = cos_k * cd - sin_k * sd;
            const double diff = sd / (cos_k * cos_next);
            const double term = sign * diff;

            const double y = term - corr;
            const double t = sum_diff + y;
            corr = (t - sum_diff) - y;
            sum_diff = t;

            sin_k = sin_next;
            cos_k = cos_next;
            sign = -sign;
            ++k;
        }
    }

    const double inner = sum_diff - tan_delta;
    const double outside_sign = ((q % 2) == 0) ? 1.0 : -1.0;
    return dp * outside_sign * inner;
}

std::vector<std::uint64_t> fibonacci_sequence(const int count) {
    std::vector<std::uint64_t> fibs;
    fibs.reserve(static_cast<std::size_t>(count));
    fibs.push_back(1ULL);
    fibs.push_back(1ULL);
    for (int i = 2; i < count; ++i) {
        fibs.push_back(fibs[static_cast<std::size_t>(i - 1)] + fibs[static_cast<std::size_t>(i - 2)]);
    }
    return fibs;
}

double solve() {
    const std::vector<std::uint64_t> fibs = fibonacci_sequence(36);  // F1..F36
    double sum = 0.0;
    double corr = 0.0;
    for (int n = 3; n <= 34; ++n) {
        const int p = static_cast<int>(fibs[static_cast<std::size_t>(n)]);
        const int q = static_cast<int>(fibs[static_cast<std::size_t>(n - 2)]);
        const double term = area_A(p, q);
        const double y = term - corr;
        const double t = sum + y;
        corr = (t - sum) - y;
        sum = t;
    }
    return sum;
}

bool run_checkpoints() {
    const double a1 = area_A(8, 3);
    if (std::fabs(a1 - 9.9411254970) > 5e-11) {
        std::cerr << "Checkpoint failed: A(8,3)" << '\n';
        return false;
    }
    const double a2 = area_A(130021, 50008);
    if (std::fabs(a2 - 10.9210371479) > 5e-11) {
        std::cerr << "Checkpoint failed: A(130021,50008)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << std::fixed << std::setprecision(10) << solve() << '\n';
    return 0;
}

Python

import math

def area_A(p, q):
    dp = float(p)
    delta = math.pi / dp
    sd = math.sin(delta)
    cd = math.cos(delta)
    tan_delta = sd / cd
    
    sum_diff = 0.0
    if q >= 2:
        k = 2
        sin_k = sd
        cos_k = cd
        sign = 1.0
        corr = 0.0  # Kahan compensation
        while k <= q:
            sin_next = sin_k * cd + cos_k * sd
            cos_next = cos_k * cd - sin_k * sd
            diff = sd / (cos_k * cos_next)
            term = sign * diff
            
            y = term - corr
            t = sum_diff + y
            corr = (t - sum_diff) - y
            sum_diff = t
            
            sin_k = sin_next
            cos_k = cos_next
            sign = -sign
            k += 1
            
    inner = sum_diff - tan_delta
    outside_sign = 1.0 if (q % 2) == 0 else -1.0
    return dp * outside_sign * inner

def fibonacci_sequence(count):
    fibs = [1, 1]
    for i in range(2, count):
        fibs.append(fibs[i - 1] + fibs[i - 2])
    return fibs

def solve():
    fibs = fibonacci_sequence(36)
    ans = 0.0
    corr = 0.0
    for n in range(3, 35):
        p = fibs[n]
        q = fibs[n - 2]
        term = area_A(p, q)
        y = term - corr
        t = ans + y
        corr = (t - ans) - y
        ans = t
        
    ans = round(ans * 10**10) / 10**10
    return f"{ans:.10f}"

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

Java

public class Euler841 {
    static double areaA(int p, int q) {
        double dp = (double) p;
        double delta = Math.PI / dp;
        double sd = Math.sin(delta);
        double cd = Math.cos(delta);
        double tanDelta = sd / cd;

        double sumDiff = 0.0;
        if (q >= 2) {
            int k = 2;
            double sinK = sd;
            double cosK = cd;
            double sign = 1.0;
            double corr = 0.0;

            while (k <= q) {
                double sinNext = sinK * cd + cosK * sd;
                double cosNext = cosK * cd - sinK * sd;
                double diff = sd / (cosK * cosNext);
                double term = sign * diff;

                double y = term - corr;
                double t = sumDiff + y;
                corr = (t - sumDiff) - y;
                sumDiff = t;

                sinK = sinNext;
                cosK = cosNext;
                sign = -sign;
                k++;
            }
        }

        double inner = sumDiff - tanDelta;
        double outsideSign = ((q % 2) == 0) ? 1.0 : -1.0;
        return dp * outsideSign * inner;
    }

    static long[] fibonacciSequence(int count) {
        long[] fibs = new long[count];
        fibs[0] = 1;
        fibs[1] = 1;
        for (int i = 2; i < count; ++i) {
            fibs[i] = fibs[i - 1] + fibs[i - 2];
        }
        return fibs;
    }

    public static String solve() {
        long[] fibs = fibonacciSequence(36);
        double sum = 0.0;
        double corr = 0.0;

        for (int n = 3; n <= 34; ++n) {
            int p = (int) fibs[n];
            int q = (int) fibs[n - 2];
            double term = areaA(p, q);
            double y = term - corr;
            double t = sum + y;
            corr = (t - sum) - y;
            sum = t;
        }

        sum = Math.round(sum * 1e10) / 1e10;
        return String.format(java.util.Locale.US, "%.10f", sum);
    }

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