def square_roots_default(n)
      return [0] if n == 0
      return [] if jacobi(n) == -1   # Step 1

      # Step 2, except we don't want to use randomness.
      b = (1...prime).find { |i| jacobi(i) == -1 }

      s, t = self.class.factor_out_twos(prime - 1)  # Step 3
      n_inv = inverse(n)  # Step 4

      # Step 5
      c = power b, t
      r = power n, (t + 1) / 2

      # Step 6
      (1...s).each do |i|
        d = power r * r * n_inv, 1 << (s - i - 1)
        r = mod r * c if d == prime - 1
        c = square c
      end

      square_roots_given_candidate n, r
    end
def square_roots_default(n)
      return [0] if n == 0
      return [] if jacobi(n) == -1   # Step 1

      # Step 2, except we don't want to use randomness.
      b = (1...prime).find { |i| jacobi(i) == -1 }

      s, t = self.class.factor_out_twos(prime - 1)  # Step 3
      n_inv = inverse(n)  # Step 4

      # Step 5
      c = power b, t
      r = power n, (t + 1) / 2

      # Step 6
      (1...s).each do |i|
        d = power r * r * n_inv, 1 << (s - i - 1)
        r = mod r * c if d == prime - 1
        c = square c
      end

      square_roots_given_candidate n, r
    end
This is Algorithm 3.34 from http://cacr.uwaterloo.ca/hac/
 It has no limitations except prime must be at least three, so we could
 remove all the other square root algorithms if we wanted to.