# runpaint/numb

Integer#ord, Integer#power_mod, Integer#coprime: Implement

1 parent 71f518c commit 4ce1003d21b7a7ada1cc21b3e7fd14237f92ce42 committed Oct 16, 2010
Showing with 183 additions and 4 deletions.
1. +22 −4 lib/numb/divisors.rb
2. +29 −0 lib/numb/ord.rb
3. +15 −0 lib/numb/power_mod.rb
4. +72 −0 spec/numb/ord_spec.rb
5. +45 −0 spec/numb/power_mod_spec.rb
 @@ -65,6 +65,22 @@ def coprime?(x) alias :⊥ :coprime? alias :stranger? :coprime? + # An enumeration of numbers coprime to `x` from `self` onward. If a + # block is given, it is yielded to with the next number in the + # sequence; otherwise, an `Enumerator` is returned. + # + # 4.coprime(3).first(5) #=> [4, 5, 7, 8, 10] + # + # @param [Integer] x each number in the enumeration is coprime with `x` + # @yield [Integer] n next number that is >= `self` and coprime with `x` + # @return [Enumerator] returned if the block is omitted + def coprime(x) + return enum_for(__method__, x) unless block_given? + (self..Float::INFINITY).each do |n| + yield n if n.coprime?(x) + end + end + # A deficient number is a number n for which σ(n) < 2n. That is, the sum of # its divisors are less than the number. (To calculate the sum of divisors # for an arbitrary integer see Integer#σ). @@ -347,18 +363,20 @@ def first_with_n_divisors limit = Prime.first(pf.size).zip(pf.reverse).map{|b,e| b**(e-1)}.reduce(:*) neighbour = ->(e) { (self/e).divisors.reject{|d| d > e} } x, div, exponents = self, divisors, {} - Prime.each do |b| - max_exponent = Math.log(limit, b).floor - d = div.reject{|d| d > max_exponent}.sort.reverse - [1] + d = begin + max_exponent = Math.log(limit, b).floor + div.reject{|d| d > max_exponent} + rescue FloatDomainError + div + end.sort.reverse - [1] unless b == 2 prev_neighbours = exponents[exponents.keys.last].values.flatten d.reject!{|e, n| not prev_neighbours.include?(e)} end exponents[b] = Hash[d.map{|e| [e, neighbour.(e)]}] break if (x /= b) < 1 end - complete_chain = ->(b, e, goal=self, chain=nil) do chain ||= {chain: [[b, e]]} return chain unless exponents.key?(b) and exponents[b].key?(e)
 @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- +class Integer + # Returns the multiplicative order of `self` % `m`, or `nil` if + # `self` is not coprime to `m` + # + # Given an integer `self` and a positive integer `m` with gcd(`a`, + # `m`) = 1, the multiplicative order of `self` modulo `m` is the + # smallest positive integer `k` with: `self`^`k` ≡ 1 (modulo `m`). + # + # @param [Integer] m the modulus + # @return [Integer, nil] the power, `k`, or `nil` if `self` and `m` are not coprime + def ord(m) + return unless coprime?(m) + m.prime_division.inject(1) do |result, f| + (p, k), r = f, 1 + pk = p ** k + # We could calculate the totient here as `(p - 1) * p ** (k - + # 1)`, but it feels cleaner to separate the logic + (t = pk.φ).prime_division.each do |q, e| + x = power_mod(t / q ** e, pk) + while x != 1 + r *= q + x = x.power_mod(q, pk) + end + end + result.lcm(r) + end + end +end
 @@ -0,0 +1,15 @@ +class Integer + # `self`^`b` mod `m` + # + # @param [Integer] b power to which `self` should be raised + # @param [Integer] m modulus + # @return [Integer] `self`^`b` % `m` + def power_mod(b, m) + result = 1 + b.to_s(2).chars.each do |bit| + result = (result * result) % m + result = (result * self) % m if bit==?1 + end + result + end +end
 @@ -0,0 +1,72 @@ +# coding: utf-8 +describe Integer, "#ord" do + @seq = { + # A002326 + 2 => [1,2,4,3,6,10,12,4,8,18,6,11,20,18,28,5,10,12,36, + 12,20,14,12,23,21,8,52,20,18,58,60,6,12,66,22,35, + 9,20,30,39,54,82,8,28,11,12,10,36,48,30,100,51,12, + 106,36,36,28,44,12,24,110,20,100,7,14,130,18,36, + 68,138,46,60,28].map.with_index{|k, i| [2*i + 1, k]}, + # A050975 + 3 => 2.coprime(3).first(74).zip([1,2,4,6,2,4,5,3,6,4,16,18,4,5, + 11,20,3,6,28,30,8,16,12,18,18,4,8,42,10,11,23,42,20,6,52,20, + 6,28,29, 10,30,16,12,22,16,12,35,12,18,18,30,78,4,8,41,16, + 42,10,88,6,22,23,36,48,42,20,100,34,6,52,53,27,20, 12]), + # A050976 + 4 => 2.coprime(4).first(78).zip([1,2,3,3,5,6,2,4,9,3,11,10,9,14,5,5,6,18,6,10,7,6, + 23,21,4,26,10,9,29,30,3,6,33,11,35,9,10,15,39,27, + 41,4,14,11,6,5,18,24,15,50,51,6,53,18,18,14,22,6, + 12,55,10,50,7,7,65,9,18,34,69,23,30,14,21,74,15, 12,10,26]), + # A050977 + 5 => 2.coprime(5).first(81).zip([1,2,1,2,6,2,6,5,2,4,6,4,16,6,9,6,5,22,2,4,18,6, + 14,3,8,10,16,6,36,9,4,20,6,42,5,22,46,4,42,16,4, + 52,18,6,18,14,29,30,3,6,16,10,22,16,22,5,6,72,36, + 9,30,4,39,54,20,82,6,42,14,10,44,12,22,6,46,8,96, 42,30,25,16]), + + # A050978 + 6 => 2.coprime(6).first(67).zip([1,2,10,12,16,9,11,5,14,6,2,4,40,3,23,14,26,10,58, + 60,12,33,35,36,10,78,82,16,88,12,9,12,10,102,106, + 108,112,11,16,110,25,126,130,18,136,23,60,14,37, + 150,6,156,22,27,83,156,43,10,178,60,4,80,19,96,14, 198,14]), + + # A050979 + 7 => 2.coprime(7).first(81).zip([1,1,2,4,1,2,3,4,10,2,12,4,2,16,3,3,4,10,22,2,4, + 12,9,7,4,15,4,10,16,6,9,3,12,4,40,6,10,12,22,23,2, + 4,16,12,26,9,20,3,7,29,4,60,15,8,12,10,66,16,22, + 70,6,24,9,4,6,12,78,4,27,40,41,16,6,7,10,88,12,22, 15,23,12]), + # A050980 + 8 => 2.coprime(8).first(75).zip([2,4,1,2,10,4,4,8,6,2,11,20,6,28,5,10,4,12,4,20, + 14,4,23,7,8,52,20,6,58,20,2,4,22,22,35,3,20,10,13, + 18,82,8,28,11,4,10,12,16,10,100,17,4,106,12,12,28, + 44,4,8,110,20,100,7,14,130,6,12,68,46,46,20,28,14, 148,5]), + + # A050981 + 9 => 2.coprime(9).first(82).zip([1,1,2,3,1,2,5,3,3,2,8,9,2,5,11,10,3,3,14,15,4,8, + 6,9,9,2,4,21,5,11,23,21,10,3,26,10,3,14,29,5,15,8, + 6,11,8,6,35,6,9,9,15,39,2,4,41,8,21,5,44,3,11,23, + 18,24,21,10,50,17,3,26,53,27,10,6,56,22,14,29,24, 5,5,15]), + + # A002329 + 10 => 2.coprime(10).first(52).zip([1,6,1,2,6,16,18,6,22,3,28,15,2,3,6,5,21,46,42,16, + 13,18,58,60,6,33,22,35,8,6,13,9,41,28,44,6,15,96, + 2,4,34,53,108,3,112,6,48,22,5,42,21,130]) + + } + + @seq.each do |a, seq| + seq.each do |n, k| + it "returns #{k} for #{a}.ord(#{n})" do + a.ord(n).should == k + end + end + + it "returns nil for #{a}.ord(#{a})" do + a.ord(a).should be_nil + end + + b = a * rand(100) + it "returns nil for #{a}.ord(#{b})" do + a.ord(b).should be_nil + end + end +end
 @@ -0,0 +1,45 @@ +# coding: utf-8 +describe Integer, "#power_mod" do + @seq = { + # A015910 + 2 => [0,0,2,0,2,4,2,0,8,4,2,4,2,4,8,0,2,10,2,16,8,4,2, + 16,7,4,26,16,2,4,2,0,8,4,18,28,2,4,8,16,2,22,2,16, + 17,4,2,16,30,24,8,16,2,28,43,32,8,4,2,16,2,4,8,0, + 32,64,2,16,8,44,2,64,2,4,68,16,18,64,2,16,80,4,2, + 64].map.with_index(1){|pm, b| [pm, b, b]}, + + # A002380 + 3 => [0,1,1,3,1,19,25,11,161,227,681,1019,3057,5075, + 15225,29291,55105,34243,233801,439259,269201, + 1856179,3471385,6219851,1882337,5647011,50495465, + 17268667,186023729,21200275,63600825,1264544299, + 3793632897,7085931395].map.with_index{|pm, b| [pm, b, 2**b]}, + + # A056969 + 10 => [0,0,1,0,0,4,3,0,1,0,10,4,10,2,10,0,10,10,10,0,13, + 12,10,16,0,22,1,4,10,10,10,0,10,32,5,28,10,24,25, + 0,10,22,10,12,10,8,10,16,31,0,31,16,10,28,10,16, + 31,42,10,40,10,38,55,0,30,34,10,4,34,60,10,64,10, + 26,25,44,54,40].map.with_index(1){|pm, b| [pm, b, b]} + [ + 0,0,1,0,0,4,4,0,1,0,1,4,3,4,10,0,4,10,9,0,4,12, + 13,16,0,16,10,4,16,10,5,0,1,4,25,28,10,28,16,0,1, + 4,31,12,10,36,27,16,11,0,4,16,28,10,45,32,28,16, + 16,40,47,36,46,0,55,34,10,4,13,60,20,64,72,10,25, + 28,67,16,67,0].map.with_index(1){|pm, b| [pm, 100, b]}, + + # A116609 + 13 => [0,1,1,1,3,1,6,1,1,9,2,1,0,1,7,1,13,1,13,1,13,15, + 13,1,18,13,1,1,13,19,13,1,19,33,27,1,13,17,13,1, + 13,1,13,5,28,31,13,1,48,49,4,13,13,1,32,1,31,53, + 13,1,13,45,55,1,13,31,13,1,58,29,13,1,13,21,7,61, + 62,13,13,1,1,5].map.with_index(1){|pm, b| [pm, b, b]} + } + + @seq.each do |a, seq| + seq.each do |pm, b, m| + it "returns #{pm} for #{a}.power_mod(#{b}, #{m})" do + a.power_mod(b, m).should == pm + end + end + end +end