123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384 |
- #!/usr/bin/ruby
- # Pythagorean Triple Occurrence
- # https://projecteuler.net/problem=827
- # For a given number n, the number of Pythagorean triples in which the number n occurs, is equivalent with:
- # (number of solutions to x^2 - y^2 = n^2, with x,y > 0) + (number of solutions to x^2 + y^2 = n^2, with x,y > 0).
- # In Sidef, this is:
- # (diff_of_squares(n**2).len - 1) + (sum_of_squares(n**2).len - 1)
- func count_of_pythagorean_triples(n) {
- # Translation of the PARI/GP program by Michel Marcus, Mar 07 2016.
- # https://oeis.org/A046081
- var n0_sigma0 = 1
- var n1_sigma0 = 1
- for p,e in (n.factor_exp) {
- if (p == 2) {
- n0_sigma0 *= (2*e - 1)
- }
- else {
- n0_sigma0 *= (2*e + 1)
- if (p % 4 == 1) {
- n1_sigma0 *= (2*e + 1)
- }
- }
- }
- (n0_sigma0 + n1_sigma0)/2 - 1
- }
- func smallest_number_with_n_divisors(threshold, upperbound = false, least_solution = Inf, k = 1, max_a = Inf, solutions = 1, n = 1) {
- if (solutions == threshold) {
- return n
- }
- if (solutions > threshold) {
- return least_solution
- }
- var p = k.prime
- for a in (1 .. max_a) {
- n *= p
- break if (n > least_solution)
- var new_solutions = count_of_pythagorean_triples(n)
- if (new_solutions <= threshold) {
- least_solution = __FUNC__(threshold, upperbound, least_solution, k+1, (upperbound ? a : Inf), new_solutions, n)
- }
- }
- return least_solution
- }
- func p827(n) {
- var sum = 0
- for k in (1..n) {
- say "Computing upper-bound for k = #{k}"
- var upperbound = smallest_number_with_n_divisors(10**k, true)
- say "Upper-bound: #{upperbound}"
- var value = smallest_number_with_n_divisors(10**k, false, upperbound)
- say "Exact value: #{value}\n"
- sum += value
- }
- return sum
- }
- say count_of_pythagorean_triples(48)
- say count_of_pythagorean_triples(8064000)
- assert_eq(
- 1000.of(count_of_pythagorean_triples).slice(1),
- 1000.of { .sqr.sum_of_squares.len.dec + .sqr.diff_of_squares.len.dec }.slice(1)
- )
- say p827(18)
|