在Julia中用Int64保存整数，
如果需要许多位的整数就用BigInt类型。

从整数值拆分各个数字

编写一个函数，
从一个正整数如12345，
拆分成各个数字的数组[1,2,3,4,5]。

输入:

  • n: 要拆分的值。
  • base: n的进制。
  • ndigits: 输出的位数，不足时在左边添加0，超过时出错。
function todigits(n, base = 10; ndigits=0)
    x = Int[]
    while n > 0
        (n, d) = divrem(n, base)
        push!(x, d)
    end
    x = reverse(x)
    nx = length(x)
    if ndigits > 0
        if nx < ndigits
            x = [zeros(Int, ndigits - nx); x]
        elseif nx > ndigits
            error("ndigits too small!")
        end
    end
    x     
end

todigits(12345, ndigits=8)'
## 0  0  0  1  2  3  4  5

这个问题Base库已经提供了digits函数，结果从个位倒序输出：

digits(12345)'
## 5  4  3  2  1

素因子分解

素数是大于等于2的整数中，
只能被自己和1整除的数。
等于两个大于1的数的乘积的数称为合数。
1既不是素数，也不是合数。
前几个素数为2,3,5,7,11,13,19，
除了2以外都是奇数。

设n为正整数，
将其分解为各个素数因子的乘积，
如

60=2×2×3×5

先写一个不利用保存的素数序列的版本。
输入正整数，输出从小到大素因子组成的数组，元素可重复。

function prime_factors(n::Integer)
    x = Int[]
    d = BigInt(2)
    while n > 1
        # 如果d是因子，则从n中都除掉
        while n > 1 && n % d == 0 
            n = n ÷ d
            push!(x, d)
        end
        
        if n == 1 # 已经没有素因子
            break
        else
            # 考虑下一个可能的因子，所有因子都>d了
            d += (d == 2 ? 1 : 2)
            if n < d^2
                # 如果n还是合数，n=ab, 则a>=d, b>=d, n >= d^2
                # 所以n < d^2时n为素数
                push!(x, n)
                break
            end
            # 这里n还有可能是合数，继续用d尝试能否除尽
        end # if - else
    end # while n > 1
    x
end
prime_factors(60)'
2  2  3  5

设已经有m个素数保存在一维数组primes中。
在素因子分解时可以利用这些已有的素数，
这样在考虑因子时，就不必考虑所有的奇数了。
程序会复杂一些，
但是速度应该能提高。

输入正整数，和素数表的开头若干个组成的数组，
输出从小到大素因子组成的数组，元素可重复。
注意下面的函数与前面的函数名称相同，
但位置参数不同，
所以算是同一函数的不同方法。

function prime_factors(n::Integer, primes=[2])
    x = Int[]
    np = length(primes)
    if np <= 8
        primes = Int[2, 3, 5, 7, 11, 13, 17, 19]
        np = 8
    end

    ip = 1
    d = primes[ip]
    
    # 循环继续的条件，是n还有大于1的因子
    while n > 1

        # 把当前考虑的因子都除去
        while n > 1 && n % d == 0
            n = n ÷ d
            push!(x, d)
        end
        
        if n == 1
            # 所有素因子都已经找到
            break
        else # n > 1
            # 考虑下一个可能的因子
            if ip < np
                ip += 1
                d = primes[ip]
            else
                d += 2
            end

            if d^2 > n 
                # n是合数要求n至少是d*d，即 n >= d^2, 
                # 所以d^2>n时n不是合数，又有n>1, 所以n是素数
                push!(x, n)
                break 
            end 

            # 进入下一轮循环尝试用d作为因子
        end # if n==1
    end # while n>1
    x
end

比较两个版本：

using BenchmarkTools
@btime prime_factors(12345678901234567890)
@btime prime_factors(12345678901234567890, [2])

发现运算速度相近。
分解结果为[2, 3, 3, 5, 101, 3541, 3607, 3803, 27961]。

在

n很大时会很慢，
如n = BigInt(2)^100 + 5。
有可能输入一个很大的素数表可以帮助。

最大公因数和最小公倍数

Julia的gcd(x, y)计算最大公因数，
lcm(x, y)计算最小公因数。
如：

gcd(30, 24)
## 6
lcm(30, 24)
## 120

求素数表

用逐个试验能否除尽的方法

前面进行素因子分解的程序，只要化简一下就可以变成判断素数的程序：
如果找到了任何一个小于n的因子，
就不是素数，
否则是素数。
程序如下：

# 用逐个除数试除的方法判断是否素数。
function isprime(n::Integer)
    @assert n >= 1
    if n==1
        return false
    elseif n==2
        return true
    end 
    
    prime = false 
    d = 2 
    while n % d != 0 && d^2 < n 
        # 只要当前除数不能整除n，就考虑下一个除数
        if d==2
            d = 3
        else 
            d += 2
        end 
    end
    # 退出while循环时，如果不是因为d^2 > n退出，就是被中途
    # 两种退出原因就决定了是素数还是合数
    if n % d == 0
        return false
    else
        return true
    end 
