12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152 |
- #!/usr/bin/ruby
- # Represent an odd prime number p as x^2 + d*y^2 using the Cornacchia–Smith method.
- # Given an odd prime p and a positive integer d not divisible by p, this algorithm
- # either reports that p = x^2 + d*y^2 has no integral solution, or returns a solution.
- # Implementation of "Algorithm 2.3.12" from the book "Prime Numbers - A Computational Perspective".
- # Algorithm extended to any composite number n, returning multiple positive solutions. (it may not return ALL solutions)
- func quadratic_form_representations(d, n) {
- if ((kronecker(-d, n) == -1) && n.is_prime) {
- return []
- }
- var solutions = []
- for x in (sqrtmod_all(-d, n)) {
- if (2*x < n) {
- x = (n - x)
- }
- var (a, b, c) = (n, x, n.isqrt)
- while (b > c) {
- (a,b) = (b, a % b)
- }
- var t = (n - b**2)
- d.divides(t) || next
- idiv(t,d).is_square || next
- solutions << [b, isqrt(idiv(t,d))]
- }
- solutions.uniq!
- solutions.map_2d {|x,y| x**2 + d*y**2 }.each {|v| assert_eq(v, n) }
- return solutions
- }
- say quadratic_form_representations(12, 97) #=> [[7, 2]]
- say quadratic_form_representations(14, 3*5*7) #=> [[7, 2]]
- say quadratic_form_representations(18, 43*97) #=> [[61, 5], [11, 15]]
- say quadratic_form_representations( 5, 3*5*7) #=> [[10, 1], [5, 4]]
- say quadratic_form_representations(1234, 1810585219) #=> [[41087, 315]]
- say quadratic_form_representations(14, 2**127 - 1) #=> [[10229855171020855649, 2162856025860368397]]
- say quadratic_form_representations( 1, 99025) #=> [[256, 183], [297, 104], [312, 41], [311, 48]]
|