Problem 580: Squarefree Hilbert Numbers

Project Euler Problem 580 Solution

EulerSolve provides an optimized solution for Project Euler Problem 580, Squarefree Hilbert Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A Hilbert number is a positive integer of the form $4k+1$. The problem asks for the number of Hilbert numbers below a large exclusive limit $X$ that are squarefree inside the Hilbert semigroup: equivalently, numbers $n\equiv1\pmod4$ such that no Hilbert number $h>1$ satisfies $h^2\mid n$. The implementations work with the inclusive bound $N=X-1$, so the target quantity is $$S(N)=\#\{n\le N:n\equiv1\pmod4,\ \nexists h\equiv1\pmod4,\ h>1,\ h^2\mid n\}.$$ A direct scan with explicit Hilbert-square testing is far too slow near $10^{16}$, so the solution rewrites the condition as a Möbius-weighted sum over odd square divisors. Mathematical Approach The main job is to understand which ordinary prime exponents are allowed in a squarefree Hilbert number, and then encode that description with inclusion-exclusion. Step 1: Translate the Hilbert-square condition into ordinary prime exponents Write an odd number $n\equiv1\pmod4$ as $$n=\prod_{p\equiv1\pmod4} p^{a_p}\prod_{q\equiv3\pmod4} q^{b_q}.$$ If some prime $p\equiv1\pmod4$ has $a_p\ge2$, then $p$ itself is a Hilbert number and $p^2\mid n$, so $n$ is not squarefree in the Hilbert sense. For primes $q\equiv3\pmod4$, one copy of $q^2$ is still allowed because $q\not\equiv1\pmod4$. But two square units are too many....

Detailed mathematical approach

Problem Summary

A Hilbert number is a positive integer of the form $4k+1$. The problem asks for the number of Hilbert numbers below a large exclusive limit $X$ that are squarefree inside the Hilbert semigroup: equivalently, numbers $n\equiv1\pmod4$ such that no Hilbert number $h>1$ satisfies $h^2\mid n$.

The implementations work with the inclusive bound $N=X-1$, so the target quantity is

$$S(N)=\#\{n\le N:n\equiv1\pmod4,\ \nexists h\equiv1\pmod4,\ h>1,\ h^2\mid n\}.$$

A direct scan with explicit Hilbert-square testing is far too slow near $10^{16}$, so the solution rewrites the condition as a Möbius-weighted sum over odd square divisors.

Mathematical Approach

The main job is to understand which ordinary prime exponents are allowed in a squarefree Hilbert number, and then encode that description with inclusion-exclusion.

Step 1: Translate the Hilbert-square condition into ordinary prime exponents

Write an odd number $n\equiv1\pmod4$ as

$$n=\prod_{p\equiv1\pmod4} p^{a_p}\prod_{q\equiv3\pmod4} q^{b_q}.$$

If some prime $p\equiv1\pmod4$ has $a_p\ge2$, then $p$ itself is a Hilbert number and $p^2\mid n$, so $n$ is not squarefree in the Hilbert sense.

For primes $q\equiv3\pmod4$, one copy of $q^2$ is still allowed because $q\not\equiv1\pmod4$. But two square units are too many. If

$$\sum_{q\equiv3\pmod4}\left\lfloor\frac{b_q}{2}\right\rfloor\ge2,$$

then we can build a nontrivial Hilbert number from those $3\bmod4$ primes: for example $q^2$ from one prime with $b_q\ge4$, or $qr$ from two different primes with $b_q,b_r\ge2$. Its square then divides $n$.

Therefore $n$ is squarefree Hilbert exactly when

$$a_p\le1\quad\text{for every }p\equiv1\pmod4,$$

$$\sum_{q\equiv3\pmod4}\left\lfloor\frac{b_q}{2}\right\rfloor\le1.$$

Step 2: Obtain a structural decomposition

The previous criterion has a concrete consequence. A squarefree Hilbert number is of one of two types:

