Problem 768: Chandelier

Project Euler Problem 768 Solution

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

Problem Summary Let $\zeta_n=e^{2\pi i/n}$. The quantity $f(n,m)$ counts the $m$-element subsets $A\subseteq\{0,1,\dots,n-1\}$ for which $$\sum_{j\in A}\zeta_n^j=0.$$ Geometrically, we choose $m$ equally weighted points among the $n$ equally spaced directions on the unit circle and ask when their vector sum is exactly zero. Directly testing all $\binom{n}{m}$ subsets is hopeless for the real target, so the implementations reorganize the problem into a smaller primitive cycle and then rebuild the full answer from a generating function. Mathematical Approach Write $$S=\operatorname{rad}(n),\qquad G=\frac{n}{S},$$ where $\operatorname{rad}(n)$ is the product of the distinct prime divisors of $n$. The core fact used by the implementations is that all repeated prime powers can be separated cleanly, so the hard combinatorics happens only on the square-free part $S$. Step 1: Encode a Selection by a Polynomial For a subset $A\subseteq\{0,\dots,n-1\}$, define the indicator polynomial $$F_A(x)=\sum_{j=0}^{n-1}\varepsilon_j x^j,\qquad \varepsilon_j\in\{0,1\},$$ with $\varepsilon_j=1$ exactly when $j\in A$. Then $|A|=\sum \varepsilon_j$, and the zero-sum condition is simply $$F_A(\zeta_n)=0.$$ So the problem is not about floating-point geometry at all; it is about deciding when a $0/1$-polynomial vanishes at a primitive $n$-th root of unity....

Detailed mathematical approach

Problem Summary

Let $\zeta_n=e^{2\pi i/n}$. The quantity $f(n,m)$ counts the $m$-element subsets $A\subseteq\{0,1,\dots,n-1\}$ for which

$$\sum_{j\in A}\zeta_n^j=0.$$

Geometrically, we choose $m$ equally weighted points among the $n$ equally spaced directions on the unit circle and ask when their vector sum is exactly zero. Directly testing all $\binom{n}{m}$ subsets is hopeless for the real target, so the implementations reorganize the problem into a smaller primitive cycle and then rebuild the full answer from a generating function.

Mathematical Approach

Write

$$S=\operatorname{rad}(n),\qquad G=\frac{n}{S},$$

where $\operatorname{rad}(n)$ is the product of the distinct prime divisors of $n$. The core fact used by the implementations is that all repeated prime powers can be separated cleanly, so the hard combinatorics happens only on the square-free part $S$.

Step 1: Encode a Selection by a Polynomial

For a subset $A\subseteq\{0,\dots,n-1\}$, define the indicator polynomial

$$F_A(x)=\sum_{j=0}^{n-1}\varepsilon_j x^j,\qquad \varepsilon_j\in\{0,1\},$$

with $\varepsilon_j=1$ exactly when $j\in A$. Then $|A|=\sum \varepsilon_j$, and the zero-sum condition is simply

$$F_A(\zeta_n)=0.$$

So the problem is not about floating-point geometry at all; it is about deciding when a $0/1$-polynomial vanishes at a primitive $n$-th root of unity.

Step 2: Split the $n$-Cycle into $G$ Primitive Blocks

Group exponents by their residue modulo $G$:

$$F_A(x)=\sum_{r=0}^{G-1}x^r B_r(x^G),$$

where

$$B_r(y)=\sum_{q=0}^{S-1}\varepsilon_{r+qG}y^q.$$

Because $\zeta_n^G=\zeta_S$, evaluation at $\zeta_n$ gives

$$F_A(\zeta_n)=\sum_{r=0}^{G-1}\zeta_n^r B_r(\zeta_S).$$

Now $n$ and $S$ have the same prime divisors, so

$$[\mathbb{Q}(\zeta_n):\mathbb{Q}(\zeta_S)]=\frac{\varphi(n)}{\varphi(S)}=\frac{n}{S}=G.$$

Therefore $1,\zeta_n,\zeta_n^2,\dots,\zeta_n^{G-1}$ are linearly independent over $\mathbb{Q}(\zeta_S)$. The sum above is zero if and only if every block vanishes separately:

