Multiplicative order

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

This task has been clarified. Its programming examples are in need of review to ensure that they still fit the requirements of the task.

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 oder of 37 relative to 1000 is 100 because 37^100 is 1 (modulo 1000), and no number smaller than 100 would do. Note that the multiplicative order is undefined if a and m are not relatively prime.

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.

J

The dyadic verb mo converts it's 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 (this line could do with a more detailed explanation). 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

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