检查 num 是否有效地反素数

Question

我正在尝试编写一个谓词来检查 num 是否是反素数。

高度合数（有时称为“反素数”）是一个正整数，其除数比任何较小的正整数都多。

我想出了这个。

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

可悲的是，它似乎真的效率低下而且速度慢。尤其是对于大数字。

我正在学习编码，所以效率有点难以理解。

感谢您的帮助。

Answer 1

我可以想出一个更简洁的方法来编写你的函数count_divisors()。请参阅下面名为 count_divisors_np() 的新版本。 （gmpy.is_square() 允许对非常大的数字进行快速准确的计算，但您可以改用divs[-1]**2 == num。）

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)

我比较了速度，在我的 PC 上，对于 100,000 以上的数字，新功能的运行速度始终比旧功能快，而对于超过 250,000 的数字，它的速度几乎是旧功能的两倍。 （请参阅下面我用来比较运行时的代码。）

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)

Answer 2

应用highly composite number 的规则，您可以将数字分解为其主要因素，从而获得良好的性能。

• 所有素数必须至少按顺序使用一次
• 最后一个素数的幂必须是 1
• 质因数的幂必须是非递增顺序
• 如果满足以上所有条件，请检查其他素因数组合，这些组合可以产生更多（或相等）除数和更小的乘积。
• 特殊情况：4 和 36 是反素数

无限素数生成器（您可以根据需要制作自己的，或者只使用前 100 个素数的硬编码列表）：

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

反素检查功能：

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

输出：

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

其他测试：

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

请注意，我没有包含任何计时基准，因为这都是瞬时的，只需使用前 100 个素数，您就可以在不到一秒的时间内处理多达 190 位的数字

寻找反素数

isAntiPrime 函数可以通过利用其 canReduce() 函数调整为返回最大的反素数

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

输出：

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

由于过程的组合性质，这对于非常大的数字会有点慢。