$$F_A(\zeta_n)=0\iff B_r(\zeta_S)=0\quad\text{for all }r=0,\dots,G-1.$$

Each residue class modulo $G$ is thus an independent rotated copy of an $S$-cycle. This is the decisive reduction used by the code.

Step 3: Count Zero-Sum Subsets on One $S$-Cycle

Let $c_k$ be the number of $k$-element subsets of $\{0,\dots,S-1\}$ whose $S$-th roots sum to zero. For one such subset $T$, define

$$F_T(x)=\sum_{q\in T}x^q.$$

Since $\zeta_S$ is a primitive $S$-th root of unity, its minimal polynomial over $\mathbb{Q}$ is $\Phi_S(x)$. Hence

$$F_T(\zeta_S)=0\iff \Phi_S(x)\mid F_T(x).$$

So the problem becomes: among all $0/1$-polynomials of degree $\lt S$, count how many become the zero element in the quotient ring

$$\mathbb{Z}[x]/(\Phi_S(x)).$$

The implementations represent each power $x^q$ by its coefficient vector modulo $\Phi_S(x)$, whose dimension is $\varphi(S)$. A subset is valid exactly when those vectors add up to the zero vector.

Step 4: Meet-in-the-Middle for All $c_k$

The $S$ exponents are split into two halves. For each left-half subset, the implementation computes its vector sum in $\mathbb{Z}[x]/(\Phi_S(x))$, records its size, and stores how many times that vector occurs. Then it enumerates right-half subsets, negates their vector sum, and looks up matches from the left side. If a left subset of size $a$ and a right subset of size $b$ have opposite vectors, their union is a zero-sum subset of size $a+b$.

One meet-in-the-middle pass therefore produces every $c_k$ for $0\le k\le m$, not just one cardinality. This is why the later polynomial stage can be exact and inexpensive.

Step 5: Rebuild the Full Count with a Generating Function

For one primitive block define the ordinary generating polynomial

$$P_S(t)=\sum_{k=0}^{S}c_k t^k.$$

Because the $G$ residue classes are independent, choosing a valid subset from the full $n$-cycle is the same as choosing, for each block, a zero-sum subset on the $S$-cycle. Sizes add, so the global generating function is

$$P_S(t)^G.$$

Therefore

$$\boxed{f(n,m)=\left[t^m\right]P_S(t)^G.}$$

This boxed formula is exactly what the C++, Python, and Java implementations evaluate.

Worked Example: $f(12,4)$

Here $n=12$, so $S=\operatorname{rad}(12)=6$ and $G=12/6=2$. On one 6-cycle, the zero-sum subsets are easy to classify geometrically:

$$c_0=1,\qquad c_2=3,\qquad c_3=2,\qquad c_4=3,\qquad c_6=1,$$

and all other $c_k$ are $0$. Thus

$$P_6(t)=1+3t^2+2t^3+3t^4+t^6.$$

We only need the coefficient of $t^4$, so terms above degree $4$ can be ignored:

$$f(12,4)=\left[t^4\right](1+3t^2+2t^3+3t^4)^2.$$

The $t^4$ contribution comes from $1\cdot 3t^4$, $3t^4\cdot 1$, and $3t^2\cdot 3t^2$. Hence

$$f(12,4)=3+3+9=15,$$

which matches the small checkpoint used by the implementations.

How the Code Works

The implementation first computes $S=\operatorname{rad}(n)$ and $G=n/S$. It then builds the cyclotomic polynomial $\Phi_S(x)$ from the identity $x^k-1=\prod_{d\mid k}\Phi_d(x)$, reducing powers of $x$ to coefficient vectors modulo $\Phi_S(x)$. With those vectors available, it performs the meet-in-the-middle count described above and obtains all coefficients $c_0,c_1,\dots,c_m$ for one primitive block.

Next it forms the polynomial $P_S(t)$ and multiplies it by itself $G$ times, truncating degrees above $m$ after every multiplication because only $\left[t^m\right]$ matters. All arithmetic is exact: the Python version relies on native arbitrary-precision integers, the Java version uses arbitrary-precision integer objects, and the C++ version stores large integers explicitly so no modulus ever hides carries or cancellations.

Complexity Analysis

