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)) .

Translation of: python

Works with: ALGOL 68 version Standard - no extensions to language used

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<lang algol> MODE LOOPINT = INT;

PROC gcd = (LONG INT in a, in b)LONG INT: (

   LONG INT a := in a, b := in b;
   WHILE b /= 0 DO LONG INT swap=a; a:=b; b:=swap %* a OD;
   a

);

PROC lcm = (LONG INT a, b)LONG INT : (a*b) OVER gcd(a, b);

MODE YIELDBOOL = PROC(BOOL)VOID;

OP ANY = (PROC(YIELDBOOL)VOID for)BOOL: (

# FOR value IN for DO #
  for((BOOL value)VOID:
      IF value THEN return true FI
  );
  else: FALSE EXIT
  return true: TRUE

);

OP ALL = (PROC(YIELDBOOL)VOID for)BOOL: (

# FOR value IN for DO #
  for((BOOL value)VOID:
      IF NOT value THEN return false FI
  );
  else: TRUE EXIT
  return false: FALSE

);

PROC is prime = (LONG INT p)BOOL:

   ( p > 1 |#ANDF# ALL((YIELDBOOL yield)VOID: factored iterator(p, (LONG INT f, LONG INT e)VOID: yield(f = p))) | FALSE );

FLEX[4]LONG INT prime list := (2,3,5,7);

OP +:= = (REF FLEX[]LONG INT lhs, LONG INT rhs)VOID: (

   [UPB lhs +1] LONG INT next lhs;
   next lhs[:UPB lhs] := lhs;
   lhs := next lhs;
   lhs[UPB lhs] := rhs 

);

PROC primes iterator = (PROC (LONG INT)VOID yield)VOID: (

   LONG INT p;
   FOR p index TO UPB prime list DO
       p:= prime list[p index];
       yield(p)
   OD;
   DO
       p +:= 2;
       WHILE NOT is prime(p) DO
           p +:= 2
       OD;
       prime list +:= p;
       yield(p)
   OD

);

PROC factored iterator = (LONG INT in a, PROC (LONG INT,LONG INT)VOID yield)VOID: (

   LONG INT a := in a;
 # FOR p IN primes iterator DO #
   primes iterator( (LONG INT p)VOID:(
       LONG INT j := 0;
       WHILE a MOD p = 0 DO
           a := a % p;
           j +:= 1
       OD;
       IF j > 0 THEN yield (p,j) FI;
       IF a < p*p THEN done FI
     )
   );
   done:
   IF a > 1 THEN yield (a,1) FI

);

PROC pow mod = (LONG INT b, in e, mod)LONG INT: (

   LONG INT 
       sq := b, 
       e := in e,
       out:= IF ODD e THEN b ELSE 1 FI;
   WHILE
       e := e OVER 2;
       e /= 0 
   DO
       sq := sq * sq %* mod;
       IF ODD e THEN out := out * sq %* mod FI
   OD;
   out

);

PROC range iterator = (LOOPINT lwb, step, upb, YIELDLONGINT yield)VOID:

   FOR i FROM lwb BY step TO upb DO yield(i) OD;

MODE SORTSTRUCT = LONG INT; PR READ "prelude/sort.a68" PR;

PROC mult0rdr1 = (LONG INT a, p, e)LONG INT: (

   LONG INT m := p ** SHORTEN e;
   LONG INT t := (p-1)*(p**SHORTEN (e-1)); #  = Phi(p**e) where p prime #
   LONG INT q;
   FLEX[0]LONG INT qs := (1);
 # FOR f0,f1 in factored(t) DO #
   factored iterator(t, (LONG INT f0,LONG INT f1)VOID: (
           FLEX[SHORTEN((f1+1)*UPB qs)]LONG INT next qs;
           FOR j TO SHORTEN f1 + 1 DO
               FOR q index TO UPB qs DO
                   q := qs[q index];
                   next qs[(j-1)*UPB qs+q index] := q * f0**(j-1)
               OD
           OD;
           qs := next qs
       )
   );
   VOID(in place shell sort(qs));

   FOR q index TO UPB qs DO
       q := qs[q index];
       IF pow mod(a,q,m)=1 THEN done FI
   OD;
   done:
   q

);

MODE YIELDLONGINT = PROC(LONG INT)VOID; MODE LONGINTITERATOR = PROC(YIELDLONGINT)VOID;

PROC int reduce iterator = (PROC (LONG INT,LONG INT)LONG INT diadic, LONGINTITERATOR iterator, LONG INT initial value)LONG INT: (

 LONG INT out := initial value;

1. FOR next IN iterator DO #

 iterator((LONG INT next)VOID:
   out := diadic(out, next)
 );
 out

);

PROC mult order = (LONG INT a, LONG INT m)LONG INT: (

   PROC mofs iterator = (YIELDLONGINT yield)VOID:(
     # FOR r IN factored(m) DO #
       factored iterator(m, (LONG INT p, LONG INT count)VOID:
           yield(mult0rdr1(a,p,count))
       )
   );
   int reduce iterator(lcm, mofs iterator, 1)

);

main:(

   FORMAT d = $g(-0)$;
   printf((d, mult order(37, 1000), $l$));        # 100 #
   LONG INT b := LENG 10**20-1;
   printf((d, mult order(2, b), $l$)); # 3748806900 #
   printf((d, mult order(17,b), $l$)); # 1499522760 #
   b := 100001;
   printf((d, mult order(54,b), $l$));
   printf((d, pow mod( 54, mult order(54,b),b), $l$));
   IF ANY( (YIELDBOOL yield)VOID: FOR r FROM 2 TO SHORTEN mult order(54,b)-1 DO yield(1=pow mod(54,r, b)) OD  )
   THEN
       printf(($g$, "Exists a power r < 9090 where pow mod(54,r,b) = 1", $l$))
   ELSE
       printf(($g$, "Everything checks.", $l$))
   FI

)</lang> Output:

100
3748806900
1499522760
9090
1
Everything checks.

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>