Check if num is anti-prime effectively

Viewed 176

I am trying to write a predicate which checks if num is anti-prime.

A highly composite number (sometimes referred to as 'anti-prime') is a positive integer with more divisors than any smaller positive integer has.

I came up with this.

def count_divisors(num):
    divisor_count = 0
    for i in range(1, int(sqrt(num)) + 1):
        if num % i == 0:
            if num / i == i:
                divisor_count += 1
                continue
            divisor_count += 2
    return divisor_count


def is_highly_composite(x):
    original_divisors = count_divisors(x)
    for i in range(x-1, 0, -1):
        if count_divisors(i) >= original_divisors:
            return False
    return True

Sadly it seems really inefficient and slow. Especially for big numbers.

I am learning coding so efficiency is kinda hard to wrap my head around.

Thanks for any help.

2 Answers

I could think of a bit neater approach to write your function count_divisors(). See the new version named count_divisors_np() below. (gmpy.is_square() allows for fast and accurate calculation for very large numbers, but you may just use divs[-1]**2 == num instead.)

import numpy as np
import gmpy2 as gp

def count_divisors_np(num):
    divs = np.arange(1, round(np.sqrt(num)) + 1)
    return 2*np.sum(num % divs == 0) - (1 if gp.is_square(num) else 0)

I compared the speeds, and on my PC, the new functions works consistently faster than the old one for numbers above 100,000, while it is almost twice as fast for numbers larger than 250,000. (See the code below that I used to compare the runtimes.)

import time

''' Other imports and definitions of count_divisors() and count_divisors_np()'''

def is_highly_composite(counter_func, x):
    original_divisors = counter_func(x)
    for i in range(x-1, 0, -1):
        if counter_func(i) >= original_divisors:
            return False
    return True


for i in range(1, 10**3):
    tic0 = time.time()
    for x in range(i*1000 + 1, (i+1)*1000+1):
        is_highly_composite(count_divisors, x)
    toc0 = time.time() - tic0
    tic1 = time.time()
    for x in range(i* 1000 + 1, (i+1) * 1000+1):
        is_highly_composite(count_divisors_np, x)
    toc1 = time.time() - tic1
    if toc1 < toc0:
        print('from', i*1000+1, 'to', (i+1)*1000, ':', toc1, '<', toc0)

You can get good performance decomposing the number into its prime factors, applying the rules for a highly composite number.

  • All primes must be used sequentially at least once
  • The last prime factor must have a power of 1
  • The powers of prime factors must be in non-increasing order
  • If all of the above, check for other combinations of prime factors that can produce a higher (or equal) number of divisors with a smaller product.
  • Special cases: 4 and 36 are anti-primes

Infinite primes generator (you can make your own as needed, or just use a hard coded list of the first 100 primes):

primes = [2,3]
skips  = {9:6}
def getPrimes(count=-1):
    yield from primes[:max(0,count) or None]
    p = primes[-1]+2
    while count:
        if p in skips:
            mult,step = p,skips.pop(p)
            mult += step
            while mult in skips: mult += step
            skips[mult] = step
        else:
            skips[p*p] = 2*p
            primes.append(p)
            yield p
            count -= 1
        p += 2

Anti-prime checking function:

from math import log
def isAntiPrime(N):
    if N in [1,4,36]: return True # special cases
    primeGen  = getPrimes()      # infinite primes generator
    factors   = []               # prime factors 
    lastPower = N                # to check non-ascending order 
    divCount  = 1                # number of divisors
    R         = N                # reduced N (by factors)
    for prime in primeGen:       # go through prime factors
        factors.append(prime)
        if R==1: break           # until number exausted 
        if R%prime: return False # unused prime factor
        power = 0                # Count power of prime factor
        while R%prime == 0:
            power += 1
            R    //= prime       # reduce number as we go
        if lastPower<power: return False # increasing order
        lastPower = power
        divCount *= power+1      # compute number of divisors
    if lastPower != 1: return False # last prime's power must be 1

    def canReduce(i=0,prod=1,count=1): # recursively look for smaller number
        if count>=divCount and prod<N: return True  # found one
        while i>=len(factors):                      # load primes list
            factors.append(next(primeGen))          # up to index 
        prime = factors[i]                          # factor at index
        for n in range(1,int(log(N/prod,prime))+1): # try prime^n
            prod *= prime
            if canReduce(i+1,prod,count*(n+1)):     # combined
                return True

    return not canReduce() # no better power combo -> anti-prime