Let $d=\varphi(S)$. Building and using the quotient-ring representation costs polynomial time in $S$ and $d$, which is small compared with the subset search. The meet-in-the-middle stage enumerates about $2^{S/2}$ subsets on each side and processes vectors of length $d$, so its practical cost is roughly $O(2^{S/2}d)$ hash-based work with $O(2^{S/2})$ stored states. The polynomial stage performs $G$ truncated convolutions up to degree $m$, giving $O(Gm^2)$ exact coefficient operations; the bit-cost depends on the size of the final integer. For the target case $n=360$, the crucial win is that the hard subset stage depends on $S=\operatorname{rad}(360)=30$, not on all $360$ positions.

Footnotes and References

Problem page: https://projecteuler.net/problem=768
Roots of unity: Wikipedia — Root of unity
Cyclotomic polynomials: Wikipedia — Cyclotomic polynomial
Generating functions: Wikipedia — Generating function
Meet-in-the-middle: Wikipedia — Meet-in-the-middle

Problem 768 source code

C++

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <thread>
#include <unordered_map>
#include <vector>

using namespace std;

namespace {

struct BigInt {
    static constexpr uint32_t kBase = 1000000000U;
    vector<uint32_t> digits;

    BigInt(uint64_t value = 0) {
        if (value == 0) {
            digits.push_back(0);
            return;
        }
        while (value > 0) {
            digits.push_back(static_cast<uint32_t>(value % kBase));
            value /= kBase;
        }
    }

    bool is_zero() const {
        return digits.size() == 1 && digits[0] == 0;
    }

    void trim() {
        while (digits.size() > 1 && digits.back() == 0) {
            digits.pop_back();
        }
    }

    void add(const BigInt& other) {
        uint64_t carry = 0;
        size_t n = max(digits.size(), other.digits.size());
        if (digits.size() < n) digits.resize(n, 0);
        for (size_t i = 0; i < n; ++i) {
            uint64_t sum = carry + digits[i];
            if (i < other.digits.size()) sum += other.digits[i];
            digits[i] = static_cast<uint32_t>(sum % kBase);
            carry = sum / kBase;
        }
        if (carry) digits.push_back(static_cast<uint32_t>(carry));
    }

    BigInt operator*(const BigInt& other) const {
        if (is_zero() || other.is_zero()) return BigInt(0);
        BigInt res;
        res.digits.assign(digits.size() + other.digits.size(), 0);
        for (size_t i = 0; i < digits.size(); ++i) {
            uint64_t carry = 0;
            for (size_t j = 0; j < other.digits.size() || carry; ++j) {
                uint64_t cur = res.digits[i + j] + carry;
                if (j < other.digits.size()) {
                    cur += static_cast<uint64_t>(digits[i]) * other.digits[j];
                }
                res.digits[i + j] = static_cast<uint32_t>(cur % kBase);
                carry = cur / kBase;
            }
        }
        res.trim();
        return res;
    }

    bool equals_uint64(uint64_t value) const {
        BigInt other(value);
        return digits == other.digits;
    }
};

ostream& operator<<(ostream& os, const BigInt& value) {
    if (value.digits.empty()) return os << 0;
    os << value.digits.back();
    for (int i = static_cast<int>(value.digits.size()) - 2; i >= 0; --i) {
        os << setw(9) << setfill('0') << value.digits[i];
    }
    return os;
}

struct Poly {
    vector<BigInt> coeffs;

    explicit Poly(int max_deg = 0) : coeffs(max_deg + 1, BigInt(0)) {}

    static Poly one() {
        Poly p(0);
        p.coeffs[0] = BigInt(1);
        return p;
    }

