Sam Roberts <sroberts / uniserve.com> wrote in message news:<20041123000430.GA12382 / ensemble.local>...
> Modular reductions, modulare inversions, etc.
> 
> I'm looking through the library docs now, but not seeing anything.
> 
> Cheers,
> Sam

It's far from complete or optimal, and it's nearly undocumented,
but maybe my Imod library might do what you need.  Ultimately, it
should have the functionality of "Mod" in Pari/gp. This generates
"integermods" -- mathematical entities equivalent to residue classes
of integers.  My rules for arithmetic of integermods with different
moduli are derived from Pari/gp.

The Imod library can create integermods in two ways:
a = Imod(7, 13)         # 7 mod 13
a = 7.imod(13)          # 7 mod 13
b = a.inv               # mod 13 reciprocal: b == Imod(2, 13)
a*b                     # => Imod(1, 13)
Imod(5, 13)*Imod(7,13)  # => Imod(9, 13)
b = Imod(4, 12)
b.inv                   # => nil (no inverse) (probably should raise exception)
7.imod(4) == -1.imod(4) # => true

 BTW, if there is some particular functionality you'd like,
feel free to ask.

 Here is what I have so far.

Put the following in a file "imod.rb"
-------------------------------------
#!/usr/bin/env ruby
# -*- ruby -*-
# file imod.rb
# Integermods for ruby
# Uses ordinary integers for modulus == 0

reqfiles = %w(English mathn)
reqfiles.each{|f| require f}

class Integer
# The lcm bundled with mathn (from rational.rb) is broken
  alias :old_lcm :lcm 
  def lcm(b)
    if self == 0 and b == 0 then return 0 end
    old_lcm(b)
  end

  def num
    self
  end

  def modulus
    0
  end

  def factors # prime factors with multiplicities
    a = self.abs; factor_hash = Hash.new
    Prime.new.each { |p|
      while a % p == 0 do
        factor_hash[p] ||= 0 # initialize if it doesn't exist
        factor_hash[p] += 1
        a = a.div(p)
      end
      break if a < p
    }
    factor_hash
  end

  def moebius_mu
    h = factors; mu = 1
    h.each_key{ |k|
      return 0 if h[k] > 1
      mu = -mu
    }
    mu
  end

  def totatives
    t = [1] # 1 is always a totative
    (2..self).each{|k| t << k if self.gcd(k) == 1}
    t
  end

  # Is this really much faster than self.totatives.length ?
  # Can this be rewritten with inject?
  def totient 
    p = 1
    self.factors.each{|k, s| p *= k**s - k**(s-1)}
    p
  end

  alias :euler_phi :totient

  def primitive_roots
    # STUB
  end

  def imod(m)
    Imod(self,m)
  end

end
    
# a & m are integers
def Imod(a, m)
  if m == 0 then a else Imod.new(a, m) end
end

# Duck-type x.is_a?(Integer) as x.responds_to?(modulus) and x.modulus==0
# @modulus.totient is used a lot.  This should be precalculated or at
# least memoized.
class Imod < Numeric
  attr_reader :num, :modulus

# Imod.new should not be used directly: use constructor Imod or n.imod(m)
# instead.  num and modulus are integers: modulus != 0
# In the constructor Imod and integer method imod, a modulus of 0 results
# in an integer

  def initialize(num, modulus)
    raise ArgumentError unless num.kind_of?(Integer)
    raise ArgumentError unless modulus.kind_of?(Integer)
    raise ArgumentError if modulus == 0
    @num = num.modulo(modulus.abs)
    @modulus = modulus
  end

  def coerce(other)
    if other.is_a?(Integer)
      [Imod(other, @modulus), self]
    end
  end

  def unit?
    @num.gcd(@modulus) == 1
  end

  def ==(other)
    other.is_a?(Imod) and  @num == other.num and @modulus == other.modulus
  end

  def +(other)
    if other.is_a?(Integer)
      self + Imod(other, @modulus)
    else # assume Imod
      g = @modulus.gcd(other.modulus)
      Imod(@num + other.num, g)
    end
  end

  def -(other)
    if other.is_a?(Integer)
      self - Imod(other, @modulus)
    else # assume Imod
      g = @modulus.gcd(other.modulus)
      Imod(@num - other.num, g)
    end
  end

  def *(other)
    if other.is_a?(Integer)
      g = @modulus
    else # assume other is Imod
      g = @modulus.gcd(other.modulus)
    end
    return Imod(@num * other.num, g)
  end

  def **(n)
    raise ArgumentError unless n.is_a?(Integer)
    if n < 0
      raise ArgumentError unless self.unit?
      return (self.inv)**(-n)
    else # TODO: better algorithm for the else branch
      n = n.modulo(@modulus.totient)
      p = 1
      n.times do p *= self end
    end
    return p
  end

  # probably should raise exception for non-unit
  def inv
    if self.unit?
      self**(@modulus.totient - 1)
    else
      nil
    end
  end

  # temporary bone-headed implementation
  def square?
    s = false
    (0...@modulus).each{|num|
      t = Imod(num, @modulus)
      if t*t == self
        s = true
        break
      end
    }
    s
  end

  alias :quadratic_residue? :square?

  def primitive_root?
    # STUB
  end
  alias :prim_root? :primitive_root?

  def discrete_log(g)
    # STUB
  end

  alias :ind :discrete_log

  def haupt_exponent
    # STUB
  end

  alias :multiplicative_order :haupt_exponent
  alias :modulo_order :haupt_exponent
  # Would "#{@num}.imod(#{@modulus})" be better?
  def to_s
    "Imod(#{@num}, #{@modulus})"
  end

# Need Legendre-Jacobi symbols and other good stuff
# Need Chinese remainder stuff

end