end
show(filter(isprime, [1:100;]))
## [2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
##  31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
##  73, 79, 83, 89, 97]

筛法求素数表

为了获得

N以内的素数表，
可以用如下的方法：
将

1:N罗列出来，
划去1，
然后划去2除本身以外的倍数，
再划去3除本身以外的倍数，
然后划去下一个没有被划去的数的2倍及以上的倍数，
如此重复一直到找不到下一个未被划去的数为止。
最后输出没有被划去的数。

用Julia实现，
输入上限n，
设置一个布尔型数组x[1:n]，
初值为true,
数k被划去，
就设x[k]=false。

md"""
筛法求素数表。

输入上限n，输出n以内的素数的数组。
"""
function primes_sieve(n)
    x = Array{Bool}(undef, n)
    x .= true
    x[1] = false
    d = BigInt(2)
    while d <= n 
        pos = d * 2
        ## 各个倍数
        while pos <= n 
            x[pos] = false 
            pos += d 
        end 
        
        ## 找到下一个未被划去的数
        while true 
            d += 1
            if d > n || x[d]
                break
            end 
        end 
    end 

    [i for i=1:n if x[i]]
end
string(primes_sieve(100))
## [2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
##  31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
##  73, 79, 83, 89, 97]

水仙花数和兰德尔数

如果十进制的

abc10等于其三个数字的立方和，
称为水仙花数。
如果十进制的

n位数等于其数字的

n次方的和，
称为兰德尔数。

下面用Julia求出

n位数的兰德尔数。

function randle(ndigits)
    res = Int[]
    for x = (10^(ndigits-1)):(10^ndigits - 1)
        ds = digits(x)
        if x == sum(ds .^ ndigits)
            push!(res, x)
        end 
    end 

    return res
end 
show(randle(3)) ## [153, 370, 371, 407]
show(randle(4)) ## [1634, 8208, 9474]
show(randle(5)) ## [54748, 92727, 93084]
show(randle(6)) ## [548834]
show(randle(7)) ## [1741725, 4210818, 9800817, 9926315]
show(randle(8)) ## [24678050, 24678051, 88593477]

勾股数

当勾股定理中三个边长都是整数的时候，
称这三个数为勾股数，
如3, 4, 5:

32+42=52

查找斜边在

n以下的勾股数，去掉三个边有大于1的公因子的情况。
对两个直角边长做二重循环查找。

md"""
斜边长在n以下的勾股数。
"""
function pythagoras_numbers(n)
    res = [(0,0,0)]
    for a in 1:(n ÷ 2)
        for b in a:(n-1)
            x = a*a + b*b 
            c = round(Int, sqrt(x))
            if x == c*c && c ≤ n && gcd(a, b, c)==1
                push!(res, (a, b, c))
            end 
        end 
    end 
    return res[2:end]
end 

print(pythagoras_numbers(100))
## [(3, 4, 5), (5, 12, 13), (7, 24, 25), (8, 15, 17), 
##  (9, 40, 41), (11, 60, 61), (12, 35, 37), (13, 84, 85), 
##  (16, 63, 65), (20, 21, 29), (28, 45, 53), (33, 56, 65), 
##  (36, 77, 85), (39, 80, 89), (48, 55, 73)]

角谷猜想

考虑如下的自然数递推序列：

xn+1={xn/2,若xn为偶数,3xn+1,若xn为奇数.

比如，

x0=5,
序列为

5,16,8,4,2,1，
从这里开始就会不停重复

4,2,1。
称这样的序列为“冰雹序列”。
角谷猜想是，
从任何正整数出发这样递推都会回到1，
现在还没有发现任何反例，也没有数学证明。
将从

x0出发到1为止的步数称为此

x0对应的冰雹序列长度。
作图显示冰雹序列长度分布。

我们对

1:n的所有

x0都计算长度，
并用直方图表示长度分布。

计算冰雹序列长度的函数：

function hailstone(n::Int)
    len = 0
    while n ≠ 1
        if iseven(n)
            n = n ÷ 2
        else
            n = 3n + 1
        end
        len += 1
    end

    return len 
end

作直方图测试的程序：

using CairoMakie
CairoMakie.activate!()
ys = [hailstone(n) for n=2:10_000_000];
hist(ys, bins=1000) |> display
knitr::include_graphics("figs/exa-numbers-hail001.png")

北京大学Julia语言入门讲义第22章： Julia编程示例–自然数处理

对前一千万个自然数计算了冰雹序列长度，
用时仅几秒钟。