    Poly multiply(const Poly& other, int max_deg) const {
        int deg_a = static_cast<int>(coeffs.size()) - 1;
        int deg_b = static_cast<int>(other.coeffs.size()) - 1;
        int new_deg = min(max_deg, deg_a + deg_b);
        Poly res(new_deg);
        for (int i = 0; i <= deg_a; ++i) {
            if (coeffs[i].is_zero()) continue;
            for (int j = 0; j <= deg_b && i + j <= new_deg; ++j) {
                if (other.coeffs[j].is_zero()) continue;
                BigInt prod = coeffs[i] * other.coeffs[j];
                res.coeffs[i + j].add(prod);
            }
        }
        return res;
    }
};

vector<long long> trim_poly(vector<long long> poly) {
    while (poly.size() > 1 && poly.back() == 0) {
        poly.pop_back();
    }
    return poly;
}

vector<long long> divide_poly(vector<long long> a, const vector<long long>& b) {
    int da = static_cast<int>(a.size()) - 1;
    int db = static_cast<int>(b.size()) - 1;
    vector<long long> q(max(0, da - db + 1), 0);
    for (int i = da; i >= db; --i) {
        long long coef = a[i];
        if (coef == 0) continue;
        q[i - db] = coef;
        for (int j = 0; j <= db; ++j) {
            a[i - db + j] -= coef * b[j];
        }
    }
    return trim_poly(q);
}

vector<long long> cyclotomic(int n) {
    vector<vector<long long>> phi(n + 1);
    phi[1] = {-1, 1};
    for (int k = 2; k <= n; ++k) {
        vector<long long> poly(k + 1, 0);
        poly[0] = -1;
        poly[k] = 1;
        for (int d = 1; d < k; ++d) {
            if (k % d == 0) {
                poly = divide_poly(poly, phi[d]);
            }
        }
        phi[k] = trim_poly(poly);
    }
    return phi[n];
}

vector<int> multiply_by_x(const vector<int>& v, const vector<int>& rel) {
    int d = static_cast<int>(v.size());
    vector<int> res(d, 0);
    for (int i = 0; i + 1 < d; ++i) {
        res[i + 1] += v[i];
    }
    int h = v[d - 1];
    if (h != 0) {
        for (int i = 0; i < d; ++i) {
            res[i] -= h * rel[i];
        }
    }
    return res;
}

struct VectorHash {
    size_t operator()(const vector<int>& v) const noexcept {
        uint64_t h = 0xcbf29ce484222325ULL;
        for (int x : v) {
            uint64_t y = static_cast<uint64_t>(static_cast<int64_t>(x) ^ (static_cast<int64_t>(x) >> 32));
            h ^= y + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        }
        return static_cast<size_t>(h);
    }
};

vector<uint64_t> count_zero_sums(int S, int m) {
    vector<long long> phi = cyclotomic(S);
    int d = static_cast<int>(phi.size()) - 1;
    vector<int> rel(d, 0);
    for (int i = 0; i < d; ++i) rel[i] = static_cast<int>(phi[i]);

    vector<vector<int>> powers(S, vector<int>(d, 0));
    powers[0][0] = 1;
    for (int i = 1; i < S; ++i) {
        powers[i] = multiply_by_x(powers[i - 1], rel);
    }

    int mid = S / 2;
    int right_n = S - mid;
    uint64_t left_masks = 1ULL << mid;
    uint64_t right_masks = 1ULL << right_n;

    unordered_map<vector<int>, vector<uint32_t>, VectorHash> left_sums;
    left_sums.reserve(left_masks * 2);

    vector<int> sum(d, 0);
    for (uint64_t mask = 0; mask < left_masks; ++mask) {
        fill(sum.begin(), sum.end(), 0);
        int cnt = 0;
        for (int k = 0; k < mid; ++k) {
            if (mask & (1ULL << k)) {
                ++cnt;
                const auto& pv = powers[k];
                for (int i = 0; i < d; ++i) sum[i] += pv[i];
            }
        }
        if (cnt > m) continue;
        auto it = left_sums.find(sum);
        if (it == left_sums.end()) {
            it = left_sums.emplace(sum, vector<uint32_t>(m + 1, 0)).first;
        }
        it->second[cnt] += 1;
    }

    unsigned hw = thread::hardware_concurrency();
    int threads = hw == 0 ? 1 : static_cast<int>(hw);
    if (right_masks < 4096) threads = 1;
    if (threads > 8) threads = 8;
    if (right_masks < static_cast<uint64_t>(threads)) {
        threads = static_cast<int>(right_masks);
    }

    vector<vector<uint64_t>> partial(threads, vector<uint64_t>(m + 1, 0));
    vector<thread> workers;
    uint64_t chunk = (right_masks + threads - 1) / threads;

    auto worker = [&](int idx, uint64_t start, uint64_t end) {
        vector<int> local_sum(d, 0);
        auto& local = partial[idx];
        for (uint64_t mask = start; mask < end; ++mask) {
            fill(local_sum.begin(), local_sum.end(), 0);
            int cnt = 0;
            for (int k = 0; k < right_n; ++k) {
                if (mask & (1ULL << k)) {
                    ++cnt;
                    const auto& pv = powers[mid + k];
                    for (int i = 0; i < d; ++i) local_sum[i] += pv[i];
                }
            }
            if (cnt > m) continue;
            for (int i = 0; i < d; ++i) local_sum[i] = -local_sum[i];
            auto it = left_sums.find(local_sum);
            if (it == left_sums.end()) continue;
            const auto& counts = it->second;
            for (int lc = 0; lc + cnt <= m; ++lc) {
                uint32_t ways = counts[lc];
                if (ways) local[lc + cnt] += ways;
            }
        }
    };

    for (int t = 0; t < threads; ++t) {
        uint64_t start = static_cast<uint64_t>(t) * chunk;
        uint64_t end = min(right_masks, start + chunk);
        workers.emplace_back(worker, t, start, end);
    }
    for (auto& th : workers) th.join();

    vector<uint64_t> total(m + 1, 0);
    for (int t = 0; t < threads; ++t) {
        for (int k = 0; k <= m; ++k) {
            total[k] += partial[t][k];
        }
    }
    return total;
}

int radical(int n) {
    int result = 1;
    int x = n;
    for (int p = 2; p * p <= x; ++p) {
        if (x % p == 0) {
            result *= p;
            while (x % p == 0) x /= p;
        }
    }
    if (x > 1) result *= x;
    return result;
}

BigInt solve_chandelier(int n, int m) {
    if (m < 0 || m > n) return BigInt(0);
    if (m == 0) return BigInt(1);

    int S = radical(n);
    int G = n / S;

    vector<uint64_t> counts = count_zero_sums(S, m);
    Poly P(m);
    for (int k = 0; k <= m; ++k) P.coeffs[k] = BigInt(counts[k]);

    Poly Final = Poly::one();
    for (int i = 0; i < G; ++i) {
        Final = Final.multiply(P, m);
    }
    return m < static_cast<int>(Final.coeffs.size()) ? Final.coeffs[m] : BigInt(0);
}

}  // namespace