$$\text{either }n\text{ is ordinary squarefree,}$$

$$\text{or }n=q^2m\text{ for exactly one prime }q\equiv3\pmod4\text{ and an ordinary squarefree }m.$$

The second case is unique. If two different $3\bmod4$ primes contributed square factors, or if one such prime contributed exponent at least $4$, then the sum of the floor terms above would be at least $2$, which is forbidden.

This explains examples such as $45=3^2\cdot5$, which is allowed, and $81=3^4$, which is not.

Step 3: Encode ordinary squarefreeness with the Möbius function

Let $\mathbf{1}_{\mathrm{sqf}}(m)$ be the indicator of ordinary squarefree integers. The standard identity is

$$\mathbf{1}_{\mathrm{sqf}}(m)=\sum_{u^2\mid m}\mu(u).$$

Using the decomposition from Step 2, the indicator of squarefree Hilbert numbers becomes

$$I(n)=\mathbf{1}_{n\equiv1\pmod4}\left(\mathbf{1}_{\mathrm{sqf}}(n)+\sum_{\substack{q\equiv3\pmod4\\ q\text{ prime}}}\mathbf{1}_{q^2\mid n}\,\mathbf{1}_{\mathrm{sqf}}\!\left(\frac{n}{q^2}\right)\right).$$

The first term counts numbers with no square factor at all. The second term restores the legal case where the only square part is one factor $q^2$ with $q\equiv3\pmod4$.

Step 4: Collect everything into one coefficient

Substitute the Möbius identity into both squarefree indicators:

$$I(n)=\mathbf{1}_{n\equiv1\pmod4}\left(\sum_{u^2\mid n}\mu(u)+\sum_{\substack{q\equiv3\pmod4\\ q\text{ prime}\\ q^2\mid n}}\ \sum_{u^2\mid n/q^2}\mu(u)\right).$$

In the second double sum, reindex with $v=qu$. Then all terms become contributions from odd squares $v^2\mid n$:

$$I(n)=\mathbf{1}_{n\equiv1\pmod4}\sum_{\substack{v^2\mid n\\ v\text{ odd}}} c(v),$$

where

$$c(v)=\mu(v)+\sum_{\substack{q\mid v\\ q\equiv3\pmod4\\ q\text{ prime}}}\mu\!\left(\frac{v}{q}\right).$$

This is the exact coefficient formula used by the implementations. The correction term cancels illegal square divisors such as $5^2$, keeps legal ones such as $3^2$, and subtracts higher powers such as $3^4$ again.

Step 5: Turn the indicator into a summatory formula

Let

$$A(t)=\#\{m\le t:m\equiv1\pmod4\}=\left\lfloor\frac{t+3}{4}\right\rfloor.$$

For any odd $v$, the condition $v^2\mid n$ with $n\equiv1\pmod4$ is equivalent to writing $n=v^2m$ with $m\equiv1\pmod4$, because $v^2\equiv1\pmod4$. Hence

$$S(N)=\sum_{\substack{v\ge1\\ v\text{ odd}}} c(v)\,A\!\left(\left\lfloor\frac{N}{v^2}\right\rfloor\right).$$

Only values $v\le\sqrt{N}$ contribute, so the sum is finite:

$$\boxed{S(N)=\sum_{\substack{1\le v\le\lfloor\sqrt{N}\rfloor\\ v\text{ odd}}} c(v)\left\lfloor\frac{\left\lfloor N/v^2\right\rfloor+3}{4}\right\rfloor.}$$

Worked Example: Why $S(99)=23$

There are exactly $25$ integers $n\le99$ with $n\equiv1\pmod4$. For this bound, only odd $v\le9$ matter.

The relevant coefficients are

$$c(1)=1,\qquad c(3)=0,\qquad c(5)=-1,\qquad c(7)=0,\qquad c(9)=-1.$$

