特殊号码计数

Question

它是一个gcd（各位数字的四次方之和，各位数字的乘积）大于1的数。 例如。 123 是一个特殊数字，因为(1+16+81, 6) 的 hcf 大于 1。

我必须找到低于输入 n 的所有这些数字的计数。 例如。 对于n=120，它们是（1和120）之间的57个特殊数字，

我已经编写了一个代码，但它非常慢，你能告诉我以某种又好又快的方式来做吗？ 有没有什么方法可以使用一些数学来做到这一点。

import math,numpy
t = int(input())

ans = []

for i in range(0,t):
    ans.append(0)
    n = int(input())
    for j in range(1, n+1):
        res = math.gcd(sum([pow(int(k),4) for k in str(j)]),numpy.prod([int(k) for k in str(j)]))
        if res>1:
            ans[i] = ans[i] + 1

for i in range(0,t):
    print(ans[i])

Answer 1

关键的观察是特殊数字的十进制表示构成了常规语言。下面是 Python 中的有限状态识别器。本质上，我们跟踪乘积的质因数（gcd > 1 相当于有一个共同的质因数）和幂总和 mod 2×3×5×7 的余数，以及一些状态处理涉及零的边缘情况。

从那里，我们可以构造一个显式自动机，然后使用动态规划来计算值小于 n 的接受字符串的数量。

def is_special(digits):
    greater_than_zero = False
    greater_than_one = False
    prod_zero = False
    prod_two = False
    mod_two = 0
    prod_three = False
    mod_three = 0
    prod_five = False
    mod_five = 0
    prod_seven = False
    mod_seven = 0
    for d in digits:
        if d > 1:
            greater_than_one = True
        elif d == 1:
            if greater_than_zero:
                greater_than_one = True
            else:
                greater_than_zero = True
        if d == 0:
            prod_zero = True
        if d % 2 == 0:
            prod_two = True
        mod_two = (mod_two + d ** 4) % 2
        if d % 3 == 0:
            prod_three = True
        mod_three = (mod_three + d ** 4) % 3
        if d % 5 == 0:
            prod_five = True
        mod_five = (mod_five + d ** 4) % 5
        if d % 7 == 0:
            prod_seven = True
        mod_seven = (mod_seven + d ** 4) % 7
    return (
        greater_than_one
        if prod_zero
        else (
            (prod_two and mod_two == 0)
            or (prod_three and mod_three == 0)
            or (prod_five and mod_five == 0)
            or (prod_seven and mod_seven == 0)
        )
    )


# Test case
import functools
import math
import operator


def digits_of(n):
    return [int(digit) for digit in str(n)]


def is_special_reference(digits):
    power_sum = sum(digit ** 4 for digit in digits)
    product = functools.reduce(operator.mul, digits, 1)
    return math.gcd(power_sum, product) > 1


def test():
    for n in range(10000):
        digits = digits_of(n)
        assert is_special(digits) == is_special_reference(digits), str(n)


if __name__ == "__main__":
    test()

Answer 2

这是一个 O(log n) 算法，用于实际计算小于或等于 n 的特殊数字。它一次构建一个数字字符串，跟踪该数字字符串的乘积是否被 2、3、5 和 7 整除，以及这些数字的四次方总和对 2、3、5 和 7 取模的余数。

根据这些质因数和这些因数的余数的整除性来测试一个数字是否特殊的逻辑与大卫的答案相同，并且在那里得到了更好的解释。由于素数除积的可能性只有 2^4 种，而余数的可能性为 2*3*5*7 种，因此存在恒定数量的组合对于 O(2^4 * 210 * log n) = O(log n) 的运行时间，这两种情况都是可能的。

