<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

200. Cracking RSA
time limit per test: 0.25 sec.
memory limit per test: 65536 KB
input: standard
output: standard



The following problem is somehow related to the final stage of many famous integer factorization algorithms involved in some cryptoanalytical problems, for example cracking well-known RSA public key system.

The most powerful of such algorithms, so called quadratic sieve descendant algorithms, utilize the fact that if n = pq where p and q are large unknown primes needed to be found out, then if v2=w2 (mod n), u ≠ v (mod n) and u ≠ -v (mod n), then gcd(v + w, n) is a factor of n (either p or q).

Not getting further in the details of these algorithms, let us consider our problem. Given m integer numbers b1, b2, ..., bm such that all their prime factors are from the set of first t primes, the task is to find such a subset S of {1, 2, ..., m} that product of bi for i from S is a perfect square i.e. equal to u2 for some integer u. Given such S we get one pair for testing (product of S elements stands for v when w is known from other steps of algorithms which are of no interest to us, testing performed is checking whether pair is nontrivial, i.e. u ≠ v (mod n) and u ≠ -v (mod n)). Since we want to factor n with maximum possible probability, we would like to get as many such sets as possible. So the interesting problem could be to calculate the number of all such sets. This is exactly your task.

Input

The first line of the input file contains two integers t and m (1 ≤ t ≤ 100, 1 ≤ m ≤ 100). The second line of the input file contains m integer numbers bi such that all their prime factors are from t first primes (for example, if t = 3 all their prime factors are from the set {2, 3, 5}). 1 ≤ bi ≤ 109 for all i.

Output

Output the number of non-empty subsets of the given set {bi}, the product of numbers from which is a perfect square


Sample test(s)

Input

3 4
9 20 500 3

Output

3
Author:	Andrew Stankevich
Resource:	Petrozavodsk Winter Trainings 2003
Date:	2003-02-06

<|response|>
1. Abridged problem statement
   Given t (≤100) and m (≤100), and m positive integers b₁…bₘ whose prime factors lie among the first t primes, count the number of non-empty subsets S of {1…m} such that the product of bᵢ for i∈S is a perfect square. Output the exact count.

2. Key observations
   • A positive integer is a perfect square if and only if every prime in its factorization has an even exponent.
   • For each bᵢ, record a t-dimensional vector over GF(2) of the parity (0 or 1) of the exponent of each of the first t primes.
   • A subset S multiplies to a square exactly when the bitwise XOR (sum mod 2) of its vectors is the zero vector.
   • If we build an m×t matrix A whose rows are those vectors, the number of solutions x∈{0,1}ᵐ to A·x=0 over GF(2) is 2^(m−rank(A)).
   • Excluding the trivial all-zero choice (empty subset), the answer is 2^(m−rank(A))−1. This can be as large as 2^100−1, so a big integer is required.

3. Full solution approach
   1. Generate the first t primes by simple trial division.
   2. For each bᵢ, factor out each of those primes and record the exponent mod 2 into a vector vᵢ of length t.
   3. Assemble these m vectors as rows of an m×t binary matrix.
   4. Perform Gaussian elimination over GF(2) on that matrix to compute its rank r.
      - Iterate columns 0…t−1, for each find a pivot row (with a 1 in that column), swap it into position “current rank”, then XOR-eliminate that column from all other rows.
   5. The null‐space dimension is m−r, so there are 2^(m−r) solutions; subtract 1 to exclude the empty subset.
   6. Compute 2^(m−r) with a big integer (here the bigint struct) by doubling 1 a total of m−r times, then subtract 1, and print the result.

