123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 27 February 2021
- # https://github.com/trizen
- # Upside down Diophantine equation
- # https://projecteuler.net/problem=748
- # The problem asks for integer solutions to:
- # 1/x^2 + 1/y^2 = 13/z^2
- # It is easy to see that:
- # (x^2 + y^2)/13 = v^4, for some integer v.
- # Multiplying both sides by 13, we have:
- # x^2 + y^2 = v^4 * 13
- # By finding integer solutions (x,y) to the above equation, we can then compute `z` as:
- # z = sqrt((x^2 * y^2 * 13)/(x^2 + y^2))
- # z = sqrt((x^2 * y^2) / v^4)
- # We need to iterate over 1 <= v <= sqrt(N/3).
- func sum_of_two_squares_solutions(n) is cached {
- n == 0 && return [[0, 0]]
- 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 [[prod2, 0]]
- 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(-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*s > 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] }
- }
- func S(N) {
- var total = 0
- for v in (1 .. floor(sqrt(N/3))) {
- var solutions = sum_of_two_squares_solutions(13 * v**4) || next
- for x,y in (solutions) {
- y <= N || next
- var z = (x*x * y*y)/(v**4)
- if (z.is_square) {
- z = z.isqrt
- z <= N || next
- if (gcd(x,y,z) == 1) {
- say [x,y,z,v]
- total += [x,y,z].sum
- assert_eq(1/x**2 + 1/y**2, 13/z**2)
- }
- }
- }
- }
- return total
- }
- assert_eq(S(1e2), 124)
- assert_eq(S(1e3), 1470)
- assert_eq(S(1e5), 2340084)
- var total = S(1e16)
- say "Total: #{total} -> #{total % 1e9}"
|