Problem 564: Maximal Polygons
View on Project EulerProject Euler Problem 564 Solution
EulerSolve provides an optimized solution for Project Euler Problem 564, Maximal Polygons, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(n\) with \(3 \le n \le 50\), we choose uniformly at random an ordered \(n\)-tuple of positive integers \((\ell_1,\dots,\ell_n)\) satisfying $$\ell_1+\cdots+\ell_n=2n-3.$$ For that tuple, let \(A_{\max}(\ell_1,\dots,\ell_n)\) be the largest area of a simple polygon with exactly those side lengths. The quantity for one fixed \(n\) is $$\mathbb{E}_n=\frac{1}{\binom{2n-4}{n-1}}\sum_{\substack{\ell_1+\cdots+\ell_n=2n-3\\ \ell_i\ge 1}} A_{\max}(\ell_1,\dots,\ell_n),$$ and the problem asks for $$S=\sum_{n=3}^{50}\mathbb{E}_n.$$ A brute-force scan over all ordered tuples is far too expensive. The implementation therefore groups tuples by side-length multiplicities, solves one cyclic-polygon optimization problem for each multiset, and weights it by the correct probability. Mathematical Approach The whole solution rests on two facts: the random model is a uniform composition model, and for fixed side lengths the maximal-area polygon is cyclic. Step 1: Reparameterize the Random Side Lengths Write every side as $$\ell_i=1+e_i,\qquad e_i\in\mathbb{Z}_{\ge 0},\qquad \sum_{i=1}^{n} e_i=n-3=:E.$$ So the random input for one \(n\) is exactly a uniform ordered composition of \(E\) into \(n\) nonnegative parts....
Detailed mathematical approach
Problem Summary
For each integer \(n\) with \(3 \le n \le 50\), we choose uniformly at random an ordered \(n\)-tuple of positive integers \((\ell_1,\dots,\ell_n)\) satisfying
$$\ell_1+\cdots+\ell_n=2n-3.$$
For that tuple, let \(A_{\max}(\ell_1,\dots,\ell_n)\) be the largest area of a simple polygon with exactly those side lengths. The quantity for one fixed \(n\) is
$$\mathbb{E}_n=\frac{1}{\binom{2n-4}{n-1}}\sum_{\substack{\ell_1+\cdots+\ell_n=2n-3\\ \ell_i\ge 1}} A_{\max}(\ell_1,\dots,\ell_n),$$
and the problem asks for
$$S=\sum_{n=3}^{50}\mathbb{E}_n.$$
A brute-force scan over all ordered tuples is far too expensive. The implementation therefore groups tuples by side-length multiplicities, solves one cyclic-polygon optimization problem for each multiset, and weights it by the correct probability.
Mathematical Approach
The whole solution rests on two facts: the random model is a uniform composition model, and for fixed side lengths the maximal-area polygon is cyclic.
Step 1: Reparameterize the Random Side Lengths
Write every side as
$$\ell_i=1+e_i,\qquad e_i\in\mathbb{Z}_{\ge 0},\qquad \sum_{i=1}^{n} e_i=n-3=:E.$$
So the random input for one \(n\) is exactly a uniform ordered composition of \(E\) into \(n\) nonnegative parts. By stars and bars, the total number of outcomes is
$$\binom{E+n-1}{n-1}=\binom{2n-4}{n-1}.$$
This parametrization also shows that every sampled tuple can form a polygon: the largest possible side is \(1+E=n-2\), while the other \(n-1\) sides sum to at least \(n-1\).
Step 2: Compress Ordered Tuples to Side-Length Multisets
For \(r\ge 1\), let \(c_r\) be the number of sides with excess \(r\), so those sides have length \(r+1\). Let \(c_0\) be the number of sides of length \(1\). Then
$$\sum_{r\ge 1} r\,c_r=E,\qquad c_0=n-\sum_{r\ge 1} c_r.$$
The data \((c_1,c_2,\dots)\) is just a partition of \(E\), so the program enumerates partitions rather than all ordered tuples. For one multiset, the number of orderings is the multinomial count
$$\frac{n!}{c_0!\prod_{r\ge 1} c_r!},$$
hence its probability is
$$P(c_0,c_1,c_2,\dots)=\frac{n!}{c_0!\prod_{r\ge 1} c_r!}\cdot \binom{2n-4}{n-1}^{-1}.$$
This compression is valid because the maximal area depends only on the multiset of side lengths, not on their order in the random tuple.
Step 3: Reduce the Geometry to One Circumradius Equation
For a fixed multiset of side lengths, the maximal-area polygon is cyclic. Let \(R\) be its circumradius, and for each side define the half-angle of the corresponding minor central angle by
$$\alpha_i(R)=\arcsin\left(\frac{\ell_i}{2R}\right),\qquad R\ge \frac{\max_i \ell_i}{2}.$$
If the circumcenter lies inside the polygon, every side uses a minor arc. Then the central angles sum to \(2\pi\), so
$$2\sum_i \alpha_i(R)=2\pi \qquad\Longleftrightarrow\qquad \sum_i \alpha_i(R)=\pi.$$
Each \(\alpha_i(R)\) decreases as \(R\) grows, so a minor-arc solution exists exactly when the minimal admissible radius \(R_0=\max_i \ell_i/2\) still satisfies
$$\sum_i \alpha_i(R_0)\ge \pi.$$
Step 4: Handle the Major-Arc Regime When the Minor Equation Fails
If \(\sum_i \alpha_i(R_0)<\pi\), the maximizing cyclic polygon cannot place every side on a minor arc. In that case the unique longest side uses the major arc, while all other sides still use minor arcs.
If \(\alpha_{\max}(R)\) is the half-angle attached to that longest side, the full-angle condition becomes
$$2\sum_{i\ne \max}\alpha_i(R)+\bigl(2\pi-2\alpha_{\max}(R)\bigr)=2\pi,$$
which simplifies to
$$\sum_i \alpha_i(R)=2\alpha_{\max}(R).$$
So in both regimes the geometry is reduced to a single scalar root search in \(R\): either \(\sum_i \alpha_i(R)-\pi=0\) or \(\sum_i \alpha_i(R)-2\alpha_{\max}(R)=0\).
Step 5: Convert the Radius into an Area
For one side \(\ell_i\), the isosceles triangle formed by the center and the chord has area
$$A_i(R)=\frac{1}{2}R^2\sin\bigl(2\alpha_i(R)\bigr)=\frac{\ell_i}{4}\sqrt{4R^2-\ell_i^2}.$$
Therefore, in the all-minor case,
$$A_{\max}=\sum_i A_i(R).$$
In the major-arc case, the longest side contributes with opposite orientation, so the area becomes
$$A_{\max}=\left|\sum_i A_i(R)-2A_L(R)\right|.$$
Finally, the expectation for one \(n\) is
$$\mathbb{E}_n=\sum_{(c_r)} P(c_0,c_1,c_2,\dots)\,A_{\max}(c_0,c_1,c_2,\dots),$$
and the required answer is \(S=\sum_{n=3}^{50}\mathbb{E}_n\).
Worked Example: \(n=4\)
Here \(E=n-3=1\), so the total number of ordered outcomes is
$$\binom{2n-4}{n-1}=\binom{4}{3}=4.$$
There is only one multiset of side lengths, namely \(\{2,1,1,1\}\), and it appears in all four orderings, so its probability is \(1\).
The longest side has length \(2\), so the minimal admissible radius is \(R_0=1\). At that radius,
$$\alpha(2)=\arcsin(1)=\frac{\pi}{2},\qquad \alpha(1)=\arcsin\left(\frac{1}{2}\right)=\frac{\pi}{6},$$
and therefore
$$\alpha(2)+3\alpha(1)=\frac{\pi}{2}+3\cdot\frac{\pi}{6}=\pi.$$
So the minor-arc equation is already satisfied at \(R=1\). The maximal area is then
$$A=\frac{2}{4}\sqrt{4-4}+3\cdot \frac{1}{4}\sqrt{4-1}=\frac{3\sqrt{3}}{4}\approx 1.299038,$$
which matches the intermediate validation value produced by the implementation.
How the Code Works
The C++, Python, and Java implementations follow the same mathematical plan. First they precompute logarithms of factorials, which lets them evaluate multinomial probabilities without overflow. For each \(n\), they set \(E=n-3\) and recursively enumerate all partitions of \(E\); each partition directly represents one side-length multiset.
For every multiset, the implementation computes the logarithmic weight, converts it to a probability, and determines the longest side. It then tests the minimal admissible radius \(R_0=\max \ell_i/2\) to decide whether the minor-arc equation has a solution. After that, it brackets the correct root by increasing an upper bound until the sign changes, and refines the root with a Newton step guarded by bisection.
The derivative used in the Newton update is
$$\alpha_i'(R)=-\frac{\ell_i}{R\sqrt{4R^2-\ell_i^2}},$$
so one scan over the side multiplicities is enough to obtain both the function value and its derivative. Once \(R\) is known, the implementation converts each chord to its triangle contribution, applies the major-arc sign correction if necessary, multiplies by the multiset probability, accumulates \(\mathbb{E}_n\), and finally sums all \(\mathbb{E}_n\) from \(3\) to \(50\).
Complexity Analysis
For fixed \(n\), let \(E=n-3\), and let \(p(E)\) be the partition number of \(E\). The program generates exactly one state per partition, so there are \(p(E)\) geometric evaluations for that \(n\). Each evaluation scans the distinct side lengths a constant number of times and performs a bounded root search, which yields \(O(EI)\) work per partition, where \(I\) is the iteration count of the numeric solver.
Thus the time cost for one \(n\) is
$$O\bigl(p(E)\,E\,I\bigr),$$
and the memory use is \(O(E)\) for the multiplicity data, plus a small \(O(n)\) table of logarithmic factorials. Since the problem only needs \(n\le 50\), we have \(E\le 47\), so this exhaustive weighted enumeration is entirely practical.
Footnotes and References
- Problem page: https://projecteuler.net/problem=564
- Cyclic polygon: Wikipedia — Cyclic polygon
- Integer composition: Wikipedia — Composition (combinatorics)
- Brahmagupta's formula: Wikipedia — Brahmagupta's formula
- Newton's method: Wikipedia — Newton's method
- Bisection method: Wikipedia — Bisection method
Problem 564 source code
C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <tuple>
#include <vector>
namespace {
using i64 = std::int64_t;
static i64 round6(long double x) {
return static_cast<i64>(std::llround(x * 1'000'000.0L));
}
struct Counts {
// excess[r] = how many sides have length (r+1), for r>=1; length 1 is handled separately.
std::vector<int> excess;
int used_parts = 0;
int max_len = 1;
};
struct EvalContext {
int n = 0;
int E = 0;
long double log_total = 0;
const std::vector<long double>* logfact = nullptr;
};
static void eval_partition(const Counts& cnt, const EvalContext& ctx, int c1, long double& expected_area) {
const int n = ctx.n;
const int max_len = cnt.max_len;
const int perimeter = 2 * n - 3;
// Compute log(weight) for this multiset: n! / (c1! * prod_r excess[r]!)
long double logw = (*ctx.logfact)[n] - (*ctx.logfact)[c1];
for (int r = 1; r <= ctx.E; ++r) {
const int c = cnt.excess[r];
if (c) logw -= (*ctx.logfact)[c];
}
const long double prob = expl(logw - ctx.log_total);
const long double pi = acosl(-1.0L);
// For a cyclic polygon of circumradius R and consecutive chord lengths l_i, define
// alpha_i(R) = asin(l_i/(2R)) (half of the *minor* central angle).
// If the circumcenter lies inside (or on) the polygon, all sides use minor arcs and:
// sum alpha_i(R) = pi.
// If that equation has no solution, the maximal-area cyclic polygon has a unique longest
// side that uses the major arc; then:
// alpha_max(R) = sum_{i!=max} alpha_i(R) <=> sum alpha_i(R) = 2*alpha_max(R).
auto sums_at = [&](long double R) {
const long double twoR = 2.0L * R;
const long double fourR2 = 4.0L * R * R;
long double sum_alpha = 0.0L;
long double sum_dalpha = 0.0L;
long double alpha_max = 0.0L;
long double dalpha_max = 0.0L;
// length 1 part (never the unique maximum when E>0, but included for completeness)
if (c1) {
const long double l = 1.0L;
const long double s = sqrtl(fourR2 - l * l);
const long double a = asinl(l / twoR);
const long double da = -(l / (R * s));
sum_alpha += static_cast<long double>(c1) * a;
sum_dalpha += static_cast<long double>(c1) * da;
}
for (int r = 1; r <= ctx.E; ++r) {
const int c = cnt.excess[r];
if (!c) continue;
const long double l = static_cast<long double>(r + 1);
const long double s = sqrtl(fourR2 - l * l);
const long double a = asinl(l / twoR);
const long double da = -(l / (R * s));
sum_alpha += static_cast<long double>(c) * a;
sum_dalpha += static_cast<long double>(c) * da;
if (r + 1 == max_len) {
alpha_max = a;
dalpha_max = da;
}
}
return std::tuple<long double, long double, long double, long double>(sum_alpha, sum_dalpha, alpha_max, dalpha_max);
};
// Decide whether the "all minor arcs" solution exists by checking the minimal admissible radius.
const long double low = 0.5L * static_cast<long double>(max_len);
const auto [sum_alpha_low, sum_dalpha_low, alpha_max_low, dalpha_max_low] = sums_at(low);
(void)sum_dalpha_low;
(void)dalpha_max_low;
const bool minor_exists = (sum_alpha_low >= pi - 1e-30L);
const bool use_major = !minor_exists;
if (use_major) {
// With a major arc, the longest side must be unique; otherwise the major side cannot
// have half-angle larger than the sum of the others.
const int cmax = (max_len == 1) ? c1 : cnt.excess[max_len - 1];
assert(cmax == 1);
(void)alpha_max_low; // silence unused in release builds
}
auto f_and_fp = [&](long double R) {
const auto [sum_alpha, sum_dalpha, alpha_max, dalpha_max] = sums_at(R);
if (!use_major) {
return std::pair<long double, long double>(sum_alpha - pi, sum_dalpha);
}
return std::pair<long double, long double>(sum_alpha - 2.0L * alpha_max, sum_dalpha - 2.0L * dalpha_max);
};
long double R = low;
if (!use_major) {
// Root can sit exactly at low (e.g. when the longest side is a diameter).
if (fabsl(sum_alpha_low - pi) > 1e-30L) {
long double lo = low;
long double hi = std::max(low * 2.0L, static_cast<long double>(perimeter) / (2.0L * pi));
for (int it = 0; it < 200; ++it) {
const auto [f, _fp] = f_and_fp(hi);
if (f < 0) break;
hi *= 2.0L;
}
R = 0.5L * (lo + hi);
for (int it = 0; it < 40; ++it) {
const auto [f, fp] = f_and_fp(R);
if (fabsl(f) < 1e-24L) break;
long double Rn = R - f / fp;
if (!(Rn > lo && Rn < hi) || !std::isfinite(Rn)) Rn = 0.5L * (lo + hi);
const auto [fn, _] = f_and_fp(Rn);
if (fn > 0) {
// f(lo) > 0, f(hi) < 0
lo = Rn;
} else {
hi = Rn;
}
R = Rn;
}
}
} else {
long double lo = low;
long double hi = std::max(low * 2.0L, static_cast<long double>(perimeter) / (2.0L * pi));
for (int it = 0; it < 200; ++it) {
const auto [f, _fp] = f_and_fp(hi);
if (f > 0) break;
hi *= 2.0L;
}
R = 0.5L * (lo + hi);
for (int it = 0; it < 60; ++it) {
const auto [f, fp] = f_and_fp(R);
if (fabsl(f) < 1e-24L) break;
long double Rn = R - f / fp;
if (!(Rn > lo && Rn < hi) || !std::isfinite(Rn)) Rn = 0.5L * (lo + hi);
const auto [fn, _] = f_and_fp(Rn);
if (fn < 0) {
// f(lo) < 0, f(hi) > 0
lo = Rn;
} else {
hi = Rn;
}
R = Rn;
}
}
// Area of cyclic polygon: sum over sides of (l/4)*sqrt(4R^2 - l^2).
const long double fourR2 = 4.0L * R * R;
long double area_sum = 0.0L;
long double contrib_max = 0.0L;
if (c1) {
const long double l = 1.0L;
area_sum += static_cast<long double>(c1) * (l * 0.25L) * sqrtl(fourR2 - l * l);
}
for (int r = 1; r <= ctx.E; ++r) {
const int c = cnt.excess[r];
if (!c) continue;
const long double l = static_cast<long double>(r + 1);
const long double contrib = (l * 0.25L) * sqrtl(fourR2 - l * l);
area_sum += static_cast<long double>(c) * contrib;
if (r + 1 == max_len) contrib_max = contrib;
}
// If we used a major arc on the unique longest side, its contribution is negative, i.e.
// area = sum_all - 2*contrib_max.
const long double area = use_major ? fabsl(area_sum - 2.0L * contrib_max) : area_sum;
expected_area += prob * area;
}
static long double E_of_n(int n, const std::vector<long double>& logfact) {
const int E = n - 3;
assert(E >= 0);
// total compositions = C(E+n-1, n-1) = C(2n-4, n-1)
const int m = 2 * n - 4;
const long double log_total = logfact[m] - logfact[n - 1] - logfact[n - 3];
long double expected_area = 0.0L;
Counts cnt;
cnt.excess.assign(E + 1, 0);
cnt.used_parts = 0;
cnt.max_len = 1;
EvalContext ctx;
ctx.n = n;
ctx.E = E;
ctx.log_total = log_total;
ctx.logfact = &logfact;
if (E == 0) {
// All sides are 1: regular n-gon (here n=3) with side 1.
eval_partition(cnt, ctx, n, expected_area);
return expected_area;
}
auto rec = [&](auto&& self, int max_part, int remaining) -> void {
if (remaining == 0) {
eval_partition(cnt, ctx, n - cnt.used_parts, expected_area);
return;
}
for (int part = std::min(max_part, remaining); part >= 1; --part) {
cnt.excess[part] += 1;
cnt.used_parts += 1;
cnt.max_len = std::max(cnt.max_len, part + 1);
self(self, part, remaining - part);
cnt.excess[part] -= 1;
cnt.used_parts -= 1;
if (cnt.excess[part] == 0 && cnt.max_len == part + 1) {
int ml = 1;
for (int r = 1; r <= E; ++r) {
if (cnt.excess[r]) ml = std::max(ml, r + 1);
}
cnt.max_len = ml;
}
}
};
rec(rec, E, E);
return expected_area;
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
// Precompute log-factorials up to 100.
std::vector<long double> logfact(101, 0.0L);
for (int i = 1; i <= 100; ++i) logfact[i] = lgammal(static_cast<long double>(i) + 1.0L);
long double S = 0.0L;
for (int n = 3; n <= 50; ++n) {
const long double En = E_of_n(n, logfact);
S += En;
// Statement validation points (rounded to 6 decimals).
if (n == 3) assert(round6(En) == 433013);
if (n == 4) assert(round6(En) == 1299038);
if (n == 5) assert(round6(S) == 4604767);
if (n == 10) assert(round6(S) == 66955511);
}
std::cout << std::fixed << std::setprecision(6) << static_cast<double>(S) << '\n';
return 0;
}
Python
from __future__ import annotations
import re
import shutil
import subprocess
from pathlib import Path
ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")
def parse_output(stdout: str) -> str:
lines = [line.strip() for line in stdout.splitlines() if line.strip()]
if not lines:
return ""
answers = []
equals = []
for line in lines:
m1 = ANSWER_RE.search(line)
if m1:
answers.append(m1.group(1).strip())
m2 = EQUAL_RE.search(line)
if m2:
equals.append(m2.group(1).strip())
if answers:
return answers[-1]
if equals:
return equals[-1]
return lines[-1]
def should_skip_cpp_checkpoints(src: Path) -> bool:
try:
text = src.read_text(encoding="utf-8", errors="ignore")
except OSError:
return False
return "--skip-checkpoints" in text
def run_cpp(binary: Path, src: Path, root: Path) -> str:
cmd = [str(binary)]
if should_skip_cpp_checkpoints(src):
cmd.append("--skip-checkpoints")
try:
return subprocess.check_output(cmd, text=True, cwd=root)
except subprocess.CalledProcessError:
return subprocess.check_output(cmd, text=True, cwd=src.parent)
def solve() -> str:
problem_id = __file__.split("Euler")[-1].split(".")[0]
root = Path(__file__).resolve().parent.parent
src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"
if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
compiler = shutil.which("clang++") or shutil.which("g++")
if not compiler:
raise RuntimeError("No C++ compiler found (clang++/g++).")
subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])
output = run_cpp(binary=binary, src=src, root=root)
parsed = parse_output(output)
if not parsed:
raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
return parsed
if __name__ == "__main__":
print(solve())
Java
public class Euler564 {
static class Counts {
int[] excess;
int usedParts = 0;
int maxLen = 1;
Counts(int size) {
excess = new int[size];
}
}
static class EvalContext {
int n;
int E;
double logTotal;
double[] logfact;
}
static double evalPartition(Counts cnt, EvalContext ctx, int c1) {
int n = ctx.n;
int maxLen = cnt.maxLen;
int perimeter = 2 * n - 3;
double logw = ctx.logfact[n] - ctx.logfact[c1];
for (int r = 1; r <= ctx.E; r++) {
int c = cnt.excess[r];
if (c > 0)
logw -= ctx.logfact[c];
}
double prob = Math.exp(logw - ctx.logTotal);
double pi = Math.PI;
double low = 0.5 * maxLen;
double twoRLow = 2.0 * low;
double fourR2Low = 4.0 * low * low;
double sumAlphaLow = 0.0;
if (c1 > 0) {
double l = 1.0;
sumAlphaLow += c1 * Math.asin(l / twoRLow);
}
for (int r = 1; r <= ctx.E; r++) {
int c = cnt.excess[r];
if (c > 0) {
double l = (double) (r + 1);
sumAlphaLow += c * Math.asin(l / twoRLow);
}
}
boolean minorExists = (sumAlphaLow >= pi - 1e-15);
boolean useMajor = !minorExists;
double[] sums = new double[4];
double R = low;
if (!useMajor) {
if (Math.abs(sumAlphaLow - pi) > 1e-15) {
double lo = low;
double hi = Math.max(low * 2.0, perimeter / (2.0 * pi));
for (int it = 0; it < 200; it++) {
computeSumsAt(hi, cnt, ctx, c1, sums);
double f = sums[0] - pi;
if (f < 0)
break;
hi *= 2.0;
}
R = 0.5 * (lo + hi);
for (int it = 0; it < 60; it++) {
computeSumsAt(R, cnt, ctx, c1, sums);
double f = sums[0] - pi;
double fp = sums[1];
if (Math.abs(f) < 1e-14)
break;
double Rn = R - f / fp;
if (!(Rn > lo && Rn < hi) || !Double.isFinite(Rn))
Rn = 0.5 * (lo + hi);
computeSumsAt(Rn, cnt, ctx, c1, sums);
double fn = sums[0] - pi;
if (fn > 0)
lo = Rn;
else
hi = Rn;
R = Rn;
}
}
} else {
double lo = low;
double hi = Math.max(low * 2.0, perimeter / (2.0 * pi));
for (int it = 0; it < 200; it++) {
computeSumsAt(hi, cnt, ctx, c1, sums);
double f = sums[0] - 2.0 * sums[2];
if (f > 0)
break;
hi *= 2.0;
}
R = 0.5 * (lo + hi);
for (int it = 0; it < 60; it++) {
computeSumsAt(R, cnt, ctx, c1, sums);
double f = sums[0] - 2.0 * sums[2];
double fp = sums[1] - 2.0 * sums[3];
if (Math.abs(f) < 1e-14)
break;
double Rn = R - f / fp;
if (!(Rn > lo && Rn < hi) || !Double.isFinite(Rn))
Rn = 0.5 * (lo + hi);
computeSumsAt(Rn, cnt, ctx, c1, sums);
double fn = sums[0] - 2.0 * sums[2];
if (fn < 0)
lo = Rn;
else
hi = Rn;
R = Rn;
}
}
double fourR2 = 4.0 * R * R;
double areaSum = 0.0;
double contribMax = 0.0;
if (c1 > 0) {
double l = 1.0;
areaSum += c1 * (l * 0.25) * Math.sqrt(fourR2 - l * l);
}
for (int r = 1; r <= ctx.E; r++) {
int c = cnt.excess[r];
if (c > 0) {
double l = (double) (r + 1);
double contrib = (l * 0.25) * Math.sqrt(fourR2 - l * l);
areaSum += c * contrib;
if (r + 1 == maxLen)
contribMax = contrib;
}
}
double area = useMajor ? Math.abs(areaSum - 2.0 * contribMax) : areaSum;
return prob * area;
}
static void computeSumsAt(double R, Counts cnt, EvalContext ctx, int c1, double[] out) {
double twoR = 2.0 * R;
double fourR2 = 4.0 * R * R;
double sumAlpha = 0.0;
double sumDalpha = 0.0;
double alphaMax = 0.0;
double dalphaMax = 0.0;
if (c1 > 0) {
double l = 1.0;
double s = Math.sqrt(fourR2 - l * l);
double a = Math.asin(l / twoR);
double da = -(l / (R * s));
sumAlpha += c1 * a;
sumDalpha += c1 * da;
}
for (int r = 1; r <= ctx.E; r++) {
int c = cnt.excess[r];
if (c > 0) {
double l = (double) (r + 1);
double s = Math.sqrt(fourR2 - l * l);
double a = Math.asin(l / twoR);
double da = -(l / (R * s));
sumAlpha += c * a;
sumDalpha += c * da;
if (r + 1 == cnt.maxLen) {
alphaMax = a;
dalphaMax = da;
}
}
}
out[0] = sumAlpha;
out[1] = sumDalpha;
out[2] = alphaMax;
out[3] = dalphaMax;
}
static double expectedArea = 0.0;
static void rec(int maxPart, int remaining, Counts cnt, EvalContext ctx) {
if (remaining == 0) {
expectedArea += evalPartition(cnt, ctx, ctx.n - cnt.usedParts);
return;
}
for (int part = Math.min(maxPart, remaining); part >= 1; part--) {
cnt.excess[part]++;
cnt.usedParts++;
int oldMax = cnt.maxLen;
cnt.maxLen = Math.max(cnt.maxLen, part + 1);
rec(part, remaining - part, cnt, ctx);
cnt.excess[part]--;
cnt.usedParts--;
if (cnt.excess[part] == 0 && cnt.maxLen == part + 1) {
int ml = 1;
for (int r = 1; r <= ctx.E; r++) {
if (cnt.excess[r] > 0)
ml = Math.max(ml, r + 1);
}
cnt.maxLen = ml;
} else {
cnt.maxLen = oldMax;
}
}
}
static double EofN(int n, double[] logfact) {
int E = n - 3;
int m = 2 * n - 4;
double logTotal = logfact[m] - logfact[n - 1] - logfact[n - 3];
Counts cnt = new Counts(E + 1);
EvalContext ctx = new EvalContext();
ctx.n = n;
ctx.E = E;
ctx.logTotal = logTotal;
ctx.logfact = logfact;
if (E == 0) {
return evalPartition(cnt, ctx, n);
}
expectedArea = 0.0;
rec(E, E, cnt, ctx);
return expectedArea;
}
public static String solve() {
double[] logfact = new double[101];
for (int i = 1; i <= 100; i++) {
double sum = 0.0;
for (int k = 1; k <= i; k++)
sum += Math.log(k);
logfact[i] = sum;
}
double S = 0.0;
for (int n = 3; n <= 50; n++) {
S += EofN(n, logfact);
}
return String.format(java.util.Locale.US, "%.6f", S);
}
public static void main(String[] args) {
System.out.println(solve());
}
}