Output:

i = 0
for n in range(1,1000000):
    if isAntiPrime(n):
        i += 1
        print(f"{i:2}",n)

 1 1
 2 2
 3 4
 4 6
 5 12
 6 24
 7 36
 8 48
 9 60
10 120
11 180
12 240
13 360
14 720
15 840
16 1260
17 1680
18 2520
19 5040
20 7560
21 10080
22 15120
23 20160
24 25200
25 27720
26 45360
27 50400
28 55440
29 83160
30 110880
31 166320
32 221760
33 277200
34 332640
35 498960
36 554400
37 665280
38 720720

Other tests:

for n in [166320,166322,498960,494851, 169733893,
          180068938933390833890254309768093,
          7, 244877512201, 51983595331405913494380649, 
          2473017995717899, 252222613252687371148099,
          7426227904486559070625807, 24247699, 3463957,
          3674876304763078242658251219757, 713928607, 70693, 16807, 
          24975008738755062001]:
    print(isAntiPrime(n),n)

True 166320
False 166322
True 498960
False 494851
False 169733893
False 180068938933390833890254309768093
False 7
False 244877512201
False 51983595331405913494380649
False 2473017995717899
False 252222613252687371148099
False 7426227904486559070625807
False 24247699
False 3463957
False 3674876304763078242658251219757
False 713928607
False 70693
False 16807
False 24975008738755062001

Note I didn't include any timing benchmarks because this is all instantaneous and just with the first 100 primes you can process numbers with up to 190 digits in less than a second

Finding anti-primes

The isAntiPrime function can be adapted to return the largest anti-prime <=N by leveraging its canReduce() function, turning it into an goal seeker from which we can get the largest divisor count with the smallest corresponding number:

def antiPrime(N):
    if N in [1,4,36]: return N   # special cases
    primeGen  = getPrimes()      # infinite primes generator
    factors   = []               # prime factors 
    divCount  = 1                # number of divisors
    R         = N                # reduced N (by factors)
    for prime in primeGen:       # go through prime factors
        factors.append(prime)
        if R==1: break           # until number exausted
        if prime*prime>R:        # when past square root,
            prime = R            #   R is last prime
        power = 0                # Count power of prime factor
        while R%prime == 0:
            power += 1
            R    //= prime       # reduce number as we go
        divCount *= power+1      # compute number of divisors

    betterN = 0 
    def reduce(i=0,prod=1,count=1,lastPower=N): # look for smaller number
        nonlocal divCount,betterN
        if prod>=N: return
        if count>divCount or count==divCount and prod<betterN:
            divCount,betterN = count,prod           # found better
        while i>=len(factors):                      # load primes list
            factors.append(next(primeGen))          # up to index 
        n,prime = 1,factors[i]                      # factor at index
        while n<=lastPower and prod*prime<N:        # try prime^n       
            prod *= prime
            reduce(i+1,prod,count*(n+1),n) # combined
            n    += 1
            
    reduce()
    return betterN or N

Output:

for n in (1000000,713928607,2473017995717899,24975008738755062001,
         180068938933390833890254309768093):
    print(n,"anti-prime:",antiPrime(n))

1000000 anti-prime: 720720
713928607 anti-prime: 698377680
2473017995717899 anti-prime: 2021649740510400
24975008738755062001 anti-prime: 18401055938125660800
180068938933390833890254309768093 anti-prime: 156839524845080402008057229856000

This one is a bit slower on very large numbers because of the combinatory nature of the process.

Related