Indeed, $c(5)=-1$ removes multiples of $25$, while $c(9)=-1$ removes multiples of $81=(3^2)^2$. The terms for $3$ and $7$ vanish because a single square $3^2$ or $7^2$ is still allowed in the Hilbert setting.

Therefore

$$S(99)=A(99)-A\!\left(\left\lfloor\frac{99}{25}\right\rfloor\right)-A\!\left(\left\lfloor\frac{99}{81}\right\rfloor\right)=25-1-1=23,$$

which matches the checkpoint count below $100$.

How the Code Works

The C++, Python, and Java implementations all use the same arithmetic plan. They first convert the exclusive input bound $X$ into the inclusive bound $N=X-1$, compute $U=\lfloor\sqrt{N}\rfloor$, and store only odd candidates up to $U$. That halves the sieve size immediately, because every relevant square root $v$ is odd.

Next, the implementation builds an odd-only prime sieve, derives the Möbius function $\mu(v)$ on those odd values, and records the primes congruent to $3\pmod4$. It then starts from the base coefficients $\mu(v)$ and adds the correction term $\mu(v/q)$ along odd multiples of each prime $q\equiv3\pmod4$, producing the final weights $c(v)$.

Finally, it evaluates

$$c(v)\left\lfloor\frac{\left\lfloor N/v^2\right\rfloor+3}{4}\right\rfloor$$

for every odd $v\le U$ and accumulates the result. The C++ version keeps the same formula but can split the final summation into several ranges when multithreading is enabled; the Python and Java versions run the same arithmetic in a single thread.

Complexity Analysis

Let $U=\lfloor\sqrt{N}\rfloor$. Building the odd-only sieve and the odd Möbius table costs $O(U\log\log U)$ time and $O(U)$ memory. Adding the correction terms over multiples of primes $q\equiv3\pmod4$ has the same harmonic-sum behavior, so it also stays within $O(U\log\log U)$ time.

The final accumulation is a single pass over the odd indices, so it is $O(U)$. Overall, the method is near-linear in $\sqrt{N}$ and uses $O(U)$ memory.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=580
Möbius function: Wikipedia — Möbius function
Squarefree integer: Wikipedia — Squarefree integer
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Modular arithmetic: Wikipedia — Modular arithmetic

Problem 580 source code

C++

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;

constexpr u64 kDefaultLimitExclusive = 10'000'000'000'000'000ULL;
constexpr u64 kSampleLimitExclusive = 10'000'000ULL;
constexpr u64 kSampleExpected = 2'327'192ULL;
constexpr u64 kBruteCheckpointLimitExclusive = 200'003ULL;

struct Options {
    u64 limit_exclusive = kDefaultLimitExclusive;
    bool allow_multithreading = true;
    bool run_checkpoints = true;
    unsigned requested_threads = 0;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& 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;
    }

    u64 parsed = 0;
    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_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }

    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        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;
        }

        u64 limit = 0;
        if (parse_u64_after_prefix(arg, "--limit=", limit)) {
            options.limit_exclusive = limit;
            continue;
        }

        unsigned threads = 0;
        if (parse_unsigned_after_prefix(arg, "--threads=", threads)) {
            options.requested_threads = threads;
            continue;
        }

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

    return true;
}

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1ULL) <= x / (r + 1ULL)) {
        ++r;
    }
    while (r > x / r) {
        --r;
    }
    return r;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2'000'000ULL) {
        return 1;
    }

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }

    return std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

struct SieveData {
    std::vector<std::int8_t> mobius_odd;
    std::vector<int> primes_3mod4;
};

