Permalink
Browse files

Integer#kronecker: Add

  • Loading branch information...
1 parent 12179fb commit b382526212bbb1aa5a4a4d98af8725c87c6be3fc @runpaint committed Oct 17, 2010
Showing with 159 additions and 0 deletions.
  1. +83 −0 lib/numb/kronecker.rb
  2. +76 −0 spec/numb/kronecker_spec.rb
View
@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+class Integer
+ # Returns the Kronecker symbol for `self`/`b`
+ #
+ # The Kronecker symbol is an extension of the Jacobi symbol to all
+ # integers.
+ #
+ # 2.kronecker 3 #=> -1
+ # -8.kronecker 3 #=> 1
+ # 6.kronecker 7 #=> -1
+ #
+ # @param [Integer] b The "denominator"
+ # @return [Integer] -1, 0, 1
+ def kronecker b
+ # The following algorithm is excerpted from Henri Cohen's A Course
+ # in Computational Algebraic Number Theory, 3rd. ed. It appears on
+ # page 48, and is entitled "Algorithm 1.4.12 (Kronecker-Binary).
+
+ # Given a, b ∈ ℤ, this algorithm computes the Kronecker symbol
+ # (a/b) (hence the Legendre symbol when b is an odd prime).
+ raise ArgumentError unless b.is_a?(Integer)
+
+ a = self
+ # The following Array is a lookup table for computing (-1)^(a^2 -
+ # 1)/8 (Cohen, p. 29)
+ tab2 = [0, 1, 0, -1, 0, -1, 0, 1]
+
+ # 1. [Test b = 0] If b = 0, then output 0 if |a| ≠ 1, 1 if |a| = 1
+ # and terminate the algorithm
+ return a.abs == 1 ? 1 : 0 if b.zero?
+
+ # 2. [Remove 2's from b] If a and b are both even, output 0 and
+ # terminate the algorithm. Otherwise, set v ← 0 and while b is
+ # even set v ← v + 1 and b ← b/2. Then if v is even set k ← 1,
+ # otherwise set k ← (-1)^(a^2 - 1)/8 (by table lookup, not by
+ # computing (a^2 - 1)/8). Finally, if b < 0 set b ← -b, and if in
+ # addition a < 0 set k ← -k.
+ return 0 if a.even? and b.even?
+ v, k = 0, nil
+ while b.even?
+ v += 1
+ b = b.quo(2).to_i
+ end
+ k = v.even? ? 1 : tab2[a&7]
+ if b < 0
+ b = -b
+ k =-k if a < 0
+ end
+
+ # 3. [Reduce size once] (Here b is odd and b > 0.) Set a ← a mod b
+ raise unless b.odd? and b > 0
+ a %= b
+
+ loop do
+ # 4. [Finished?] If a = 0, output 0 if b > 1, k if b = 1, and
+ # terminate the algorithm
+ return (b > 1) ? 0 : k if a.zero?
+
+ # 5. [Remove powers of 2] Set v ← 0 and, while a is even, set v ←
+ # v + 1 and a ← a/2. If v is odd, set k ← (-1)^((b^2 - 1)/8) * k
+ v = 0
+ while a.even?
+ v += 1
+ a = a.quo(2).to_i
+ end
+ k = ((-1)**((b**2).pred.quo 8) * k).to_i if v.odd?
+
+ # 6. [Subtract and apply reciprocity] (Here a and b are odd.) Set
+ # r ← b - a. If r > 0, then set k ← (-1)^((a-1)(b-1)/4) * k (using
+ # if statements), b ← a and a ← r, else set a ← -r. Go to step 4.
+ raise unless a.odd? and b.odd?
+ r = b - a
+ if r > 0
+ # The optimisation below is from Cohen, p. 29
+ k = -k unless (a & b & 2).zero?
+ b = a
+ a = r
+ else
+ a = -r
+ end
+ end
+ end
+end
@@ -0,0 +1,76 @@
+describe Integer, "#kronecker" do
+ @seq = {
+ # A109017
+ -6 => [0,1,0,0,0,1,0,1,0,0,0,1,0,-1,0,0,0,-1,0,-1,0,0,0,
+ -1,0,1,0,0,0,1,0,1,0,0,0,1,0,-1,0,0,0,-1,0,-1,0,0,
+ 0,-1,0,1,0,0,0,1,0,1,0,0,0,1,0,-1,0,0,0,-1,0,-1,0,
+ 0,0,-1,0,1,0,0,0,1,0,1,0,0,0,1,0,-1,0,0,0,-1,0,-1,
+ 0,0,0,-1,0,1,0,0,0,1,0,1,0].map.with_index,
+ # -5..-2 are extracted from <http://mathworld.wolfram.com/KroneckerSymbol.html>
+ -5 => [1, -1, 1, 1, 0, -1, 1, -1, 1, 0, -1, 1, -1, -1, 0, 1, -1,
+ -1, -1, 0].map.with_index(1),
+ -4 => [1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0,
+ -1, 0].map.with_index(1),
+ -3 => [1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0, 1, -1, 0,
+ 1, -1].map.with_index(1),
+ -2 => [1, 0, 1, 0, -1, 0, -1, 0, 1, 0, 1, 0, -1, 0, -1, 0, 1, 0,
+ 1, 0].map.with_index(1),
+ # A034947
+ -1 => [1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,-1,-1,1,1,1,-1,1,
+ 1,-1,-1,-1,1,1,-1,-1,1,-1,-1,1,1,1,-1,1,1,-1,-1,1,
+ 1,1,-1,-1,1,-1,-1,-1,1,1,-1,1,1,-1,-1,-1,1,1,-1,
+ -1,1,-1,-1,1,1,1,-1,1,1,-1,-1,1,1,1,-1,-1,1,-1,-1,
+ 1,1].map.with_index(1),
+ # 0..1 are extracted from <http://mathworld.wolfram.com/KroneckerSymbol.html>
+ 0 => [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0].map.with_index(1),
+ 1 => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1].map.with_index(1),
+ # A091337
+ 2 => [1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,
+ 0,1,0,1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1,0,-1,
+ 0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1,
+ 0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,
+ 0,1,0,-1,0,-1,0,1,0,1,0,-1,0,-1,0,1,0,1].map.with_index(1),
+ # A091338
+ 3 => [1,-1,0,1,-1,0,-1,-1,0,1,1,0,1,1,0,1,-1,0,-1,-1,0,
+ -1,1,0,1,-1,0,-1,-1,0,-1,-1,0,1,1,0,1,1,0,1,-1,0,
+ -1,1,0,-1,1,0,1,-1,0,1,-1,0,-1,1,0,1,1,0,1,1,0,1,
+ -1,0,-1,-1,0,-1,1,0,1,-1,0,-1,-1,0,-1,-1,0,1,1,0,
+ 1,1,0,-1,-1,0,-1,1,0,-1,1,0,1,-1,0,1,-1,0].map.with_index(1),
+ # A080891
+ 5 => [0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,
+ 0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,
+ 1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,
+ -1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,
+ -1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0,1,-1,-1,1,0].map.with_index,
+ # The below is extracted from <http://mathworld.wolfram.com/KroneckerSymbol.html>
+ 6 => [1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, -1, 0, 0, 0, -1, 0, 1,
+ 0].map.with_index(1)
+ }
+
+ @seq.each do |a, bs|
+ bs.each do |k, b|
+ it "returns #{k} for #{a}.kronecker(#{b})" do
+ a.kronecker(b).should == k
+ end
+ end
+ end
+
+ # A071935, A071936
+ [1,-1,1,1,1,-1,1,1,1,-1,-1,1,1,-1,1,1,1,-1,1,1,1,
+ -1,-1,1,1,-1,-1,1,1,-1,1,1,1,-1,1,1,1,-1,1,1,1,-1,
+ -1,1,1,-1,-1,1,1,-1,1,1,1,-1,-1,1,1,-1,-1,1,1,-1,
+ 1,1,1,-1,1,1,1,-1,1,1,1,-1,-1,1,1,-1,1,1,1,-1,1,1,
+ 1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,1,1,1,-1,1,1,1,
+ -1,-1].map.with_index(1){|k, n| [n, n.succ, k]} + \
+ [1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,
+ 1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,
+ 1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,
+ 1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,
+ 1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,-1,1,
+ 1,1,1,1,1].map.with_index(1){|k, n| [n.succ, n, k]}.each do |a, b, k|
+ it "returns #{k} for #{a}.kronecker(#{b})" do
+ a.kronecker(b).should == k
+ end
+ end
+end

0 comments on commit b382526

Please sign in to comment.