int main() {
    BigInt v12 = solve_chandelier(12, 4);
    cout << "f(12, 4)   = " << v12;
    cout << (v12.equals_uint64(15) ? " [PASS]" : " [FAIL]") << '\n';

    BigInt v36 = solve_chandelier(36, 6);
    cout << "f(36, 6)   = " << v36;
    cout << (v36.equals_uint64(876) ? " [PASS]" : " [FAIL]") << '\n';

    cout << "f(360, 20) = " << solve_chandelier(360, 20) << '\n';
    return 0;
}

Python

def trim_poly(poly):
    while len(poly) > 1 and poly[-1] == 0:
        poly.pop()
    return poly

def divide_poly(a, b):
    da = len(a) - 1
    db = len(b) - 1
    q = [0] * max(0, da - db + 1)
    
    a = list(a) 
    for i in range(da, db - 1, -1):
        coef = a[i]
        if coef == 0:
            continue
        q[i - db] = coef
        for j in range(db + 1):
            a[i - db + j] -= coef * b[j]
    return trim_poly(q)

def cyclotomic(n):
    phi = [[] for _ in range(n + 1)]
    phi[1] = [-1, 1]
    for k in range(2, n + 1):
        poly = [0] * (k + 1)
        poly[0] = -1
        poly[k] = 1
        for d in range(1, k):
            if k % d == 0:
                poly = divide_poly(poly, phi[d])
        phi[k] = trim_poly(poly)
    return phi[n]

def multiply_by_x(v, rel):
    d = len(v)
    res = [0] * d
    for i in range(d - 1):
        res[i + 1] += v[i]
    h = v[d - 1]
    if h != 0:
        for i in range(d):
            res[i] -= h * rel[i]
    return res