SieveData build_sieve_data(u64 max_u) {
    const std::size_t odd_count = static_cast<std::size_t>(max_u / 2ULL + 1ULL);
    std::vector<std::uint8_t> composite(odd_count, 0U);
    std::vector<int> odd_primes;
    if (max_u >= 3ULL) {
        odd_primes.reserve(static_cast<std::size_t>(max_u / std::max<u64>(
                                                                1ULL,
                                                                static_cast<u64>(
                                                                    std::log(static_cast<long double>(max_u))))));
    }

    for (std::size_t idx = 1; idx < odd_count; ++idx) {
        if (composite[idx] != 0U) {
            continue;
        }

        const u64 p = 2ULL * static_cast<u64>(idx) + 1ULL;
        odd_primes.push_back(static_cast<int>(p));

        if (p > max_u / p) {
            continue;
        }

        std::size_t mark = static_cast<std::size_t>((p * p - 1ULL) / 2ULL);
        while (mark < odd_count) {
            composite[mark] = 1U;
            mark += static_cast<std::size_t>(p);
        }
    }

    SieveData out;
    out.mobius_odd.assign(odd_count, static_cast<std::int8_t>(1));
    out.primes_3mod4.reserve(odd_primes.size() / 2ULL + 1ULL);

    for (const int p_int : odd_primes) {
        const u64 p = static_cast<u64>(p_int);
        if ((p & 3ULL) == 3ULL) {
            out.primes_3mod4.push_back(p_int);
        }

        std::size_t idx = static_cast<std::size_t>((p - 1ULL) / 2ULL);
        while (idx < odd_count) {
            out.mobius_odd[idx] = static_cast<std::int8_t>(-out.mobius_odd[idx]);
            idx += static_cast<std::size_t>(p);
        }

        if (p > max_u / p) {
            continue;
        }
        const u64 p2 = p * p;
        idx = static_cast<std::size_t>((p2 - 1ULL) / 2ULL);
        while (idx < odd_count) {
            out.mobius_odd[idx] = static_cast<std::int8_t>(0);
            idx += static_cast<std::size_t>(p2);
        }
    }

    return out;
}

std::vector<std::int16_t> build_coefficients(const std::vector<std::int8_t>& mobius_odd,
                                             const std::vector<int>& primes_3mod4) {
    std::vector<std::int16_t> coeff(mobius_odd.size(), 0);
    for (std::size_t i = 0; i < mobius_odd.size(); ++i) {
        coeff[i] = static_cast<std::int16_t>(mobius_odd[i]);
    }

    // c(u) = mu(u) + sum_{q|u, q==3(mod4) prime} mu(u/q), for odd u.
    for (const int q_int : primes_3mod4) {
        const std::size_t q = static_cast<std::size_t>(q_int);
        std::size_t n_idx = (q - 1ULL) / 2ULL;
        std::size_t m_idx = 0;
        while (n_idx < coeff.size()) {
            coeff[n_idx] = static_cast<std::int16_t>(coeff[n_idx] + mobius_odd[m_idx]);
            ++m_idx;
            n_idx += q;
        }
    }

    return coeff;
}

i128 accumulate_range(const std::vector<std::int16_t>& coeff,
                     u64 inclusive_limit,
                     std::size_t begin,
                     std::size_t end) {
    i128 local = 0;
    for (std::size_t idx = begin; idx < end; ++idx) {
        const int c = static_cast<int>(coeff[idx]);
        if (c == 0) {
            continue;
        }
        const u64 u = 2ULL * static_cast<u64>(idx) + 1ULL;
        const u64 u2 = u * u;
        // Number of integers <= inclusive_limit / u^2 that are 1 (mod 4).
        const u64 terms = (inclusive_limit / u2 + 3ULL) / 4ULL;
        local += static_cast<i128>(c) * static_cast<i128>(terms);
    }
    return local;
}

