特殊号码计数

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\xc3\x973\xc3\x975\xc3\x977 的余数，以及一点状态来处理涉及零的边缘情况。

\n

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

\n

def is_special(digits):\n    greater_than_zero = False\n    greater_than_one = False\n    prod_zero = False\n    prod_two = False\n    mod_two = 0\n    prod_three = False\n    mod_three = 0\n    prod_five = False\n    mod_five = 0\n    prod_seven = False\n    mod_seven = 0\n    for d in digits:\n        if d > 1:\n            greater_than_one = True\n        elif d == 1:\n            if greater_than_zero:\n                greater_than_one = True\n            else:\n                greater_than_zero = True\n        if d == 0:\n            prod_zero = True\n        if d % 2 == 0:\n            prod_two = True\n        mod_two = (mod_two + d ** 4) % 2\n        if d % 3 == 0:\n            prod_three = True\n        mod_three = (mod_three + d ** 4) % 3\n        if d % 5 == 0:\n            prod_five = True\n        mod_five = (mod_five + d ** 4) % 5\n        if d % 7 == 0:\n            prod_seven = True\n        mod_seven = (mod_seven + d ** 4) % 7\n    return (\n        greater_than_one\n        if prod_zero\n        else (\n            (prod_two and mod_two == 0)\n            or (prod_three and mod_three == 0)\n            or (prod_five and mod_five == 0)\n            or (prod_seven and mod_seven == 0)\n        )\n    )\n\n\n# Test case\nimport functools\nimport math\nimport operator\n\n\ndef digits_of(n):\n    return [int(digit) for digit in str(n)]\n\n\ndef is_special_reference(digits):\n    power_sum = sum(digit ** 4 for digit in digits)\n    product = functools.reduce(operator.mul, digits, 1)\n    return math.gcd(power_sum, product) > 1\n\n\ndef test():\n    for n in range(10000):\n        digits = digits_of(n)\n        assert is_special(digits) == is_special_reference(digits), str(n)\n\n\nif __name__ == "__main__":\n    test()\n

\n

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