"""
Martin Fuller, Sep 08 2023.
Program to calculate A111233: "Number of nonempty subsets of 
{1, 1/2, 1/3, ..., 1/n} that sum to an integer."

Theorem 1:
Let p be prime, e >= 1, and let x & y be fractions in lowest terms, where the
denominator of y is not divisible by p^e.
Then p^e divides the denominator of x iff p^e divides the denominator of x+y.

Proof:
Suppose p^e does not divide the denominator of x. Then p^e does not divide
lcm(denominator of x, denominator of y), and therefore does not divide the
denominator of x+y.

Suppose p^e does not divide the denominator of x+y. The denominator of (-y)
is the same as the denominator of y, and therefore is not divisible by p^e.
Therefore we can use the above result for (x+y) + (-y) = x.

Theorem 2:
Let S be a subset of {1/1, 1/2, ..., 1/n} such that sum(S) is an integer.
Let X be any subset of S such that S\X = Y contains no fractions whose
denominator is divisible by p^e, for some fixed prime p and e >= 1.
Then the denominator of sum(X) is not divisible by p^e.

Proof:
Let x = sum(X), y = sum(Y).
p^e does not divide the denominator of y, since it does not divide any
denominator in the set Y, and we can apply Theorem 1 repeatedly.
p^e does not divide the denominator of x+y, because x+y = sum(S) is integer.
Therefore p^e does not divide the denominator of x, by Theorem 1.

Theorem 3:
Let n = k*p^e, with 1 <= k <= floor(p/2), p prime, e >= 1. Working in Z[p],
if 0 is not the sum of any subset of {1/1, 1/2, ..., 1/k} which contains 1/k,
then a(n) = a(n-1).

Proof:
Let S be a subset of {1/1, 1/2, ..., 1/n} such that sum(S) is an integer.
Let X = {x in S | p^e divides denominator of x}. Let Y = S\X.
By Theorem 2 we know that p^e does not divide the denominator of sum(X).
Therefore p divides the numerator of sum(X)*p^e.
Let W = {(x*p^e) in Z[p] | x in X}.  W has the following properties:
  * W is well-defined and |W| = |X|, because n < p^(e+1).
  * W is a subset of {1/1, 1/2, ..., 1/k} in Z[p].
  * sum(W) = 0, because p divides the numerator of sum(X)*p^e. 
Using the initial assumption we conclude that W does not contain 1/k.
Therefore X and S do not contain 1/n, and all subsets counted by a(n)
are also counted by a(n-1).

Question: If the conditions for Theorem 3 do not hold, is a(n) > a(n-1)?
"""

from fractions import Fraction
from collections import Counter
from more_itertools import map_reduce, powerset

import itertools

def main():
    for n in itertools.count(1):
        if not is_a111233_fixed(n):
            a = a111233(n)
        print(f'{n} {a}')

def is_a111233_fixed(n:int) -> bool:
    """Checks the conditions for Theorem 3"""
    if n==1: return False
    p,e = largest_prime_factor(n)
    k = n // p**e
    if k > p//2: return False
    s = {pow(k, -1, p)}
    for i in (pow(i, -1, p) for i in range(1, k)):
        s |= {(i + j) % p for j in s}
    return 0 not in s

def a111233(n:int) -> int:
    """A111233: Number of nonempty subsets of {1, 1/2, 1/3, ..., 1/n}
    that sum to an integer.

    This implementation uses Counter dictionaries to combine equal sums and
    keep track of the corresponding number of subsets.
    Integer parts of each sum are discarded.

    Method:
    Group the numbers {2, 3, ... n} in descending order of their largest
    prime factor and its exponent.
    n=12: [11], [7], [5,10], [9], [3,6,12], [8], [4], [2]
    Note: Numbers divisible by p^e cannot be after the corresponding group.
    
    For each group in order:
      * Generate all subset sums from reciprocals in the group.
      * Partition above group sums according to (s*p^e) in Z[p].
      * Generate new partial sums from the old partial sums (initially {0})
        plus each group sum. Use Z[p] to ensure the denominators are not
        divisible by p^e, because only these sums can lead to an integer.
    
    Finally, double the total (subsets with {1}) and subtract 1 (empty set).
    """
    groups = map_reduce(range(2,n+1), largest_prime_factor)
    groups = sorted(groups.items(), reverse=True)

    fraction_modp = FractionModp()

    partial_sums = Counter([Fraction(0)])
    for (p,e),g in groups:
        power = p**e
        fraction_modp.setp(p)

        group_sums = Counter([0])
        for k in g:
            fk = Fraction(1,k)
            group_sums += {(f+fk)%1:c for f,c in group_sums.items()}

        group_sums_modp = [Counter() for _ in range(p)]
        for s,c in group_sums.items():
            modp = fraction_modp.convert(s * power)
            group_sums_modp[modp][s] += c

        new_sums = Counter()
        for f,c in partial_sums.items():
            modp = fraction_modp.convert(f * -power)
            for fg,cg in group_sums_modp[modp].items():
                new_sums[(f+fg)%1] += c*cg

        partial_sums = new_sums
    return partial_sums.total() * 2 - 1

class FractionModp:
    """Class to convert a Fraction object into an integer in Z[p].

    Example: 2/5 in Z[7] == 6, because (5*6) % 7 == 2.
    """
    def __init__(self):
        self.p = 0

    def setp(self, p:int):
        if self.p != p:
            self.p = p
            self.inverse = [0] * p

    def convert(self, f:Fraction) -> int:
        return (f.numerator * self.invert(f.denominator)) % self.p

    def invert(self, n:int) -> int:
        n = n % self.p
        if self.inverse[n] == 0:
            self.inverse[n] = pow(n, -1, self.p)
        return self.inverse[n]
    
def largest_prime_factor(n:int)->tuple[int,int]:
    """Return (p,e) for the largest prime p dividing n."""
    return factorization(n)[-1]

def factorization(n:int)->list[tuple[int,int]]:
    """Returns prime factorization of n as [(p,e),...]"""
    if n <= 0:
        raise ValueError(n)
    f = []
    if n == 1: return f
    for p in itertools.chain([2], range(3,n+1,2)):
        if p*p > n:
            f.append((n,1))
            return f
        e = 0
        while n % p == 0:
            e += 1
            n = n // p
        if  e:
            f.append((p,e))
            if n==1:
                return f

main()