Engee 文档
Notebook

确定友好数字(友好对)的算法

我们正在考虑实现一种算法,用于搜索高达给定限制的友好数字对。

导言

什么是友好号码?

友好数字是两个自然数,使得一个数字的适当除数的总和等于另一个数字,反之亦然。 例如,一对夫妇 它很友好:

\是它自己的除数的总和

\是它自己的除数的总和

让我们开发一个算法,找到所有这些对达到指定的限制。

主要部分

在使用外部软件包之前,您需要安装它。:

In [ ]:
import Pkg; Pkg.add("Primes")
   Resolving package versions...
    Updating `~/.project/Project.toml`
  [27ebfcd6] + Primes v0.5.7
    Updating `~/.project/Manifest.toml`
  [18e54dd8] + IntegerMathUtils v0.1.3
  [27ebfcd6] + Primes v0.5.7

连接必要的库:

In [ ]:
using Primes, Printf  # Primes для проверки простоты чисел, Printf для форматированного вывода

计算素数乘数贡献的函数

此函数计算素数的幂的总和。 从度 以前 .

In [ ]:
function pcontrib(p::Int64, a::Int64)
    n = one(p)      # Начинаем с p^0 = 1
    pcon = one(p)   # Инициализация суммы значением 1 (p^0)
    for i in 1:a    # Для каждой степени от 1 до a
        n *= p      # Увеличиваем степень p: p^i
        pcon += n   # Добавляем p^i к сумме
    end
    return pcon     # Возвращаем сумму: 1 + p + p^2 + ... + p^a
end
Out[0]:
pcontrib (generic function with 1 method)

函数操作说明 pcontrib:

要计算一个数字的所有除数的总和,必须考虑到如果数字具有规范分解 ,那么它的所有除数的总和将等于产品:

功能 pcontrib 为特定素因子计算一个这样的括号。

计算所有除数之和的函数

使用其素因式分解计算数的所有除数的总和。

In [ ]:
function divisorsum(n::Int64)
    dsum = one(n)           # Начальное значение произведения - 1
    for (p, a) in factor(n) # Для каждого простого множителя p степени a в разложении n
        dsum *= pcontrib(p, a) # Домножаем на сумму (1 + p + p^2 + ... + p^a)
    end
    dsum -= n               # Вычитаем само число, чтобы получить сумму только собственных делителей
end
Out[0]:
divisorsum (generic function with 1 method)

功能说明 divisorsum:

使用算术的基本定理-每个自然数都可以表示为素数幂的乘积(规范分解)。 通过系数的乘积计算所有除数的总和 pcontrib. 然后从这个数量中减去数字本身。 以获得仅其自身除数的总和(小于数字本身)。

寻找友好夫妇的主要功能

主要功能,搜索所有对友好的数字达到极限。

In [ ]:
function amicables(L = 2*10^7)  # Поиск пар до значения L (по умолчанию 20 миллионов)
    acnt = 0                    # Счетчик найденных пар
    
    # Вывод информационного сообщения
    println("Amicable pairs not greater than ", L)
    
    # Перебираем все числа от 2 до L
    for i in 2:L
        !isprime(i) || continue # Пропускаем простые числа (у них сумма собственных делителей всегда 1)
        
        j = divisorsum(i)       # Вычисляем сумму собственных делителей числа i
        
        # Проверяем условие дружественности:
        # - j должен быть меньше i (чтобы не повторяться)
        # - сумма собственных делителей j должна быть равна i
        j < i && divisorsum(j) == i || continue 
        
        # Если пара найдена, увеличиваем счетчик и выводим результат
        acnt += 1
        println(@sprintf("%4d", acnt), " => ", j, ", ", i)
    end
end
Out[0]:
amicables (generic function with 2 methods)

函数中优化的解释 amicables:

循环从 以前 按升序排列。

支票 确保每对只找到一次,即如果一对已经见过面。 然后它将不再输出为 . 素数也被跳过,因为素数的适当除数之和是 ,而 永远无法形成友好的一对。