def count_special_less_equal(digits: List[int]) -> int:
    """Return the count of numbers less than or equal to the represented
    number, with the property that
    gcd(product(digits), sum(fourth powers of digits)) > 1"""

    # Count all digit strings with zeroes
    total_non_special = len(digits)

    primes = (2, 3, 5, 7)
    prime_product = functools.reduce(operator.mul, primes, 1)

    digit_to_remainders = [pow(x, 4, prime_product) for x in range(10)]

    # Map each digit 1-9 to prime factors
    # 2: 2**0, 3: 2**1, 5: 2**2, 7: 2**3
    factor_masks = [0, 0, 1, 2, 1, 4, 3, 8, 1, 2]

    def is_fac_mask_mod_special(factor_mask: int,
                                remainder: int) -> bool:
        """Return true if any of the prime factors represented in factor_mask
        have corresponding remainder 0 (i.e., divide the sum of fourth powers)"""

        return any((factor_mask & (1 << i) != 0
                    and remainder % primes[i] == 0)
                   for i in range(4))

    prefix_less_than = [Counter() for _ in range(16)]

    # Empty string
    prefix_equal = (0, 0)

    for digit_pos, digit in enumerate(digits):

        new_prefix_less_than = [Counter() for _ in range(16)]

        # Old "lesser than" prefixes stay lesser
        for fac_mask, fac_mask_counts in enumerate(prefix_less_than):
            for new_digit in range(1, 10):
                new_mask = fac_mask | factor_masks[new_digit]
                remainder_change = digit_to_remainders[new_digit]
                for old_remainder, old_count in fac_mask_counts.items():
                    new_remainder = (remainder_change + old_remainder) % prime_product
                    new_prefix_less_than[new_mask][new_remainder] += old_count

        if digit == 0:
            prefix_equal = None

        if prefix_equal is not None:
            equal_fac_mask, equal_remainder = prefix_equal
            for new_digit in range(1, digit):

                new_mask = equal_fac_mask | factor_masks[new_digit]

                remainder_change = digit_to_remainders[new_digit]
                new_remainder = (remainder_change + equal_remainder) % prime_product

                new_prefix_less_than[new_mask][new_remainder] += 1

            new_mask = equal_fac_mask | factor_masks[digit]
            remainder_change = digit_to_remainders[digit]

            new_remainder = (remainder_change + equal_remainder) % prime_product
            prefix_equal = (new_mask, new_remainder)

        prefix_less_than = new_prefix_less_than

        if digit_pos == len(digits) - 1:
            break

        # Empty string
        prefix_less_than[0][0] += 1

    for fac_mask, fac_mask_counts in enumerate(prefix_less_than):
        for remainder, rem_count in fac_mask_counts.items():
            if not is_fac_mask_mod_special(factor_mask=fac_mask,
                                           remainder=remainder):
                total_non_special += rem_count

    if prefix_equal is not None:
        if not is_fac_mask_mod_special(*prefix_equal):
            total_non_special += 1

    return 1 + int(''.join(map(str, digits))) - total_non_special

用法示例：

print(f"{count_special_less_equal(digits_of(120))}")

打印

57

并

for exponent in range(1, 19):
    print(f"Count up to 10^{exponent}: {count_special_less_equal(digits_of(10**exponent))}")

给出：

Count up to 10^1: 8
Count up to 10^2: 42
Count up to 10^3: 592
Count up to 10^4: 7400
Count up to 10^5: 79118
Count up to 10^6: 854190
Count up to 10^7: 8595966
Count up to 10^8: 86010590
Count up to 10^9: 866103492
Count up to 10^10: 8811619132
Count up to 10^11: 92967009216
Count up to 10^12: 929455398976
Count up to 10^13: 9268803096820
Count up to 10^14: 92838342330554
Count up to 10^15: 933105194955392
Count up to 10^16: 9557298732021784
Count up to 10^17: 96089228976983058
Count up to 10^18: 960712913414545906

Done in 0.3783 seconds

这将在大约三分之一秒内找到 10 到 10^18 的所有幂的频率。可以使用 numpy 数组或其他技巧（例如预先计算具有固定位数的所有数字的计数）在常数因子中进一步优化。