[欧拉计划]Julia尝试 1:10

nesteiner · 2020 年3 月 29 日 10:43

我打算作死用Julia去写欧拉计划了，能更新多少是多少

如果有更好的解请在评论里留言

辅助模块

习题

solution 1

答案:234168

Problem 1 | Project Euler | 欧拉计划

# Problem 1
# 3的倍数和5的倍数
# 如果我们列出10以内所有3或5的倍数，我们将得到3、5、6和9，这些数的和是23。

# 求1000以内所有3或5的倍数的和。
function solution_1()
    fn = x->x % 3==0 || x % 5 ==0
    reduce((r,x)->r+x,
           filter(fn,1:1000);init=0)
end

solution2

ans: 4613732

Problem 2 | Project Euler | 欧拉计划

# Problem 2
# 偶斐波那契数
# 斐波那契数列中的每一项都是前两项的和。由1和2开始生成的斐波那契数列前10项为：

# 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …
# 考虑该斐波那契数列中不超过四百万的项，求其中为偶数的项之和。
function solution_2()
    reduce(+,
           Iterators.filter(x->x % 2==0,
                            Iterators.filter(x->x<4000000,
                                             Iterators.take(fib,40))))
end

solution3

ans: 6857

Problem 3 | Project Euler | 欧拉计划

divsible(a,b) = a % b == 0
function factor(n,__primes=collect(Iterators.take(primes,1000)),res=[])
    div = first(__primes)
    if divsible(n,div)
        if n == div
            return push!(res,n)
        else
            return factor(Int(n/div),__primes,push!(res,div))
        end
    else
        return factor(n,Iterators.drop(__primes,1),res)
    end
end

function solution_3()
    last(factor(600851475143))
end

solution4

ans: 580085

Problem 4 | Project Euler | 欧拉计划

Clojure 实现

  (let [isPalindrome? (fn [n]
                        (let [s (seq (str n))]
                          (= s (reverse s))))]
    (first (for [num1 (range 999 99 -1)
                 num2 (range 999 99 -1)
                 :let [n (* num1 num2)]
                 :when (isPalindrome? n)]
             n))))

转成Julia

function solution_4()
    first(Iterators.filter(isPalindrome,
                           Base.Generator(x->reduce(*,x;init=1),
                                          Iterators.product(999:-1:100,999:-1:100))))
end

solution5

ans: 232792560

Problem 5 | Project Euler | 欧拉计划

function gcd(a,b)
    temp = a % b
    if temp == 0
        b
    else
        gcd(b,temp)
    end
end


function lcm(a,b)
    g = gcd(a,b)
    d = a / g
    Int(b * d)
end

function solution_5()
    reduce((r,x)->lcm(r,x),1:20;init=1)
end

solution 6

ans: 25164150

Problem 6 | Project Euler | 欧拉计划

function solution_6()
    nums = 1:100
    a = reduce(+,nums;init=0)^2
    b = reduce((r,x)->r+x^2,nums;init=0)
    a-b
end

solution7

ans: 104743

Problem 7 | Project Euler | 欧拉计划

function solution_7()
    last(collect(Iterators.take(primes,10001)))
end

solution9

ans: 31875000
纯数学方法

solution10

ans: 我已经运行十几分钟了，就差程序错误，还是得不到答案
142913828922

Problem 10 | Project Euler | 欧拉计划

用了一次break，虽然不想用这个东西

function solution_10()
    sum = 0
    for i in primes
        if i < 2000000
            sum += i
        else
            break
        end
    end
    sum
end

Sukanka · 2020 年4 月 11 日 06:12

习题

solution 1

ans: 234168

function solution1(n=1000)
sum(x for x in 3:n if x%3==0 || x%5==0)
end

solution 2

ans: 4613732

function solution2(n=4000000)
    φ=(1-√5)/2
    fib(n,a1=1,a2=2)=round(Int,((5+3*√5)a1-2*√5a2)/10*φ^n + ((5-3*√5)a1+2*√5a2)/10 *(1-φ)^n)
    s=0
    i=1
    f=fib(i)
    while f<n
        if f%2==0
            s+=f
        end
        i+=1
        f=fib(i)
    end
    s
end

solution 3

ans: 6857

function solution3(n=600851475143) #Pollard Rho因数分解
    k=2
    while n!=k
        if n %k==0
            n=round(Int,n/k)
        else
            k+= k==2 ? 1 : 2
        end
    end
    return k
end

solution 4

ans: 580085