def count_zero_sums(S, m):
    phi = cyclotomic(S)
    d = len(phi) - 1
    rel = phi[:d]
    
    powers = [[0] * d for _ in range(S)]
    powers[0][0] = 1
    for i in range(1, S):
        powers[i] = multiply_by_x(powers[i - 1], rel)
        
    mid = S // 2
    right_n = S - mid
    
    left_sums = {}
    
    for mask in range(1 << mid):
        cnt = bin(mask).count('1')
        if cnt > m:
            continue
        
        sum_arr = [0] * d
        for k in range(mid):
            if (mask >> k) & 1:
                pv = powers[k]
                for i in range(d):
                    sum_arr[i] += pv[i]
                    
        key = tuple(sum_arr)
        if key not in left_sums:
            left_sums[key] = [0] * (m + 1)
        left_sums[key][cnt] += 1
        
    total = [0] * (m + 1)
    
    for mask in range(1 << right_n):
        cnt = bin(mask).count('1')
        if cnt > m:
            continue
            
        sum_arr = [0] * d
        for k in range(right_n):
            if (mask >> k) & 1:
                pv = powers[mid + k]
                for i in range(d):
                    sum_arr[i] -= pv[i]
                    
        key = tuple(sum_arr)
        counts = left_sums.get(key)
        if counts:
            for lc in range(m - cnt + 1):
                ways = counts[lc]
                if ways:
                    total[lc + cnt] += ways
                    
    return total

def radical(n):
    res = 1
    x = n
    p = 2
    while p * p <= x:
        if x % p == 0:
            res *= p
            while x % p == 0:
                x //= p
        p += 1
    if x > 1:
        res *= x
    return res

def multiply_poly(A, B, max_deg):
    res = [0] * (min(max_deg, len(A) + len(B) - 2) + 1)
    for i, a in enumerate(A):
        if a == 0: continue
        for j, b in enumerate(B):
            if b == 0: continue
            if i + j <= max_deg:
                res[i + j] += a * b
    return res

def solve_chandelier(n, m):
    if m < 0 or m > n: return 0
    if m == 0: return 1
    
    S = radical(n)
    G = n // S
    
    counts = count_zero_sums(S, m)
    P = counts[:(m + 1)]
    
    Final = [1]
    for _ in range(G):
        Final = multiply_poly(Final, P, m)
        
    return Final[m] if m < len(Final) else 0

def solve():
    ans = solve_chandelier(360, 20)
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;

public class Euler768 {

    static ArrayList<Long> trimPoly(ArrayList<Long> poly) {
        while (poly.size() > 1 && poly.get(poly.size() - 1) == 0L) {
            poly.remove(poly.size() - 1);
        }
        return poly;
    }

    static ArrayList<Long> dividePoly(ArrayList<Long> aList, ArrayList<Long> bList) {
        long[] a = new long[aList.size()];
        for (int i = 0; i < a.length; i++)
            a[i] = aList.get(i);
        long[] b = new long[bList.size()];
        for (int i = 0; i < b.length; i++)
            b[i] = bList.get(i);

        int da = a.length - 1;
        int db = b.length - 1;
        long[] q = new long[Math.max(0, da - db + 1)];

        for (int i = da; i >= db; --i) {
            long coef = a[i];
            if (coef == 0)
                continue;
            q[i - db] = coef;
            for (int j = 0; j <= db; ++j) {
                a[i - db + j] -= coef * b[j];
            }
        }

        ArrayList<Long> res = new ArrayList<>();
        for (long val : q)
            res.add(val);
        return trimPoly(res);
    }

    static ArrayList<Long> cyclotomic(int n) {
        ArrayList<ArrayList<Long>> phi = new ArrayList<>(n + 1);
        for (int i = 0; i <= n; i++)
            phi.add(new ArrayList<>());

        phi.get(1).add(-1L);
        phi.get(1).add(1L);

        for (int k = 2; k <= n; ++k) {
            ArrayList<Long> poly = new ArrayList<>();
            for (int i = 0; i <= k; i++)
                poly.add(0L);
            poly.set(0, -1L);
            poly.set(k, 1L);

            for (int d = 1; d < k; ++d) {
                if (k % d == 0) {
                    poly = dividePoly(poly, phi.get(d));
                }
            }
            phi.set(k, trimPoly(poly));
        }
        return phi.get(n);
    }

    static int[] multiplyByX(int[] v, int[] rel) {
        int d = v.length;
        int[] res = new int[d];
        for (int i = 0; i + 1 < d; ++i) {
            res[i + 1] += v[i];
        }
        int h = v[d - 1];
        if (h != 0) {
            for (int i = 0; i < d; ++i) {
                res[i] -= h * rel[i];
            }
        }
        return res;
    }

