Special Number Count

The critical observation is that the decimal representations of special numbers constitute a regular language. Below is a finite-state recognizer in Python. Essentially we track the prime factors of the product (gcd > 1 being equivalent to having a prime factor in common) and the residue of the sum of powers mod 2×3×5×7, as well as a little bit of state to handle edge cases involving zeros.

From there, we can construct an explicit automaton and then count the number of accepting strings whose value is less than n using dynamic programming.

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()

Here's an O(log n) algorithm for actually counting special numbers less than or equal to n. It builds digit strings one at a time, keeping track of whether 2, 3, 5 and 7 divide that digit string's product, and the remainder modulo 2, 3, 5, and 7 of the sum of fourth powers of those digits.

The logic for testing whether a number is special based on divisibility by those prime factors and remainder of powers under those factors is the same as in David's answer, and is explained better there. Since there are only 2^4 possibilities for which primes divide the product, and 2*3*5*7 possibilities for the remainder, there are a constant number of combinations of both that are possible, for a runtime of 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

Example usage:

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

prints

57

and

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

gives:

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

This finds the frequencies for all powers of 10 up to 10^18 in about a third of a second. It's possible to optimize this further in the constant factors, using numpy arrays or other tricks (like precomputing the counts for all numbers with a fixed number of digits).