u64 count_squarefree_hilbert_upto(u64 inclusive_limit,
                                  bool allow_multithreading,
                                  unsigned requested_threads) {
    const u64 max_u = isqrt_u64(inclusive_limit);
    const SieveData sieve = build_sieve_data(max_u);
    const std::vector<std::int16_t> coeff =
        build_coefficients(sieve.mobius_odd, sieve.primes_3mod4);

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, coeff.size());
    if (threads == 1) {
        const i128 total = accumulate_range(coeff, inclusive_limit, 0, coeff.size());
        return static_cast<u64>(total);
    }

    std::vector<std::thread> pool;
    std::vector<i128> partial(threads, 0);
    pool.reserve(threads);

    for (unsigned t = 0; t < threads; ++t) {
        const std::size_t begin = coeff.size() * t / threads;
        const std::size_t end = coeff.size() * (t + 1ULL) / threads;
        pool.emplace_back([&, t, begin, end]() {
            partial[t] = accumulate_range(coeff, inclusive_limit, begin, end);
        });
    }

    for (auto& th : pool) {
        th.join();
    }

    i128 total = 0;
    for (const i128 x : partial) {
        total += x;
    }

    return static_cast<u64>(total);
}

u64 count_squarefree_hilbert_below(u64 exclusive_limit,
                                   bool allow_multithreading,
                                   unsigned requested_threads) {
    if (exclusive_limit == 0ULL) {
        return 0ULL;
    }
    return count_squarefree_hilbert_upto(exclusive_limit - 1ULL,
                                         allow_multithreading,
                                         requested_threads);
}

u64 brute_force_count_below(u64 exclusive_limit) {
    if (exclusive_limit == 0ULL) {
        return 0ULL;
    }

    const u64 inclusive_limit = exclusive_limit - 1ULL;
    const u64 root = isqrt_u64(inclusive_limit);

    std::vector<u64> hilbert_squares;
    for (u64 h = 5ULL; h <= root; h += 4ULL) {
        hilbert_squares.push_back(h * h);
    }

    u64 count = 0;
    for (u64 n = 1ULL; n <= inclusive_limit; n += 4ULL) {
        bool ok = true;
        for (const u64 sq : hilbert_squares) {
            if (sq > n) {
                break;
            }
            if (n % sq == 0ULL) {
                ok = false;
                break;
            }
        }
        if (ok) {
            ++count;
        }
    }

    return count;
}

bool run_checkpoints(const Options& options) {
    struct Point {
        u64 limit_exclusive;
        u64 expected;
        const char* label;
    };

    const std::vector<Point> fixed = {
        {10ULL, 3ULL, "count(<10)"},
        {13ULL, 3ULL, "count(<13)"},
        {100ULL, 23ULL, "count(<100)"},
        {1'000ULL, 232ULL, "count(<1000)"},
        {kSampleLimitExclusive, kSampleExpected, "count(<10^7)"}};

    for (const Point& cp : fixed) {
        const u64 got =
            count_squarefree_hilbert_below(cp.limit_exclusive, false, 1U);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed: " << cp.label << ", expected " << cp.expected
                      << ", got " << got << '\n';
            return false;
        }
        std::cout << "Checkpoint OK: " << cp.label << " = " << got << '\n';
    }

    const u64 brute_expected = brute_force_count_below(kBruteCheckpointLimitExclusive);
    const u64 brute_fast = count_squarefree_hilbert_below(
        kBruteCheckpointLimitExclusive,
        options.allow_multithreading,
        options.requested_threads);
    if (brute_expected != brute_fast) {
        std::cerr << "Brute checkpoint failed at <" << kBruteCheckpointLimitExclusive
                  << ": expected " << brute_expected << ", got " << brute_fast << '\n';
        return false;
    }
    std::cout << "Checkpoint OK: brute cross-check <" << kBruteCheckpointLimitExclusive
              << " = " << brute_fast << '\n';

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const auto t0 = std::chrono::steady_clock::now();

    if (options.run_checkpoints) {
        if (!run_checkpoints(options)) {
            return 1;
        }
    }

    const u64 answer = count_squarefree_hilbert_below(options.limit_exclusive,
                                                      options.allow_multithreading,
                                                      options.requested_threads);

    const auto t1 = std::chrono::steady_clock::now();
    const std::chrono::duration<long double> elapsed = t1 - t0;

    std::cout << "Answer: " << answer << '\n';
    std::cout << "count(<" << options.limit_exclusive << ") = " << answer << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " s\n";

    return 0;
}