    static long[] countZeroSums(int S, int m) {
        ArrayList<Long> phi = cyclotomic(S);
        int d = phi.size() - 1;
        int[] rel = new int[d];
        for (int i = 0; i < d; ++i)
            rel[i] = phi.get(i).intValue();

        int[][] powers = new int[S][d];
        powers[0][0] = 1;
        for (int i = 1; i < S; ++i) {
            powers[i] = multiplyByX(powers[i - 1], rel);
        }

        int mid = S / 2;
        int rightN = S - mid;

        HashMap<Long, int[]> leftSums = new HashMap<>(1 << mid);

        for (int mask = 0; mask < (1 << mid); ++mask) {
            int cnt = Integer.bitCount(mask);
            if (cnt > m)
                continue;

            long key = 0;
            int[] sumArr = new int[d];
            for (int k = 0; k < mid; ++k) {
                if (((mask >> k) & 1) != 0) {
                    for (int i = 0; i < d; ++i)
                        sumArr[i] += powers[k][i];
                }
            }

            for (int i = 0; i < d; ++i) {
                long val = sumArr[i] + 30; // offset by 30 to keep it positive
                key = (key << 8) | val;
            }

            int[] counts = leftSums.get(key);
            if (counts == null) {
                counts = new int[m + 1];
                leftSums.put(key, counts);
            }
            counts[cnt]++;
        }

        long[] total = new long[m + 1];

        for (int mask = 0; mask < (1 << rightN); ++mask) {
            int cnt = Integer.bitCount(mask);
            if (cnt > m)
                continue;

            long key = 0;
            int[] sumArr = new int[d];
            for (int k = 0; k < rightN; ++k) {
                if (((mask >> k) & 1) != 0) {
                    for (int i = 0; i < d; ++i)
                        sumArr[i] -= powers[mid + k][i];
                }
            }

            for (int i = 0; i < d; ++i) {
                long val = sumArr[i] + 30;
                key = (key << 8) | val;
            }

            int[] counts = leftSums.get(key);
            if (counts != null) {
                for (int lc = 0; lc + cnt <= m; ++lc) {
                    int ways = counts[lc];
                    if (ways > 0) {
                        total[lc + cnt] += ways;
                    }
                }
            }
        }
        return total;
    }

    static int radical(int n) {
        int res = 1;
        int x = n;
        for (int p = 2; p * p <= x; ++p) {
            if (x % p == 0) {
                res *= p;
                while (x % p == 0)
                    x /= p;
            }
        }
        if (x > 1)
            res *= x;
        return res;
    }

    static BigInteger[] multiplyPoly(BigInteger[] a, BigInteger[] b, int maxDeg) {
        int newDeg = Math.min(maxDeg, a.length + b.length - 2);
        BigInteger[] res = new BigInteger[newDeg + 1];
        Arrays.fill(res, BigInteger.ZERO);

        for (int i = 0; i < a.length; ++i) {
            if (a[i].equals(BigInteger.ZERO))
                continue;
            for (int j = 0; j < b.length; ++j) {
                if (b[j].equals(BigInteger.ZERO))
                    continue;
                if (i + j <= maxDeg) {
                    res[i + j] = res[i + j].add(a[i].multiply(b[j]));
                }
            }
        }
        return res;
    }

    static BigInteger solveChandelier(int n, int m) {
        if (m < 0 || m > n)
            return BigInteger.ZERO;
        if (m == 0)
            return BigInteger.ONE;

        int S = radical(n);
        int G = n / S;

        long[] counts = countZeroSums(S, m);
        BigInteger[] P = new BigInteger[m + 1];
        for (int k = 0; k <= m; ++k)
            P[k] = BigInteger.valueOf(counts[k]);

        BigInteger[] Final = new BigInteger[] { BigInteger.ONE };
        for (int i = 0; i < G; ++i) {
            Final = multiplyPoly(Final, P, m);
        }

        return (m < Final.length) ? Final[m] : BigInteger.ZERO;
    }

    public static String solve() {
        BigInteger ans = solveChandelier(360, 20);
        return ans.toString();
    }

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