In [ ]:
amicables()  # Запуск поиска с параметрами по умолчанию (до 20,000,000)
Amicable pairs not greater than 20000000
   1 => 220, 284
   2 => 1184, 1210
   3 => 2620, 2924
   4 => 5020, 5564
   5 => 6232, 6368
   6 => 10744, 10856
   7 => 12285, 14595
   8 => 17296, 18416
   9 => 66928, 66992
  10 => 67095, 71145
  11 => 63020, 76084
  12 => 69615, 87633
  13 => 79750, 88730
  14 => 122368, 123152
  15 => 100485, 124155
  16 => 122265, 139815
  17 => 141664, 153176
  18 => 142310, 168730
  19 => 171856, 176336
  20 => 176272, 180848
  21 => 196724, 202444
  22 => 185368, 203432
  23 => 280540, 365084
  24 => 308620, 389924
  25 => 356408, 399592
  26 => 319550, 430402
  27 => 437456, 455344
  28 => 469028, 486178
  29 => 503056, 514736
  30 => 522405, 525915
  31 => 643336, 652664
  32 => 600392, 669688
  33 => 609928, 686072
  34 => 624184, 691256
  35 => 635624, 712216
  36 => 667964, 783556
  37 => 726104, 796696
  38 => 802725, 863835
  39 => 879712, 901424
  40 => 898216, 980984
  41 => 998104, 1043096
  42 => 1077890, 1099390
  43 => 947835, 1125765
  44 => 1154450, 1189150
  45 => 1185376, 1286744
  46 => 1156870, 1292570
  47 => 1280565, 1340235
  48 => 1175265, 1438983
  49 => 1392368, 1464592
  50 => 1328470, 1483850
  51 => 1358595, 1486845
  52 => 1511930, 1598470
  53 => 1466150, 1747930
  54 => 1468324, 1749212
  55 => 1798875, 1870245
  56 => 1669910, 2062570
  57 => 2082464, 2090656
  58 => 2236570, 2429030
  59 => 2723792, 2874064
  60 => 2739704, 2928136
  61 => 2652728, 2941672
  62 => 2802416, 2947216
  63 => 2728726, 3077354
  64 => 2803580, 3716164
  65 => 3276856, 3721544
  66 => 3606850, 3892670
  67 => 3805264, 4006736
  68 => 3786904, 4300136
  69 => 4238984, 4314616
  70 => 4259750, 4445050
  71 => 4246130, 4488910
  72 => 4482765, 5120595
  73 => 4604776, 5162744
  74 => 5459176, 5495264
  75 => 5123090, 5504110
  76 => 5357625, 5684679
  77 => 5232010, 5799542
  78 => 5385310, 5812130
  79 => 5147032, 5843048
  80 => 5730615, 6088905
  81 => 4532710, 6135962
  82 => 5726072, 6369928
  83 => 6329416, 6371384
  84 => 6377175, 6680025
  85 => 6993610, 7158710
  86 => 6955216, 7418864
  87 => 7275532, 7471508
  88 => 5864660, 7489324
  89 => 7489112, 7674088
  90 => 7677248, 7684672
  91 => 7800544, 7916696
  92 => 7850512, 8052488
  93 => 7288930, 8221598
  94 => 8262136, 8369864
  95 => 7577350, 8493050
  96 => 9363584, 9437056
  97 => 9071685, 9498555
  98 => 9199496, 9592504
  99 => 8619765, 9627915
 100 => 9339704, 9892936
 101 => 9660950, 10025290
 102 => 8826070, 10043690
 103 => 10254970, 10273670
 104 => 8666860, 10638356
 105 => 9206925, 10791795
 106 => 10572550, 10854650
 107 => 8754130, 10893230
 108 => 10533296, 10949704
 109 => 9491625, 10950615
 110 => 9478910, 11049730
 111 => 10596368, 11199112
 112 => 9773505, 11791935
 113 => 11498355, 12024045
 114 => 10992735, 12070305
 115 => 11252648, 12101272
 116 => 11545616, 12247504
 117 => 11693290, 12361622
 118 => 12397552, 13136528
 119 => 11173460, 13212076
 120 => 11905504, 13337336
 121 => 13921528, 13985672
 122 => 10634085, 14084763
 123 => 12707704, 14236136
 124 => 13813150, 14310050
 125 => 14311688, 14718712
 126 => 15002464, 15334304
 127 => 13671735, 15877065
 128 => 14443730, 15882670
 129 => 16137628, 16150628
 130 => 15363832, 16517768
 131 => 14654150, 16817050
 132 => 15938055, 17308665
 133 => 17257695, 17578785
 134 => 17908064, 18017056
 135 => 14426230, 18087818
 136 => 18056312, 18166888
 137 => 17041010, 19150222
 138 => 18655744, 19154336
 139 => 16871582, 19325698
 140 => 17844255, 19895265
 141 => 17754165, 19985355

程序的结果将包含一个升序排列的友好数字对列表,其中第一个数字小于第二个。

输出格式: <порядковый_номер> => <меньшее_число>, <большее_число>

结论

我们回顾了友好数字搜索算法在Julia编程语言中的实现。

该算法利用将数字分解为素数因子的数学属性来有效地计算它们自己的除数之和。 效率是通过使用数论而不是直接通过所有除数来实现的。 该程序查找到指定限制的所有友好对,并用编号输出它们。 这种方法在数论和娱乐数学中有实际应用。 有趣的是,古希腊人知道友好的数字,最小的一对(220,284)是毕达哥拉斯发现的。

该示例是使用[Rosetta代码]的材料开发的(https://rosettacode.org/wiki/Amicable_pairs