Python

import math

def solve():
    LIMIT = 10**16

    def isqrt(x):
        r=int(math.isqrt(x))
        while (r+1)*(r+1)<=x: r+=1
        while r*r>x: r-=1
        return r

    inc = LIMIT - 1; mu = isqrt(inc)
    odd_cnt = mu//2+1
    # Sieve odd primes
    composite = bytearray(odd_cnt)
    odd_primes = []
    for idx in range(1, odd_cnt):
        if composite[idx]: continue
        p = 2*idx+1; odd_primes.append(p)
        if p<=mu//p:
            mark=(p*p-1)//2
            while mark<odd_cnt: composite[mark]=1; mark+=p

    # Mobius on odd numbers
    mob = [1]*odd_cnt
    for p in odd_primes:
        idx=(p-1)//2
        while idx<odd_cnt: mob[idx]=-mob[idx]; idx+=p
        if p<=mu//p:
            p2=p*p; idx=(p2-1)//2
            while idx<odd_cnt: mob[idx]=0; idx+=p2

    # Primes 3 mod 4
    p3m4 = [p for p in odd_primes if p%4==3]

    # Build coefficients c(u) = mu(u) + sum_{q|u, q==3(4)} mu(u/q)
    coeff = list(mob)
    for q in p3m4:
        n_idx=(q-1)//2; m_idx=0
        while n_idx<odd_cnt: coeff[n_idx]+=mob[m_idx]; m_idx+=1; n_idx+=q

    # Accumulate
    total = 0
    for idx in range(odd_cnt):
        c = coeff[idx]
        if c==0: continue
        u = 2*idx+1; u2=u*u
        terms = (inc//u2+3)//4
        total += c*terms
    return str(total)

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

Java

public class Euler580 {
    public static String solve() {
        long LIMIT = (long) 1e16;
        long inc = LIMIT - 1;
        int mu = (int) Math.sqrt((double) inc);
        int oddCnt = mu / 2 + 1;
        boolean[] composite = new boolean[oddCnt];
        java.util.List<Integer> oddPrimes = new java.util.ArrayList<>();
        for (int idx = 1; idx < oddCnt; idx++) {
            if (composite[idx])
                continue;
            int p = 2 * idx + 1;
            oddPrimes.add(p);
            if (p <= mu / p) {
                int mark = (p * p - 1) / 2;
                while (mark < oddCnt) {
                    composite[mark] = true;
                    mark += p;
                }
            }
        }
        int[] mob = new int[oddCnt];
        java.util.Arrays.fill(mob, 1);
        for (int p : oddPrimes) {
            int idx = (p - 1) / 2;
            while (idx < oddCnt) {
                mob[idx] = -mob[idx];
                idx += p;
            }
            if (p <= mu / p) {
                int p2 = p * p;
                idx = (p2 - 1) / 2;
                while (idx < oddCnt) {
                    mob[idx] = 0;
                    idx += p2;
                }
            }
        }
        int[] coeff = mob.clone();
        for (int q : oddPrimes) {
            if (q % 4 != 3)
                continue;
            int nIdx = (q - 1) / 2, mIdx = 0;
            while (nIdx < oddCnt) {
                coeff[nIdx] += mob[mIdx];
                mIdx++;
                nIdx += q;
            }
        }
        long total = 0;
        for (int idx = 0; idx < oddCnt; idx++) {
            int c = coeff[idx];
            if (c == 0)
                continue;
            long u = 2L * idx + 1, u2 = u * u;
            long terms = (inc / u2 + 3) / 4;
            total += c * terms;
        }
        return String.valueOf(total);
    }

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