算法-数论

1、最大公约数

def gcd(a,b):
    if b == 0:
        return a
    return gcd(b, a%b) # a和b的最大公约数等于b与a mod b 的最大公约数

def gcd(a,b):
    while b != 0:
        cur = a
        a = b
        b = cur%b
        pass
    return a

欧几里得算法

a可以表示成a = kb + r（a，b，k，r皆为正整数，且r不为0）

假设d是a,b的一个公约数，记作d|a,d|b，即a和b都可以被d整除。

而r = a - kb，两边同时除以d，r/d=a/d-kb/d，由等式右边可知m=r/d为整数，因此d|r

因此d也是b,a mod b的公约数。

因(a,b)和(b,a mod b)的公约数相等，则其最大公约数也相等，得证。

例题：

1、蓝桥杯2019年第十届省赛真题-等差数列 - C语言网 (dotcpp.com)

import math
n=int(input())
arr=list(map(int,input().split()))


if arr[0] > arr[1]:
    max_ = arr[0]
    min_ = arr[1]
else:
    max_ = arr[1]
    min_ = arr[0]
sub = max_ -min_
for i in range(2, n):
    sub = math.gcd(sub, abs(arr[i]-arr[i-1]))
    min_ = min(min_, arr[i])
    max_ = max(max_, arr[i])
