puzzle276.py

# -*- coding: utf-8 -*-
from __future__ import division, print_function
import itertools as it
import six

from pyprimes.factors import factorise, factors  # noqa

try:
    from math import gcd  # noqa
except ImportError:
    from fractions import gcd  # noqa

map, range, zip = six.moves.map, six.moves.range, six.moves.zip


def powerset(iterable):
    """powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"""
    s = list(iterable)
    return it.chain.from_iterable(it.combinations(s, r)
                                  for r in range(len(s) + 1))


def work(limit, c):
    total = 0
    for b in range(c // 2, min(c, limit - c - 1) + 1):
        gcd_bc = gcd(b, c)
        max_a = min((b, limit - b - c, limit // 3)) + 1
        min_a = max(1, c - b + 1)
        if gcd_bc == 1:
            total += max(0, max_a - min_a)
        for p in powerset(f for f, _ in factors(gcd_bc)):
            sign = -1 if len(p) % 2 else 1
            assert sign in (0, 1)
    return total


for x in range(1, 1000 // 2 + 1):
    work(1000, x)