ProjectEuler/P448.py at master · mcleu/ProjectEuler · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# -*- coding: utf-8 -*-
#==============================================================================
# The function lcm(a,b) denotes the least common multiple of a and b.
# Let A(n) be the average of the values of lcm(n,i) for 1<=i<=n.
# E.g: A(2)=(2+2)/2=2 and A(10)=(10+10+30+20+10+30+70+40+90+10)/10=32.
#
# Let S(n)=SUM A(k) for 1<=k<=n.
# S(100) 122726
# Find S(99999999019) mod 999999017.
#==============================================================================


from time import time
import numpy as np

def lcm1(a, b):
    A = max(a,b)
    B = min(a,b)
    LCM = B
    while True:
        #print LCM, A, B, LCM%A, LCM%B
        if (LCM % A == 0) and (LCM % B == 0):
            return LCM
        LCM += B

def gcd(a, b):
    """Return greatest common divisor using Euclid's Algorithm."""
    while b:
        a, b = b, a % b
    return a

def npgcd(a, b):
    a, b = np.broadcast_arrays(a, b)
    a = a.copy()
    b = b.copy()
    pos = np.nonzero(b)[0]
    while len(pos) > 0:
        b2 = b[pos]
        a[pos], b[pos] = b2, a[pos] % b2
        pos = pos[b[pos]!=0]
    return a

def lcm(a, b):
    """Return lowest common multiple."""
    return a * b // gcd(a, b)


def A(n):
    lcmsum = 0
    i = n
    while i > 0:
        lcmres = lcm(n, i)
        if i%100000==0:
            print i, lcmsum, (999999017*n)
        lcmsum = (lcmsum + lcmres) % (999999017*n)
        i -= 1
    return lcmsum / n


def S(nn):
    ssum = 0
    kk = nn
    while kk > 0:
        #if k%500==0:
        #    print k, ssum
        ssum = (ssum + A(kk) ) % 999999017
        #print 'S', k, 'Done!', ssum
        kk -= 1
    return ssum


def Sfast(nn):
    newsum = 0
    jj = nn
    while (jj > 0):
        #if jj%10==0:
        #    print '10 down'
        #print jj
        ii = jj
        while (ii > 0):
            newsum = (newsum + ii / gcd(ii,jj) ) % 999999017
            #print ii,jj #,gcd(ii,jj), ii / gcd(ii,jj)
            ii -= 1
        jj -= 1
    return newsum

def S2(nn):
    sum2 = 0
    jj = nn
    while (jj > 0):
        minisum = 0
        ii = nn
        while ( ii >= jj ):
            #print ii,jj
            minisum += 1.0 / gcd(ii,jj)
            ii -= 1
        sum2 += minisum * jj
        jj -= 1
    return int(sum2)

def Svec(nn):
    sum2 = 0
    jj = nn
    while (jj > 0):
        ii = np.array(range(jj,nn+1))
        ms = 1.0 / npgcd(ii,[jj])
        sum2 = (sum2 + np.sum(ms) * jj)  % 999999017
        jj -= 1
    return int(sum2)


NN = 1200

t1=time()
print Svec(NN)
print time() - t1

t2=time()
print Sfast(NN)
print time() - t2

#Snew(100)