4. C++ implementation
```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> z;
    int sign;

    bigint() : sign(1) {}

    bigint(long long v) { *this = v; }

    bigint(const string& s) { read(s); }

    void operator=(const bigint& v) {
        sign = v.sign;
        z = v.z;
    }

    void operator=(long long v) {
        sign = 1;
        if(v < 0) {
            sign = -1, v = -v;
        }
        z.clear();
        for(; v > 0; v = v / base) {
            z.push_back(v % base);
        }
    }

    bigint operator+(const bigint& v) const {
        if(sign == v.sign) {
            bigint res = v;

            for(int i = 0, carry = 0;
                i < (int)max(z.size(), v.z.size()) || carry; ++i) {
                if(i == (int)res.z.size()) {
                    res.z.push_back(0);
                }
                res.z[i] += carry + (i < (int)z.size() ? z[i] : 0);
                carry = res.z[i] >= base;
                if(carry) {
                    res.z[i] -= base;
                }
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint& v) const {
        if(sign == v.sign) {
            if(abs() >= v.abs()) {
                bigint res = *this;
                for(int i = 0, carry = 0; i < (int)v.z.size() || carry; ++i) {
                    res.z[i] -= carry + (i < (int)v.z.size() ? v.z[i] : 0);
                    carry = res.z[i] < 0;
                    if(carry) {
                        res.z[i] += base;
                    }
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = 0, carry = 0; i < (int)z.size() || carry; ++i) {
            if(i == (int)z.size()) {
                z.push_back(0);
            }
            long long cur = z[i] * (long long)v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
            // asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur),
            // "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) {
        int norm = base / (b1.z.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.z.resize(a.z.size());

        for(int i = a.z.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.z[i];
            int s1 = b.z.size() < r.z.size() ? r.z[b.z.size()] : 0;
            int s2 = b.z.size() - 1 < r.z.size() ? r.z[b.z.size() - 1] : 0;
            int d = ((long long)s1 * base + s2) / b.z.back();
            r -= b * d;
            while(r < 0) {
                r += b, --d;
            }
            q.z[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    friend bigint sqrt(const bigint& a1) {
        bigint a = a1;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        int n = a.z.size();

        int firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int norm = base / (firstDigit + 1);
        a *= norm;
        a *= norm;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        bigint r = (long long)a.z[n - 1] * base + a.z[n - 2];
        firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int q = firstDigit;
        bigint res;

        for(int j = n / 2 - 1; j >= 0; j--) {
            for(;; --q) {
                bigint r1 =
                    (r - (res * 2 * base + q) * q) * base * base +
                    (j > 0 ? (long long)a.z[2 * j - 1] * base + a.z[2 * j - 2]
                           : 0);
                if(r1 >= 0) {
                    r = r1;
                    break;
                }
            }
            res *= base;
            res += q;

            if(j > 0) {
                int d1 =
                    res.z.size() + 2 < r.z.size() ? r.z[res.z.size() + 2] : 0;
                int d2 =
                    res.z.size() + 1 < r.z.size() ? r.z[res.z.size() + 1] : 0;
                int d3 = res.z.size() < r.z.size() ? r.z[res.z.size()] : 0;
                q = ((long long)d1 * base * base + (long long)d2 * base + d3) /
                    (firstDigit * 2);
            }
        }

        res.trim();
        return res / norm;
    }

    bigint operator/(const bigint& v) const { return divmod(*this, v).first; }

    bigint operator%(const bigint& v) const { return divmod(*this, v).second; }

    void operator/=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = (int)z.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = z[i] + rem * (long long)base;
            z[i] = (int)(cur / v);
            rem = (int)(cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if(v < 0) {
            v = -v;
        }
        int m = 0;
        for(int i = z.size() - 1; i >= 0; --i) {
            m = (z[i] + m * (long long)base) % v;
        }
        return m * sign;
    }

    void operator+=(const bigint& v) { *this = *this + v; }
    void operator-=(const bigint& v) { *this = *this - v; }
    void operator*=(const bigint& v) { *this = *this * v; }
    void operator/=(const bigint& v) { *this = *this / v; }

    bool operator<(const bigint& v) const {
        if(sign != v.sign) {
            return sign < v.sign;
        }
        if(z.size() != v.z.size()) {
            return z.size() * sign < v.z.size() * v.sign;
        }
        for(int i = z.size() - 1; i >= 0; i--) {
            if(z[i] != v.z[i]) {
                return z[i] * sign < v.z[i] * sign;
            }
        }
        return false;
    }

    bool operator>(const bigint& v) const { return v < *this; }
    bool operator<=(const bigint& v) const { return !(v < *this); }
    bool operator>=(const bigint& v) const { return !(*this < v); }
    bool operator==(const bigint& v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint& v) const { return *this < v || v < *this; }

    void trim() {
        while(!z.empty() && z.back() == 0) {
            z.pop_back();
        }
        if(z.empty()) {
            sign = 1;
        }
    }

    bool isZero() const { return z.empty() || (z.size() == 1 && !z[0]); }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for(int i = z.size() - 1; i >= 0; i--) {
            res = res * base + z[i];
        }
        return res * sign;
    }

    friend bigint gcd(const bigint& a, const bigint& b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint& a, const bigint& b) {
        return a / gcd(a, b) * b;
    }

    void read(const string& s) {
        sign = 1;
        z.clear();
        int pos = 0;
        while(pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if(s[pos] == '-') {
                sign = -sign;
            }
            ++pos;
        }
        for(int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for(int j = max(pos, i - base_digits + 1); j <= i; j++) {
                x = x * 10 + s[j] - '0';
            }
            z.push_back(x);
        }
        trim();
    }

    friend istream& operator>>(istream& stream, bigint& v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream& stream, const bigint& v) {
        if(v.sign == -1) {
            stream << '-';
        }
        stream << (v.z.empty() ? 0 : v.z.back());
        for(int i = (int)v.z.size() - 2; i >= 0; --i) {
            stream << setw(base_digits) << setfill('0') << v.z[i];
        }
        return stream;
    }

    static vector<int> convert_base(
        const vector<int>& a, int old_digits, int new_digits
    ) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for(int i = 1; i < (int)p.size(); i++) {
            p[i] = p[i - 1] * 10;
        }
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for(int i = 0; i < (int)a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while(cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int)cur);
        while(!res.empty() && res.back() == 0) {
            res.pop_back();
        }
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll& a, const vll& b) {
        int n = a.size();
        vll res(n + n);
        if(n <= 32) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    res[i + j] += a[i] * b[j];
                }
            }
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for(int i = 0; i < k; i++) {
            a2[i] += a1[i];
        }
        for(int i = 0; i < k; i++) {
            b2[i] += b1[i];
        }

        vll r = karatsubaMultiply(a2, b2);
        for(int i = 0; i < (int)a1b1.size(); i++) {
            r[i] -= a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            r[i] -= a2b2[i];
        }

        for(int i = 0; i < (int)r.size(); i++) {
            res[i + k] += r[i];
        }
        for(int i = 0; i < (int)a1b1.size(); i++) {
            res[i] += a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            res[i + n] += a2b2[i];
        }
        return res;
    }

    bigint operator*(const bigint& v) const {
        vector<int> a6 = convert_base(this->z, base_digits, 6);
        vector<int> b6 = convert_base(v.z, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while(a.size() < b.size()) {
            a.push_back(0);
        }
        while(b.size() < a.size()) {
            b.push_back(0);
        }
        while(a.size() & (a.size() - 1)) {
            a.push_back(0), b.push_back(0);
        }
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for(int i = 0, carry = 0; i < (int)c.size(); i++) {
            long long cur = c[i] + carry;
            res.z.push_back((int)(cur % 1000000));
            carry = (int)(cur / 1000000);
        }
        res.z = convert_base(res.z, 6, base_digits);
        res.trim();
        return res;
    }
};

int t, m;
vector<int> b;

void read() {
    cin >> t >> m;
    b.resize(m);
    cin >> b;
}

void solve() {
    // - A product of a subset is a perfect square iff every prime exponent in
    //   the product is even. With prime factors drawn from the first t primes,
    //   each number maps to a t-bit vector over GF(2) recording the parity of
    //   each prime's exponent; a subset is a square iff its vectors XOR to zero.
    //
    // - The number of zero-XOR subsets is 2^(m - r), where r is the rank of the
    //   m vectors over GF(2), computed by Gaussian elimination. Subtract one to
    //   drop the empty subset.
    //
    // - m can be 100, so 2^(m - r) is produced as a big integer by repeated
    //   doubling before subtracting one.

    vector<int> primes;
    for(int p = 2; (int)primes.size() < t; p += (p == 2 ? 1 : 2)) {
        bool is_prime = true;
        for(int x: primes) {
            if((int64_t)x * x > p) {
                break;
            }
            if(p % x == 0) {
                is_prime = false;
                break;
            }
        }

        if(is_prime) {
            primes.push_back(p);
        }
    }

    vector<vector<int>> vectors(m, vector<int>(t, 0));
    for(int idx = 0; idx < m; idx++) {
        int tmp = b[idx];
        for(int i = 0; i < t; i++) {
            while(tmp % primes[i] == 0) {
                vectors[idx][i] ^= 1;
                tmp /= primes[i];
            }
            if(tmp == 1) {
                break;
            }
        }
    }

    int r = 0;
    for(int i = 0; i < t; i++) {
        int pivot = -1;
        for(int j = r; j < m; j++) {
            if(vectors[j][i] == 1) {
                pivot = j;
                break;
            }
        }

        if(pivot < 0) {
            continue;
        }

        swap(vectors[r], vectors[pivot]);
        for(int k = 0; k < m; k++) {
            if(k != r && vectors[k][i] == 1) {
                for(int c = 0; c < t; c++) {
                    vectors[k][c] ^= vectors[r][c];
                }
            }
        }

        r++;
    }

    bigint ans = 1;
    for(int i = 0; i < m - r; i++) {
        ans *= 2;
    }

    ans -= bigint(1);
    cout << ans << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

5. Python implementation with detailed comments
```python
import sys