if max_ == min_ : print(n)
else:
    print((max_ - min_) // sub + 1)

2、1223. 最大比例 - AcWing题库

辗转相减法

n = int(input())
arr = list(map(int, input().split()))
arr.sort()

def gcd(a,b):
    if b == 0:
        return a
    return gcd(b, a%b)
top, bot = 0,0
def gcd_sub(a, b):
    if a < b:
        a,b = b,a
    if b==1:
        return a
    return gcd_sub(b, a//b)
i = 1
while i < n:
    if arr[i] != arr[0]:
        gcd_ = gcd(arr[i], arr[0])
        top = arr[i]//gcd_
        bot = arr[0]//gcd_
        break
    i += 1
if i == n:
    print(1)
else:
    for i in range(i, n):
        gcd_ = gcd(arr[i], arr[0])
        a = arr[i]//gcd_
        b = arr[0]//gcd_
        top = gcd_sub(a, top)
        bot = gcd_sub(b, bot)
        
print(f'{top}/{bot}')
        

1.1 扩展欧几里得定律

def ext_euclid(a, b):     
    if b == 0:         
        return 1, 0, a     
    else:         
        x, y, q = ext_euclid(b, a % b) 
        # q = gcd(a, b) = gcd(b, a%b)         
        x, y = y, (x - (a // b) * y)         
        return x, y, q

扩展欧几里得算法求系数x，y_哔哩哔哩_bilibili

例题：1299. 五指山 - AcWing题库

def olai(a, b):
    if b == 0:
        return 1,0,a
    x, y, gcd = olai(b, a%b)
    x, y = y, x-(a//b)*y
    return x, y, gcd
    
    
def funt(n, d, x, y):
    x1, y1, gcd = olai(n, d)
    if (y-x)%gcd != 0:
        print('Impossible')
    else:
        d = y1*(y-x)//gcd
        n //= gcd # ax+by = n;   x = x+k*(b/gcd(a,b)); y = y+k*(a/gcd(a,b))
        print((d%n+n)%n)
    
for _ in range(int(input())):
    n, d, x, y = map(int, input().split())
    funt(n, d, x, y)

2、最小公倍数

def funt(a,b):
    return a*b//gcd(a,b)

最小公倍数 * 最大公约数 == a*b

3、质数筛

def getPrimes(n):
    is_primes = [True]*(n+1)  # 初始化一个布尔数组，表示从2到n的所有数都是质数
    primes = [] # 存储质数
    for i in range(2, int(n**(1/2))+1):
        if is_primes[i]:
            primes.append(i)
            # 将当前质数的所有倍数标记为非质数
            for j in range(i*i, n+1, i): 
                is_primes[j] = False
    # 后面的质数的倍数一定会超
    for i in range(int(n**(1/2))+1, n+1):
        if is_primes[i]:
            primes.append(i)
    return primes

4、区间质数筛

import math

def simple_sieve(limit):
    primes = []
    sieve = [True] * (limit + 1)
    # 0和1不是质数，所以标记为False
    sieve[0] = sieve[1] = False

    # 从2到根号limit遍历数字
    for num in range(2, int(math.sqrt(limit)) + 1):
        if sieve[num]:
            primes.append(num)
            for multiple in range(num * num, limit + 1, num):
                sieve[multiple] = False

    for num in range(int(math.sqrt(limit)) + 1, limit + 1):
        if sieve[num]:
            primes.append(num)

    return primes

def segmented_sieve(start, end):
    # 计算质数的上限
    limit = int(math.sqrt(end)) + 1
    primes = simple_sieve(limit)
    primes_count = len(primes)
    # 创建一个布尔数组，用于标记区间内的数字是否为质数，初始化为True
    sieve = [True] * (end - start + 1)

    # 对于每一个小于等于根号end的质数
    for i in range(primes_count):
        # 计算刚好大于等于start的primes[i]的数
        base = max(primes[i]*primes[i], ((start + primes[i] - 1) // primes[i]) * primes[i])

        # 将当前质数的倍数标记为非质数
        for j in range(base, end + 1, primes[i]):
            sieve[j - start] = False

    # 生成区间内的质数列表
    segmented_primes = [start + i for i in range(end - start + 1) if sieve[i]]
    return segmented_primes

start = 10
end = 50
segmented_sieve(start, end)

5、欧拉函数

def funt(n):
    count = n
    p = 2
    while p*p <= n:
        # 找到质因子
        if n % p == 0:
            while n % p == 0:
                n //= p
            count -= count//p  # n*(1-p)
        p += 1
    if n > 1:
        count -= count//n
    return count
funt(20)

欧拉函数公式证明_哔哩哔哩_bilibili

6、计算质因子个数

def funt(n):
    count = 0
    p = 2
    while p * p <= n:
        if n % p == 0:
            count += 1
            while n % p == 0:
                n //= p
        p += 1
    if n > 1:
        count += 1
    return count
funt(12) 

7、计算所有约数个数

约数，又称因数。整数a除以整数b(b≠0) 除得的商正好是整数而没有余数，就说a能被b整除，或b能整除a。a称为b的倍数，b称为a的约数。

def funt(n):
    count = 1
    p = 2
    while p*p <= n:
        cnt = 0
        while n%p == 0:
            cnt += 1
            n //= p
        count *= (cnt+1)
    if n>1:
        count *= 2
    return count
funt(12)

如12：12可以写成2 ** 2 *3，所以 对于2的选择有3种，幂为0、1、2；3的选择有两种，幂为0、1.

相乘即为约数和

8、计算所有约数和

def funt(n):
    cnt = 0
    p = 1
    while p*p <= n:
        # 一次算两个约数
        if n%p == 0:
            cnt += p
            # 平方就只用算一次
            if p*p != n:
                cnt += n//p
        p += 1
    return cnt
funt(12)

例题

1295. X的因子链 - AcWing题库

# 方法一 动规  
def funt(n):
    primes = [n]
    p = 2
    while p*p <= n:
        if n%p == 0:
            primes.append(p)
            if p*p != n:
                primes.append(n//p)
        p += 1
    if n>1 and n != primes[0]:
        primes.append(n)
    primes.sort()
    dp = [[1]*len(primes) for _ in range(2)]

    for i in range(1, len(primes)):
        max_ = cnt = 0
        for j in range(i):
            if primes[i] % primes[j] == 0:
                if dp[0][j] > max_:
                    cnt = dp[1][j]
                    max_ = dp[0][j]
                elif dp[0][j] == max_:
                    cnt += dp[1][j]
        dp[0][i],dp[1][i] = max_+1, cnt if cnt else 1

    print(dp[0][-1], dp[1][-1])
while True:
    try:
        n = int(input())
        funt(n)
    except EOFError:
        break

        
# 数学！！！ https://www.acwing.com/solution/content/97535/ 具体请看题解，我实在不知道怎么表达o(´^｀)o
def A(a, b):
    cnt = 1
    while b > 0:
        cnt *= a
        a -= 1
        b -= 1
    return cnt

def funt(n):
    cnt = []
    p = 2
    while p*p <= n:
        if n%p == 0:
            cnt_ = 0
            while n%p == 0:
                cnt_ += 1
                n //= p
            cnt.append(cnt_)
        p += 1
    if n>1:
        cnt.append(1)
    sum_ = sum(cnt)

    a = 1
    for i in cnt:
        if i > 1:
            a *= A(i,i)
    times = A(sum_, sum_)//a
    print(sum_, times)
while True:
    try:
        n = int(input())
        funt(n)
    except EOFError:
        break

9、快速幂

# n为负数则最终结果返回 1/pow(x,-n)
def pow(x, n):
    if n < 2:
        return x**n
    ret = pow(x, n//2)
    if n%2:
        return ret*ret*x
    return ret*ret
pow(2, 10)

1:
a *= A(i,i)
times = A(sum_, sum_)//a
print(sum_, times)
while True:
try:
n = int(input())
funt(n)
except EOFError:
break

## 9、快速幂

```Python
# n为负数则最终结果返回 1/pow(x,-n)
def pow(x, n):
    if n < 2:
        return x**n
    ret = pow(x, n//2)
    if n%2:
        return ret*ret*x
    return ret*ret
pow(2, 10)