Multiplicative order

Multiplicative order
You are encouraged to solve this task according to the task description, using any language you may know.

The multiplicative order of a relative to m is the least positive integer n such that a^n is 1 (modulo m). For example, the multiplicative order of 37 relative to 1000 is 100 because 37^100 is 1 (modulo 1000), and no number smaller than 100 would do.

One possible algorithm that is efficient also for large numbers is the following: By the Chinese Remainder Theorem, it's enough to calculate the multiplicative order for each prime exponent p^k of m, and combine the results with the least common multiple operation. Now the order of a wrt. to p^k must divide Φ(p^k). Call this number t, and determine it's factors q^e. Since each multiple of the order will also yield 1 when used as exponent for a, it's enough to find the least d such that (q^d)*(t/(q^e)) yields 1 when used as exponent.

Implement a routine to calculate the multiplicative order along these lines. You may assume that routines to determine the factorization into prime powers are available in some library.

---

An algorithm for the multiplicative order can be found in Bach & Shallit, Algorithmic Number Theory, Volume I: Efficient Algorithms, The MIT Press, 1996:

Exercise 5.8, page 115:

Suppose you are given a prime p and a complete factorization of p-1 . Show how to compute the order of an element a in (Z/(p))^* using O((lg p)⁴/(lg lg p)) bit operations.

Solution, page 337:

Let the prime factorization of p-1 be q1^e1q2^e2...qk^ek . We use the following observation: if x^((p-1)/qi^fi) = 1 (mod p) , and fi=ei or x^((p-1)/qi^fi+1) != 1 (mod p) , then qi^ei-fi||ord_p x . (This follows by combining Exercises 5.1 and 2.10.) Hence it suffices to find, for each i , the exponent fi such that the condition above holds.

This can be done as follows: first compute q1^e1, q2^e2, ... , qk^ek . This can be done using O((lg p)²) bit operations. Next, compute y1=(p-1)/q1^e1, ... , yk=(p-1)/qk^ek . This can be done using O((lg p)²) bit operations. Now, using the binary method, compute x1=a^y1(mod p), ... , xk=a^yk(mod p) . This can be done using O(k(lg p)³) bit operations, and k=O((lg p)/(lg lg p)) by Theorem 8.8.10. Finally, for each i , repeatedly raise xi to the qi-th power (mod p) (as many as ei-1 times), checking to see when 1 is obtained. This can be done using O((lg p)³) steps. The total cost is dominated by O(k(lg p)³) , which is O(k(lg p)⁴/(lg lg p)) .

Haskell

Assuming a function to calculate prime power factors,

primeFacsExp :: Integer -> [(Integer, Int)]

and another function

powerMod :: (Integral a, Integral b) => a -> a -> b -> a
powerMod m _ 0 =  1
powerMod m x n | n > 0 = f x' (n-1) x' where
  x' = x `rem` m
  f _ 0 y = y
  f a d y = g a d where
    g b i | even i    = g (b*b `rem` m) (i `quot` 2)
          | otherwise = f b (i-1) (b*y `rem` m)
powerMod m _ _  = error "powerMod: negative exponent"

to efficiently calculate powers modulo some integral, we get

 
 multOrder a m 
   | gcd a m /= 1  = error "Arguments not coprime"
   | otherwise     = foldl1' lcm $ map (multOrder' a) $ primeFacsExp m
 
 multOrder' a (p,k) = r where
   pk = p^k
   t = (p-1)*p^(k-1) -- totient \Phi(p^k)
   r = product $ map find_qd $ primeFacsExp $ t
   find_qd (q,e) = q^d where
     x = powerMod pk a (t `div` (q^e))
     d = length $ takeWhile (/= 1) $ iterate (\y -> powerMod pk y q) x

J

The dyadic verb mo converts its arguments to exact numbers a and m, executes mopk on the factorization of m, and combines the result with the least common multiple operation.