def first_primes(t):
    """Return a list of the first t prime numbers."""
    primes = []
    p = 2
    while len(primes) < t:
        is_p = True
        for q in primes:
            if q*q > p: break
            if p % q == 0:
                is_p = False
                break
        if is_p:
            primes.append(p)
        p = 3 if p==2 else p+2
    return primes

def rank_gf2(mat, t):
    """
    Perform in-place Gaussian elimination on mat (list of m rows, each length t)
    over GF(2). Return the rank (number of pivot rows).
    """
    m = len(mat)
    r = 0  # current pivot row
    for col in range(t):
        # find a row >= r with a 1 in this column
        pivot = -1
        for row in range(r, m):
            if mat[row][col] == 1:
                pivot = row
                break
        if pivot < 0:
            continue
        # swap pivot row into position r
        mat[r], mat[pivot] = mat[pivot], mat[r]
        # eliminate this column from all other rows
        for row in range(m):
            if row != r and mat[row][col] == 1:
                # row_i ^= pivot_row
                mat[row] = [ (a ^ b) for a,b in zip(mat[row], mat[r]) ]
        r += 1
        if r == m:
            break
    return r

def main():
    data = sys.stdin.read().split()
    t, m = map(int, data[:2])
    b_list = list(map(int, data[2:]))

    # 1. generate primes
    primes = first_primes(t)

    # 2. build exponent-parity matrix
    mat = []
    for x in b_list:
        row = [0]*t
        tmp = x
        for i,p in enumerate(primes):
            while tmp % p == 0:
                row[i] ^= 1
                tmp //= p
            if tmp == 1:
                break
        mat.append(row)

    # 3. compute rank over GF(2)
    r = rank_gf2(mat, t)

    # 4. answer = 2^(m-r) - 1
    # Python has big ints, so use shift
    free_vars = m - r
    result = (1 << free_vars) - 1
    print(result)

if __name__ == "__main__":
    main()
```
