-
Notifications
You must be signed in to change notification settings - Fork 6
/
009 Special Pythagorean triplet -- v2.sf
94 lines (74 loc) · 2.19 KB
/
009 Special Pythagorean triplet -- v2.sf
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
#!/usr/bin/ruby
# Daniel "Trizen" Șuteu
# Date: 09 May 2018
# https://github.com/trizen
# https://projecteuler.net/problem=9
# Runtime: 0.387s
func sum_of_two_squares_solutions(n) is cached {
var prod1 = 1
var prod2 = 1
var prime_powers = []
for p,e in (n.factor_exp) {
if (p % 4 == 3) { # p = 3 (mod 4)
e.is_even || return [] # power must be even
prod2 *= p**(e >> 1)
}
elsif (p == 2) { # p = 2
if (e.is_even) { # power is even
prod2 *= p**(e >> 1)
}
else { # power is odd
prod1 *= p
prod2 *= p**((e - 1) >> 1)
prime_powers.append([p, 1])
}
}
else { # p = 1 (mod 4)
prod1 *= p**e
prime_powers.append([p, e])
}
}
prod1 == 1 && return []
prod1 == 2 && return [[prod2, prod2]]
# All the solutions to the congruence: x^2 = -1 (mod prod1)
var square_roots = gather {
gather {
for p,e in (prime_powers) {
var pp = p**e
var r = sqrtmod(pp-1, pp)
take([[r, pp], [pp - r, pp]])
}
}.cartesian { |*a|
take(Math.chinese(a...))
}
}
var solutions = []
for r in (square_roots) {
var s = r
var q = prod1
while (s.sqr > prod1) {
(s, q) = (q % s, s)
}
solutions.append([prod2 * s, prod2 * (q % s)])
}
for p,e in (prime_powers) {
for (var i = e%2; i < e; i += 2) {
var sq = p**((e - i) >> 1)
var pp = p**(e - i)
solutions += (
__FUNC__(prod1 / pp).map { |pair|
pair.map {|r| sq * prod2 * r }
}
)
}
}
solutions.map {|pair| pair.sort } \
.uniq_by {|pair| pair[0] } \
.sort_by {|pair| pair[0] }
}
for n in (1..Inf) {
var a = sum_of_two_squares_solutions(n**2) || next
var s = a.first { .sum + n == 1000 } || next
say (s.prod * n)
break
}