mo=: 4 : 0
 a=. x: x
 m=. x: y
 assert. 1=a+.m
 *./ a mopk"1 |: __ q: m
)

The dyadic verb mopk expects a pair of prime and exponent in the second argument. It sets up a verb pm to calculate powers module p^k. Then calculates Φ(p^k) as t, factorizes it again into q and e, and calculates a^(t/(q^e)) as x. Now, it finds the least d such that subsequent application of pm yields 1. Finally, it combines the exponents q^d into a product.

mopk=: 4 : 0
 a=. x: x
 'p k'=. x: y
 pm=. (p^k)&|@^
 t=. (p-1)*p^k-1  NB. totient
 'q e'=. __ q: t
 x=. a pm t%q^e
 d=. (1<x)+x (pm i. 1:)&> (e-1) */\@$&.> q
 */q^d
)

For example:

   37 mo 1000
100
   2 mo _1+10^80x
190174169488577769580266953193403101748804183400400

Python

<lang python>def gcd(a, b):

   while b != 0:
       a, b = b, a % b
   return a

def lcm(a, b):

   return (a*b) / gcd(a, b)

def isPrime(p):

   return (p > 1) and all(f == p for f,e in factored(p))

primeList = [2,3,5,7] def primes():

   for p in primeList:
       yield p
   while 1:
       p += 2
       while not isPrime(p):
           p += 2
       primeList.append(p)
       yield p

def factored( a):

   for p in primes():
       j = 0
       while a%p == 0:
           a /= p
           j += 1
       if j > 0:
           yield (p,j)
       if a < p*p: break
   if a > 1:
       yield (a,1)
       

def multOrdr1(a,(p,e) ):

   m = p**e
   t = (p-1)*(p**(e-1)) #  = Phi(p**e) where p prime
   qs = [1,]
   for f in factored(t):
       qs = [ q * f[0]**j for j in range(1+f[1]) for q in qs ]
   qs.sort()

   for q in qs:
       if pow( a, q, m )==1: break
   return q

def multOrder(a,m):

   assert gcd(a,m) == 1
   mofs = (multOrdr1(a,r) for r in factored(m))
   return reduce(lcm, mofs, 1)

if __name__ == "__main__":

   print multOrder(37, 1000)        # 100
   b = 10**20-1
   print multOrder(2, b) # 3748806900
   print multOrder(17,b) # 1499522760
   b = 100001
   print multOrder(54,b)
   print pow( 54, multOrder(54,b),b)
   if any( (1==pow(54,r, b)) for r in range(1,multOrder(54,b))):
       print 'Exists a power r < 9090 where pow(54,r,b)==1'
   else:
       print 'Everything checks.'</lang>

Ruby

<lang ruby>require 'rational' # for lcm require 'mathn' # for prime_division

def powerMod(b, p, m)

 result = 1
 bits = p.to_s(2)
 for bit in bits.split()
   result = (result * result) % m
   if bit == '1'
     result = (result * b) % m
   end
 end
 result

end

def multOrder_(a, p, k)

 pk = p ** k
 t = (p - 1) * p ** (k - 1)
 r = 1
 for q, e in t.prime_division
   x = powerMod(a, t / q ** e, pk)
   while x != 1
     r *= q
     x = powerMod(x, q, pk)
   end
 end      
 r

end

def multOrder(a, m)

 m.prime_division.inject(1) {|result, f|
   result.lcm(multOrder_(a, *f))
 }

end

puts multOrder(37, 1000) # 100 b = 10**20-1 puts multOrder(2, b) # 3748806900 puts multOrder(17,b) # 1499522760 b = 100001 puts multOrder(54,b) puts powerMod(54, multOrder(54,b), b) if (1...multOrder(54,b)).any? {|r| powerMod(54, r, b) == 1}

 puts 'Exists a power r < 9090 where powerMod(54,r,b)==1'

else

 puts 'Everything checks.'

end</lang>