Problem 547: Distance of Random Points Within Hollow Square Laminae
View on Project EulerProject Euler Problem 547 Solution
EulerSolve provides an optimized solution for Project Euler Problem 547, Distance of Random Points Within Hollow Square Laminae, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Fix an outer square \(O=[0,n]\times[0,n]\). A valid hollow lamina is obtained by removing an interior axis-aligned rectangle $$H=[a,a+x]\times[b,b+y],\qquad 1\le x,y\le n-2,\qquad 1\le a\le n-x-1,\qquad 1\le b\le n-y-1.$$ The remaining region is \(R=O\setminus H\), whose area is \(A=n^2-xy\). If two points are chosen independently and uniformly from \(R\), their expected distance is $$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$ The quantity \(S(n)\) is the sum of \(E(R)\) over all valid choices of \(x,y,a,b\). The program computes this total and rounds the final value to four decimal places. Mathematical Approach The implementation converts the continuous geometry into a finite collection of precomputed interaction tables. The key object is a distance kernel for two unit cells, and everything else is built from that kernel by counting how many cell pairs realize each offset. Step 1: Build the Unit-Cell Distance Kernel Partition every rectangle into unit cells....
Detailed mathematical approach
Problem Summary
Fix an outer square \(O=[0,n]\times[0,n]\). A valid hollow lamina is obtained by removing an interior axis-aligned rectangle
$$H=[a,a+x]\times[b,b+y],\qquad 1\le x,y\le n-2,\qquad 1\le a\le n-x-1,\qquad 1\le b\le n-y-1.$$
The remaining region is \(R=O\setminus H\), whose area is \(A=n^2-xy\). If two points are chosen independently and uniformly from \(R\), their expected distance is
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
The quantity \(S(n)\) is the sum of \(E(R)\) over all valid choices of \(x,y,a,b\). The program computes this total and rounds the final value to four decimal places.
Mathematical Approach
The implementation converts the continuous geometry into a finite collection of precomputed interaction tables. The key object is a distance kernel for two unit cells, and everything else is built from that kernel by counting how many cell pairs realize each offset.
Step 1: Build the Unit-Cell Distance Kernel
Partition every rectangle into unit cells. For nonnegative integer offsets \((d_x,d_y)\), define \(K(d_x,d_y)\) as the total distance integral for two random points drawn from two unit cells whose lower-left corners differ by \((d_x,d_y)\):
$$K(d_x,d_y)=\int_0^1\int_0^1\int_0^1\int_0^1 \sqrt{(d_x+u_1-u_2)^2+(d_y+v_1-v_2)^2}\,du_1\,du_2\,dv_1\,dv_2.$$
After introducing the differences \(s=u_1-u_2\) and \(t=v_1-v_2\), the density becomes triangular, so the same quantity can be written as
$$K(d_x,d_y)=\int_{-1}^{1}\int_{-1}^{1}(1-|s|)(1-|t|)\sqrt{(d_x+s)^2+(d_y+t)^2}\,ds\,dt.$$
The implementations exploit symmetry and evaluate this integral numerically on \([0,1]^2\) with Gauss-Legendre quadrature:
$$\begin{aligned} K(d_x,d_y)=\int_0^1\int_0^1 &(1-u)(1-v)\Bigl( \sqrt{(d_x+u)^2+(d_y+v)^2} +\sqrt{(d_x+u)^2+(d_y-v)^2}\\ &+\sqrt{(d_x-u)^2+(d_y+v)^2} +\sqrt{(d_x-u)^2+(d_y-v)^2} \Bigr)\,du\,dv. \end{aligned}$$
A useful anchor is the special case \(K(0,0)\), the classical mean distance between two random points in the unit square:
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
Step 2: Aggregate a Whole Rectangle from Cell Offsets
Let \(F(w,h)\) denote the total double integral of distance over a full \(w\times h\) rectangle:
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
If two unit cells have relative offset \((d_x,d_y)\), then there are exactly \((w-|d_x|)(h-|d_y|)\) ordered pairs of cells with that offset. Therefore
$$F(w,h)=\sum_{d_x=-(w-1)}^{w-1}\sum_{d_y=-(h-1)}^{h-1}(w-|d_x|)(h-|d_y|)\,K(|d_x|,|d_y|).$$
This formula is the backbone of the precomputed rectangle table. It gives the full distance integral, not yet the expected distance. To get the mean distance inside a solid rectangle alone, one would divide by \((wh)^2\).
Step 3: Use Inclusion-Exclusion for the Hollow Region
For any two regions \(A\) and \(B\), write
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
Since the lamina is \(R=O\setminus H\), the indicator identity \(1_R=1_O-1_H\) implies
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
Distance is symmetric, so \(F_{OH}=F_{HO}\). Hence the usable formula is
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
Here \(F_{OO}=F(n,n)\) comes from the outer square table, \(F_{HH}=F(x,y)\) comes from the hole-size table, and only the cross term \(F_{OH}\) depends on the placement \((a,b)\). Once \(F_{RR}\) is known, the expected distance for that lamina is
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
Step 4: Compute the Outer-Hole Cross Term with Prefix Sums
For each unit cell \((i,j)\) inside the \(n\times n\) outer square, define its interaction with the whole outer square by
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
If the hole occupies the block of cells
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
then the cross term is simply the sum of \(W_n(i,j)\) over that block:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
A two-dimensional prefix sum over the \(W_n\) table reduces each such rectangular query to \(O(1)\) time, which is essential because the same outer square must be combined with many hole sizes and many placements.
Step 5: Sum Over All Valid Holes
Combining the previous steps yields
$$S(n)=\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}\sum_{a=1}^{n-x-1}\sum_{b=1}^{n-y-1} \frac{F(n,n)-2F_{OH}(n;x,y,a,b)+F(x,y)}{(n^2-xy)^2}.$$
This formula is exactly what the implementations evaluate. The only numerical approximation occurs in the kernel table \(K(d_x,d_y)\), and the rest is deterministic table lookup, prefix-sum extraction, and accumulation.
Worked Example: \(n=3\) and \(n=4\)
When \(n=3\), there is only one valid hole: \(x=y=1\) placed at \(a=b=1\). The lamina area is \(3^2-1=8\), so
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
This matches the checkpoint used by the implementations.
For \(n=4\), the number of admissible hollow laminae is
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
That count is another checkpoint: there are \(4\) placements for a \(1\times 1\) hole, \(2\) for \(1\times 2\), \(2\) for \(2\times 1\), and \(1\) for \(2\times 2\).
How the Code Works
The C++, Python, and Java implementations all follow the same pipeline. First they generate Gauss-Legendre nodes and weights on \([0,1]\) from Legendre polynomial roots and derivatives. With those nodes they fill the symmetric kernel table \(K(d_x,d_y)\) for all offsets \(0\le d_x,d_y\le n-1\).
Next they build a table of rectangle self-interactions \(F(w,h)\) for every \(1\le w,h\le n\). For the requested \(n\), they then compute the outer-to-cell interaction grid \(W_n(i,j)\) and turn it into a two-dimensional prefix sum so that each placement-dependent cross term \(F_{OH}\) is recovered by one rectangle query.
Finally they iterate over every valid hole size and every legal position, combine the three terms \(F_{OO}\), \(F_{OH}\), and \(F_{HH}\) by inclusion-exclusion, divide by the squared lamina area, and add the result to the running total. The C++ implementation also distributes independent hole-size pairs across worker threads, while preserving the same mathematics as the Python and Java implementations.
Complexity Analysis
Let \(q\) be the Gauss-Legendre order. Building the kernel table requires \(O(n^2q^2)\) time and \(O(n^2)\) memory. Building the rectangle self-interaction table costs \(O(n^4)\) time, and forming the outer interaction grid together with its prefix sum also costs \(O(n^4)\) time. The final sweep over all hole sizes and placements is
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4),$$
with \(O(1)\) work per placement after preprocessing. Therefore, for fixed quadrature order, the overall algorithm runs in \(O(n^4)\) time and uses \(O(n^2)\) memory.
Footnotes and References
- Problem page: https://projecteuler.net/problem=547
- Gauss-Legendre quadrature: Wikipedia — Gauss-Legendre quadrature
- Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
- Two-dimensional prefix sums: Wikipedia — Summed-area table
- The exact value \(K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}\) is the classical mean distance between two random points in a unit square.
Problem 547 source code
C++
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <utility>
#include <vector>
namespace {
constexpr int kDefaultN = 40;
constexpr int kDefaultQuadratureOrder = 64;
struct Options {
int n = kDefaultN;
bool allow_multithreading = true;
bool run_checkpoints = true;
unsigned requested_threads = 0;
int quadrature_order = kDefaultQuadratureOrder;
};
bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
std::int64_t parsed = 0;
for (char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const int digit = c - '0';
if (parsed > (std::numeric_limits<std::int64_t>::max() - digit) / 10) {
return false;
}
parsed = parsed * 10 + digit;
}
if (parsed > std::numeric_limits<int>::max()) {
return false;
}
value = static_cast<int>(parsed);
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
int parsed = 0;
if (!parse_int_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed < 0) {
return false;
}
value = static_cast<unsigned>(parsed);
return 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 == "--single-thread") {
options.allow_multithreading = false;
continue;
}
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
int n = 0;
if (parse_int_after_prefix(arg, "--n=", n)) {
options.n = n;
continue;
}
int q = 0;
if (parse_int_after_prefix(arg, "--quad=", q)) {
options.quadrature_order = q;
continue;
}
unsigned t = 0;
if (parse_unsigned_after_prefix(arg, "--threads=", t)) {
options.requested_threads = t;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.n < 3) {
std::cerr << "n must be at least 3.\n";
return false;
}
if (options.quadrature_order < 8 || (options.quadrature_order % 2) != 0) {
std::cerr << "--quad must be an even integer >= 8.\n";
return false;
}
return true;
}
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload) {
if (!allow_multithreading || workload < 2) {
return 1;
}
unsigned threads = requested_threads;
if (threads == 0) {
threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
}
threads = std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
return threads;
}
void legendre_polynomial_and_derivative(int n, long double x, long double& pn, long double& dpn) {
long double pnm1 = 1.0L;
long double pn_local = x;
if (n == 0) {
pn = 1.0L;
dpn = 0.0L;
return;
}
if (n == 1) {
pn = x;
dpn = 1.0L;
return;
}
for (int k = 2; k <= n; ++k) {
const long double pk =
((2.0L * k - 1.0L) * x * pn_local - (k - 1.0L) * pnm1) / static_cast<long double>(k);
pnm1 = pn_local;
pn_local = pk;
}
pn = pn_local;
dpn = static_cast<long double>(n) * (x * pn_local - pnm1) / (x * x - 1.0L);
}
void build_gauss_legendre_01(int order,
std::vector<long double>& nodes,
std::vector<long double>& weights) {
nodes.assign(order, 0.0L);
weights.assign(order, 0.0L);
const int half = (order + 1) / 2;
const long double pi = std::acos(-1.0L);
for (int i = 0; i < half; ++i) {
long double z =
std::cos(pi * (static_cast<long double>(i) + 0.75L) / (static_cast<long double>(order) + 0.5L));
for (int iter = 0; iter < 80; ++iter) {
long double pn = 0.0L;
long double dpn = 0.0L;
legendre_polynomial_and_derivative(order, z, pn, dpn);
const long double dz = pn / dpn;
z -= dz;
if (std::abs(dz) < 1e-25L) {
break;
}
}
long double pn = 0.0L;
long double dpn = 0.0L;
legendre_polynomial_and_derivative(order, z, pn, dpn);
const long double w = 2.0L / ((1.0L - z * z) * dpn * dpn);
const int j = order - 1 - i;
nodes[i] = (1.0L - z) * 0.5L;
nodes[j] = (1.0L + z) * 0.5L;
weights[i] = w * 0.5L;
weights[j] = w * 0.5L;
}
}
std::vector<std::vector<long double>> build_kernel_table(int max_n, int quadrature_order) {
const int max_diff = max_n - 1;
std::vector<std::vector<long double>> kernel(max_diff + 1,
std::vector<long double>(max_diff + 1, 0.0L));
std::vector<long double> nodes;
std::vector<long double> weights;
build_gauss_legendre_01(quadrature_order, nodes, weights);
for (int dx = 0; dx <= max_diff; ++dx) {
for (int dy = 0; dy <= max_diff; ++dy) {
long double sum = 0.0L;
for (int i = 0; i < quadrature_order; ++i) {
const long double u = nodes[i];
const long double wu = weights[i];
const long double ux = 1.0L - u;
for (int j = 0; j < quadrature_order; ++j) {
const long double v = nodes[j];
const long double wv = weights[j];
const long double vy = 1.0L - v;
const long double p1 = std::hypotl(static_cast<long double>(dx) + u,
static_cast<long double>(dy) + v);
const long double p2 = std::hypotl(static_cast<long double>(dx) + u,
static_cast<long double>(dy) - v);
const long double p3 = std::hypotl(static_cast<long double>(dx) - u,
static_cast<long double>(dy) + v);
const long double p4 = std::hypotl(static_cast<long double>(dx) - u,
static_cast<long double>(dy) - v);
sum += wu * wv * ux * vy * (p1 + p2 + p3 + p4);
}
}
kernel[dx][dy] = sum;
}
}
return kernel;
}
inline long double kernel_at(const std::vector<std::vector<long double>>& kernel, int dx, int dy) {
return kernel[std::abs(dx)][std::abs(dy)];
}
std::vector<std::vector<long double>> build_rect_self_table(
int max_n,
const std::vector<std::vector<long double>>& kernel) {
std::vector<std::vector<long double>> self(max_n + 1,
std::vector<long double>(max_n + 1, 0.0L));
for (int w = 1; w <= max_n; ++w) {
for (int h = 1; h <= max_n; ++h) {
long double sum = 0.0L;
for (int dx = -(w - 1); dx <= (w - 1); ++dx) {
const long double cx = static_cast<long double>(w - std::abs(dx));
for (int dy = -(h - 1); dy <= (h - 1); ++dy) {
const long double cy = static_cast<long double>(h - std::abs(dy));
sum += cx * cy * kernel_at(kernel, dx, dy);
}
}
self[w][h] = sum;
}
}
return self;
}
std::vector<std::vector<long double>> build_outer_prefix(
int n,
const std::vector<std::vector<long double>>& kernel) {
std::vector<std::vector<long double>> w(n, std::vector<long double>(n, 0.0L));
for (int hx = 0; hx < n; ++hx) {
for (int hy = 0; hy < n; ++hy) {
long double sum = 0.0L;
for (int ox = 0; ox < n; ++ox) {
for (int oy = 0; oy < n; ++oy) {
sum += kernel_at(kernel, ox - hx, oy - hy);
}
}
w[hx][hy] = sum;
}
}
std::vector<std::vector<long double>> prefix(n + 1, std::vector<long double>(n + 1, 0.0L));
for (int x = 0; x < n; ++x) {
long double row_acc = 0.0L;
for (int y = 0; y < n; ++y) {
row_acc += w[x][y];
prefix[x + 1][y + 1] = prefix[x][y + 1] + row_acc;
}
}
return prefix;
}
inline long double rect_sum(const std::vector<std::vector<long double>>& prefix,
int x0,
int y0,
int w,
int h) {
const int x1 = x0 + w;
const int y1 = y0 + h;
return prefix[x1][y1] - prefix[x0][y1] - prefix[x1][y0] + prefix[x0][y0];
}
long double compute_s_for_n(int n,
const std::vector<std::vector<long double>>& kernel,
const std::vector<std::vector<long double>>& rect_self,
bool allow_multithreading,
unsigned requested_threads) {
if (n < 3) {
return 0.0L;
}
const std::vector<std::vector<long double>> prefix = build_outer_prefix(n, kernel);
const long double f_outer_outer = rect_self[n][n];
std::vector<std::pair<int, int>> tasks;
tasks.reserve(static_cast<std::size_t>(n - 2) * static_cast<std::size_t>(n - 2));
for (int x = 1; x <= n - 2; ++x) {
for (int y = 1; y <= n - 2; ++y) {
tasks.emplace_back(x, y);
}
}
const unsigned thread_count =
choose_thread_count(allow_multithreading, requested_threads, tasks.size());
std::vector<long double> partial(thread_count, 0.0L);
std::atomic<std::size_t> next_task(0);
auto worker = [&](unsigned tid) {
long double local_sum = 0.0L;
while (true) {
const std::size_t idx = next_task.fetch_add(1, std::memory_order_relaxed);
if (idx >= tasks.size()) {
break;
}
const int x = tasks[idx].first;
const int y = tasks[idx].second;
const long double f_hole_hole = rect_self[x][y];
const long double area_ring =
static_cast<long double>(n) * static_cast<long double>(n) -
static_cast<long double>(x) * static_cast<long double>(y);
const long double inv_area_sq = 1.0L / (area_ring * area_ring);
for (int a = 1; a <= n - x - 1; ++a) {
for (int b = 1; b <= n - y - 1; ++b) {
const long double f_outer_hole = rect_sum(prefix, a, b, x, y);
const long double f_ring_ring =
f_outer_outer - 2.0L * f_outer_hole + f_hole_hole;
local_sum += f_ring_ring * inv_area_sq;
}
}
}
partial[tid] = local_sum;
};
std::vector<std::thread> threads;
threads.reserve(thread_count > 0 ? thread_count - 1 : 0);
for (unsigned t = 1; t < thread_count; ++t) {
threads.emplace_back(worker, t);
}
worker(0);
for (auto& th : threads) {
th.join();
}
long double total = 0.0L;
for (long double v : partial) {
total += v;
}
return total;
}
long long rounded4(long double x) { return std::llround(x * 10000.0L); }
bool run_checkpoints(const std::vector<std::vector<long double>>& kernel,
const std::vector<std::vector<long double>>& rect_self) {
// Checkpoint 1: unit-square mean distance.
const long double sqrt2 = std::sqrt(2.0L);
const long double expected_k00 =
(2.0L + sqrt2 + 5.0L * std::log(1.0L + sqrt2)) / 15.0L;
const long double got_k00 = kernel[0][0];
if (std::abs(got_k00 - expected_k00) > 1e-10L) {
std::cerr << "Checkpoint failed: kernel(0,0) mismatch. got=" << std::setprecision(16)
<< got_k00 << " expected=" << expected_k00 << '\n';
return false;
}
// Checkpoint 2: n=4 has 9 hollow laminae.
long long laminae_count = 0;
for (int x = 1; x <= 2; ++x) {
for (int y = 1; y <= 2; ++y) {
laminae_count += static_cast<long long>(4 - x - 1) * static_cast<long long>(4 - y - 1);
}
}
if (laminae_count != 9LL) {
std::cerr << "Checkpoint failed: laminae count for n=4 is " << laminae_count
<< ", expected 9\n";
return false;
}
// Checkpoint 3: statement samples.
const long double s3 = compute_s_for_n(3, kernel, rect_self, false, 1);
const long double s4 = compute_s_for_n(4, kernel, rect_self, false, 1);
if (rounded4(s3) != 16514LL) {
std::cerr << "Checkpoint failed: S(3) rounded to 4 decimals is " << std::fixed
<< std::setprecision(4) << s3 << ", expected 1.6514\n";
return false;
}
if (rounded4(s4) != 196564LL) {
std::cerr << "Checkpoint failed: S(4) rounded to 4 decimals is " << std::fixed
<< std::setprecision(4) << s4 << ", expected 19.6564\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
const int max_n = std::max(options.n, 4);
const std::vector<std::vector<long double>> kernel =
build_kernel_table(max_n, options.quadrature_order);
const std::vector<std::vector<long double>> rect_self = build_rect_self_table(max_n, kernel);
if (options.run_checkpoints && !run_checkpoints(kernel, rect_self)) {
return 1;
}
const long double answer = compute_s_for_n(
options.n, kernel, rect_self, options.allow_multithreading, options.requested_threads);
std::cout << std::fixed << std::setprecision(4) << answer << '\n';
return 0;
}
Python
import math
def legendre_polynomial_and_derivative(n, x):
pnm1 = 1.0
pn_local = x
if n == 0:
return 1.0, 0.0
if n == 1:
return x, 1.0
for k in range(2, n + 1):
pk = ((2.0 * k - 1.0) * x * pn_local - (k - 1.0) * pnm1) / k
pnm1 = pn_local
pn_local = pk
dpn = n * (x * pn_local - pnm1) / (x * x - 1.0)
return pn_local, dpn
def build_gauss_legendre_01(order):
nodes = [0.0] * order
weights = [0.0] * order
half = (order + 1) // 2
for i in range(half):
z = math.cos(math.pi * (i + 0.75) / (order + 0.5))
for _ in range(80):
pn, dpn = legendre_polynomial_and_derivative(order, z)
dz = pn / dpn
z -= dz
if abs(dz) < 1e-15:
break
pn, dpn = legendre_polynomial_and_derivative(order, z)
w = 2.0 / ((1.0 - z * z) * dpn * dpn)
j = order - 1 - i
nodes[i] = (1.0 - z) * 0.5
nodes[j] = (1.0 + z) * 0.5
weights[i] = w * 0.5
weights[j] = w * 0.5
return nodes, weights
def build_kernel_table(max_n, q_order):
max_diff = max_n - 1
kernel = [[0.0] * (max_diff + 1) for _ in range(max_diff + 1)]
nodes, weights = build_gauss_legendre_01(q_order)
for dx in range(max_diff + 1):
for dy in range(max_diff + 1):
s = 0.0
for i in range(q_order):
u = nodes[i]
wu = weights[i]
ux = 1.0 - u
for j in range(q_order):
v = nodes[j]
wv = weights[j]
vy = 1.0 - v
p1 = math.hypot(dx + u, dy + v)
p2 = math.hypot(dx + u, dy - v)
p3 = math.hypot(dx - u, dy + v)
p4 = math.hypot(dx - u, dy - v)
s += wu * wv * ux * vy * (p1 + p2 + p3 + p4)
kernel[dx][dy] = s
return kernel
def build_rect_self_table(max_n, kernel):
self_table = [[0.0] * (max_n + 1) for _ in range(max_n + 1)]
for w in range(1, max_n + 1):
for h in range(1, max_n + 1):
s = 0.0
for dx in range(-(w - 1), w):
cx = w - abs(dx)
for dy in range(-(h - 1), h):
cy = h - abs(dy)
s += cx * cy * kernel[abs(dx)][abs(dy)]
self_table[w][h] = s
return self_table
def build_outer_prefix(n, kernel):
w = [[0.0] * n for _ in range(n)]
for hx in range(n):
for hy in range(n):
s = 0.0
for ox in range(n):
for oy in range(n):
s += kernel[abs(ox - hx)][abs(oy - hy)]
w[hx][hy] = s
prefix = [[0.0] * (n + 1) for _ in range(n + 1)]
for x in range(n):
row_acc = 0.0
for y in range(n):
row_acc += w[x][y]
prefix[x + 1][y + 1] = prefix[x][y + 1] + row_acc
return prefix
def compute_s_for_n(n, kernel, rect_self):
if n < 3: return 0.0
prefix = build_outer_prefix(n, kernel)
f_outer_outer = rect_self[n][n]
total = 0.0
for x in range(1, n - 1):
for y in range(1, n - 1):
f_hole_hole = rect_self[x][y]
area_ring = n * n - x * y
inv_area_sq = 1.0 / (area_ring * area_ring)
local_sum = 0.0
for a in range(1, n - x):
for b in range(1, n - y):
f_outer_hole = prefix[a + x][b + y] - prefix[a][b + y] - prefix[a + x][b] + prefix[a][b]
f_ring_ring = f_outer_outer - 2.0 * f_outer_hole + f_hole_hole
local_sum += f_ring_ring * inv_area_sq
total += local_sum
return total
def solve():
n = 40
q_order = 64
kernel = build_kernel_table(n, q_order)
rect_self = build_rect_self_table(n, kernel)
ans = compute_s_for_n(n, kernel, rect_self)
return f"{ans:.4f}"
if __name__ == '__main__':
print(solve())
Java
public class Euler547 {
static void legendrePolyAndDeriv(int n, double x, double[] out) {
double pnm1 = 1.0;
double pnLocal = x;
if (n == 0) {
out[0] = 1.0;
out[1] = 0.0;
return;
}
if (n == 1) {
out[0] = x;
out[1] = 1.0;
return;
}
for (int k = 2; k <= n; k++) {
double pk = ((2.0 * k - 1.0) * x * pnLocal - (k - 1.0) * pnm1) / k;
pnm1 = pnLocal;
pnLocal = pk;
}
out[0] = pnLocal;
out[1] = n * (x * pnLocal - pnm1) / (x * x - 1.0);
}
static void buildGaussLegendre(int order, double[] nodes, double[] weights) {
int half = (order + 1) / 2;
double pi = Math.PI;
double[] out = new double[2];
for (int i = 0; i < half; i++) {
double z = Math.cos(pi * (i + 0.75) / (order + 0.5));
for (int iter = 0; iter < 80; iter++) {
legendrePolyAndDeriv(order, z, out);
double pn = out[0], dpn = out[1];
double dz = pn / dpn;
z -= dz;
if (Math.abs(dz) < 1e-15) {
break;
}
}
legendrePolyAndDeriv(order, z, out);
double pn = out[0], dpn = out[1];
double w = 2.0 / ((1.0 - z * z) * dpn * dpn);
int j = order - 1 - i;
nodes[i] = (1.0 - z) * 0.5;
nodes[j] = (1.0 + z) * 0.5;
weights[i] = w * 0.5;
weights[j] = w * 0.5;
}
}
static double[][] buildKernelTable(int maxN, int qOrder) {
int maxDiff = maxN - 1;
double[][] kernel = new double[maxDiff + 1][maxDiff + 1];
double[] nodes = new double[qOrder];
double[] weights = new double[qOrder];
buildGaussLegendre(qOrder, nodes, weights);
for (int dx = 0; dx <= maxDiff; dx++) {
for (int dy = 0; dy <= maxDiff; dy++) {
double s = 0.0;
for (int i = 0; i < qOrder; i++) {
double u = nodes[i];
double wu = weights[i];
double ux = 1.0 - u;
for (int j = 0; j < qOrder; j++) {
double v = nodes[j];
double wv = weights[j];
double vy = 1.0 - v;
double p1 = Math.hypot(dx + u, dy + v);
double p2 = Math.hypot(dx + u, dy - v);
double p3 = Math.hypot(dx - u, dy + v);
double p4 = Math.hypot(dx - u, dy - v);
s += wu * wv * ux * vy * (p1 + p2 + p3 + p4);
}
}
kernel[dx][dy] = s;
}
}
return kernel;
}
static double[][] buildRectSelfTable(int maxN, double[][] kernel) {
double[][] self = new double[maxN + 1][maxN + 1];
for (int w = 1; w <= maxN; w++) {
for (int h = 1; h <= maxN; h++) {
double s = 0.0;
for (int dx = -(w - 1); dx <= w - 1; dx++) {
double cx = w - Math.abs(dx);
for (int dy = -(h - 1); dy <= h - 1; dy++) {
double cy = h - Math.abs(dy);
s += cx * cy * kernel[Math.abs(dx)][Math.abs(dy)];
}
}
self[w][h] = s;
}
}
return self;
}
static double[][] buildOuterPrefix(int n, double[][] kernel) {
double[][] w = new double[n][n];
for (int hx = 0; hx < n; hx++) {
for (int hy = 0; hy < n; hy++) {
double s = 0.0;
for (int ox = 0; ox < n; ox++) {
for (int oy = 0; oy < n; oy++) {
s += kernel[Math.abs(ox - hx)][Math.abs(oy - hy)];
}
}
w[hx][hy] = s;
}
}
double[][] prefix = new double[n + 1][n + 1];
for (int x = 0; x < n; x++) {
double rowAcc = 0.0;
for (int y = 0; y < n; y++) {
rowAcc += w[x][y];
prefix[x + 1][y + 1] = prefix[x][y + 1] + rowAcc;
}
}
return prefix;
}
static double computeSForN(int n, double[][] kernel, double[][] rectSelf) {
if (n < 3)
return 0.0;
double[][] prefix = buildOuterPrefix(n, kernel);
double fOuterOuter = rectSelf[n][n];
double total = 0.0;
for (int x = 1; x <= n - 2; x++) {
for (int y = 1; y <= n - 2; y++) {
double fHoleHole = rectSelf[x][y];
double areaRing = (double) n * n - (double) x * y;
double invAreaSq = 1.0 / (areaRing * areaRing);
double localSum = 0.0;
for (int a = 1; a <= n - x - 1; a++) {
for (int b = 1; b <= n - y - 1; b++) {
double fOuterHole = prefix[a + x][b + y] - prefix[a][b + y] - prefix[a + x][b] + prefix[a][b];
double fRingRing = fOuterOuter - 2.0 * fOuterHole + fHoleHole;
localSum += fRingRing * invAreaSq;
}
}
total += localSum;
}
}
return total;
}
public static String solve() {
int n = 40;
int qOrder = 64;
double[][] kernel = buildKernelTable(n, qOrder);
double[][] rectSelf = buildRectSelfTable(n, kernel);
double ans = computeSForN(n, kernel, rectSelf);
return String.format(java.util.Locale.US, "%.4f", ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}