function solution4(n=3) 
    for i in 10^n-1:-1:10^(n-1)+1
        for j in i:-1:10^(n-1)
            product=i*j
            if "$product"==reverse("$product")
                return product
            end
        end
    end
end

solution 5

ans: 232792560

function solution5(n=20)
    function gcd(a::T1,b::T2) where {T1,T2 <:Integer}
        if a> b
            return gcd(b,a)
        else
            c=b%a
            if c==0
                return a
                else return gcd(a,c)
            end
        end
    end
    function gcd(a::T) where T<: AbstractArray{<:Integer,1}
        m=a[1]
        for i in 1:length(a)
            m=gcd(m,a[i])
        end
        m
    end          
    lcm(a::T1,b::T2) where {T1,T2 <:Integer} =round(Int,a*b/gcd(a,b))
    function lcm(a::T) where T<: AbstractArray{<:Integer,1}
        m=a[1]
        for i in 1:length(a)
            m=lcm(m,a[i])
        end
        m
    end
    lcm(1:n) 
end

solution 6

ans: 25164150

solution6(n=100)=sum(1:n)^2-sum((1:n).^2)

solution 7

ans: 104743

function solution7(n=10001)
    primes=zeros(Int,n)
    primes[1:2].=[2,3]
    k=3
    counter=2
    while counter < n
        k+=2
        is_prime=true
        for i in primes[2:counter]
            if k%i ==0
                is_prime=false
                break
            end
        end
        if is_prime
            counter+=1
            primes[counter]=k
        end
    end
    k
end

solution 8

ans: 23514624000

function solution8(cont=13)
    str="""
    73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450
    """
    str=strip(replace(str,"\n"=>""))
    max_prod=0
    for i in 1:length(str)-cont
        s=prod(parse(Int,num) for num in str[i:i+cont-1] )
        if s> max_prod
            max_prod=s
        end
    end
    max_prod
end

solution 9

ans: 31875000

function solution9(n=1000)
    for i in 1:floor(Int,n/3)
        for j in i:floor(Int,(n-i)/2)
            if  i^2+j^2==(n-j-i)^2
                return i*j*(n-i-j)
            end
        end
    end
end

solution 10

ans: 142913828922

function solution10(n=2000000)
    function sieve(n)
        is_prime=trues(n)
            for i in 2:ceil(Int,√n)
                if is_prime[i]
                    for j in i^2:i:n
                        if j %i==0
                            is_prime[j]=false
                        end
                    end
                end
            end
        return (x for x in 2:n if is_prime[x])
    end
    return sum(sieve(n))
end

solution 10 的性能情况

julia> @benchmark solution10()
BenchmarkTools.Trial: 
  memory estimate:  244.38 KiB
  allocs estimate:  7
  --------------
  minimum time:     24.380 ms (0.00% GC)
  median time:      24.453 ms (0.00% GC)
  mean time:        24.489 ms (0.04% GC)
  maximum time:     26.572 ms (7.93% GC)
  --------------
  samples:          205
  evals/sample:     1

wushaohua5 · 2020 年4 月 11 日 06:18

刚开始追楼主主题，solution_2有点自己的想法。个人觉得问题应该是数列中不大于40000000前偶数和，楼主求的是前40000000项偶数和；解题方法用数组存储数列然后过滤求和，比较合适。试了下，运行时间不到一秒。顺便问一下，怎么解除Julia的递归深度限制，用线性递归求斐波那契数列报错，程序没问题

nesteiner · 2020 年4 月 11 日 06:32

不知道诶，函数式编程里有一个概念叫惰性，在命令式编程里的体现是迭代器，叫iterator，你可以去查查Julia的iterate文档
递归这东西可以体现惰性的好处，比如我需要一个无限的ones数组，那么在Clojure中是这样的

(def ones
    (lazy-seq
        (cons 1 ones))

如果我需要５个1，那么我只用(take 5 ones)得到(1 1 1 1 1)

待会我把无线fib数列的说明补全，你可以参考

nesteiner · 2020 年4 月 11 日 08:01

补好了

nesteiner · 2020 年4 月 11 日 08:06

第二题算出来了

nesteiner · 2020 年4 月 27 日 11:27

来这里吧，如果有梯子的话，加好友的key是1635154_Yzd5dbO3hjRiRbLl1OfOj4uEs0AaHum0

AquaIndigo · 2020 年7 月 18 日 03:07

第十题

julia> using Primes

julia> @btime sum(primes(20_00_000))
  5.040 ms (5 allocations: 1.65 MiB)
142913828922