Question: It is a number whose gcd of (sum of quartic power of its digits, the product of its digits) is more than 1. eg. 123 is a special number because hcf of(1+16+81, 6) is more than 1. I have to find the count of all these numbers that are below input n. eg. for n=120 their are 57 special numbers between (1 and 120) I have done a code but its very slow can you please tell me to do it in some good and fast way. Is there is any way to do it using some maths. 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: 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).