123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173 |
- #!/usr/bin/ruby
- # Simple implementation of quadratic integers.
- # See also:
- # https://en.wikipedia.org/wiki/Quadratic_integer
- class QuadraticInteger(a, b, w = 2) { # represents: a + b*sqrt(w)
- method to_s {
- "QuadraticInteger(#{a}, #{b}, #{w})"
- }
- method ==(QuadraticInteger c) {
- (a == c.a) && (b == c.b) && (w == c.w)
- }
- method conjugate {
- QuadraticInteger(a, -b, w)
- }
- method norm {
- a*a - w*b*b
- }
- method add (Number c) {
- QuadraticInteger(a+c, b, w)
- }
- method add (QuadraticInteger z) {
- var (c,d) = (z.a, z.b)
- QuadraticInteger(a+c, b+d, w)
- }
- __CLASS__.alias_method(:add, '+')
- method sub (Number c) {
- QuadraticInteger(a-c, b, w)
- }
- method sub (QuadraticInteger z) {
- var (c,d) = (z.a, z.b)
- QuadraticInteger(a-c, b-d, w)
- }
- __CLASS__.alias_method(:sub, '-')
- method mul (Number c) {
- QuadraticInteger(a*c, b*c, w)
- }
- method mul (QuadraticInteger z) {
- var (c,d) = (z.a, z.b)
- QuadraticInteger(a*c + b*d*w, a*d + b*c, w)
- }
- __CLASS__.alias_method(:mul, '*')
- method mod (Number m) {
- QuadraticInteger(a % m, b % m, w)
- }
- __CLASS__.alias_method(:mod, '%')
- method pow(Number n) {
- var x = self
- var c = QuadraticInteger(1, 0, w)
- for bit in (n.digits(2)) {
- c *= x if bit
- x *= x
- }
- return c
- }
- __CLASS__.alias_method(:pow, '**')
- method powmod(Number n, Number m) {
- var x = self
- var c = QuadraticInteger(1, 0, w)
- for bit in (n.digits(2)) {
- (c *= x) %= m if bit #=
- (x *= x) %= m #=
- }
- return c
- }
- }
- #---------------------------------------------------------------------
- # Determine if a given number is probably a prime number.
- func is_quadratic_pseudoprime (n, r=2) {
- return false if (n <= 1)
- return true if (n <= 3)
- return true if (r <= 0)
- var x = QuadraticInteger(r, 1, r+2).powmod(n, n)
- x.a == r || return false
- var y = QuadraticInteger(r, -1, r+2).powmod(n, n)
- y.a == r || return false
- (x.b + y.b == n) && __FUNC__(n, r-1)
- }
- say is_quadratic_pseudoprime(43) #=> true
- say is_quadratic_pseudoprime(97) #=> true
- with (QuadraticInteger(1, 1, 2)) {|q|
- say 15.of { q.pow(_).a } #=> A001333
- say 15.of { q.pow(_).b } #=> A000129
- }
- with (QuadraticInteger(1, 1, 3)) {|q|
- say 15.of { q.pow(_).a } #=> A026150
- say 15.of { q.pow(_).b } #=> A002605
- }
- var n = (274177-1)
- var m = (2**64 + 1)
- with (QuadraticInteger(3, 4, 2)) {|q|
- var r = q.powmod(n, m)
- say gcd(r.a-1, m) #=> 2741177
- say gcd(r.b, m) #=> 2741177
- }
- #---------------------------------------------------------------------
- func is_gaussian_quadratic_pseudoprime(n) {
- return false if (n <= 1)
- return true if (n <= 3)
- static x = QuadraticInteger(1, -1, -2)
- given (n%8) {
- when ([1,3]) {
- var t = x.powmod(n-1, n)
- (t.a==1 && t.b==0)
- }
- when ([5, 7]) {
- var t = x.powmod(n+1, n)
- (t.a==3 && t.b==0)
- }
- else {
- false
- }
- }
- }
- assert([88561,107185,162401,221761,226801,334153,410041,665281,825265,1569457,1615681,2727649].all(is_gaussian_quadratic_pseudoprime))
- assert([80375707,154287451,267559627,326266051,478614067,573183451,643767931,2433943891,4297753027].all(is_gaussian_quadratic_pseudoprime))
- for n in (1..1e2) {
- if (is_gaussian_quadratic_pseudoprime(n)) {
- if (n.is_composite) {
- say "Counter-example: #{n}"
- }
- }
- elsif (n.is_prime) {
- say "Missed prime: #{n}"
- }
- }
|