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

ALGOL 68

Translation of: python

Works with: ALGOL 68 version Standard - with preludes manually inserted

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

<lang algol>MODE LOOPINT = INT;

MODE POWMODSTRUCT = LONG INT; PR READ "prelude/pow_mod.a68" PR;

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

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; MODE FORBOOL = PROC(YIELDBOOL)VOID;

OP ANY = (FORBOOL sequence)BOOL: (

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

);

OP ALL = (FORBOOL sequence)BOOL: (

# FOR value IN sequence DO #
  sequence((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: for factored(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 for primes = (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 for factored = (LONG INT in a, PROC (LONG INT,LONG INT)VOID yield)VOID: (

   LONG INT a := in a;
 # FOR p IN for primes DO #
   for primes( (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 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 #
   for factored(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 FORLONGINT = PROC(YIELDLONGINT)VOID;

PROC for reduce = (PROC (LONG INT,LONG INT)LONG INT diadic, FORLONGINT 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 for mofs = (YIELDLONGINT yield)VOID:(
     # FOR r IN factored(m) DO #
       for factored(m, (LONG INT p, LONG INT count)VOID:
           yield(mult0rdr1(a,p,count))
       )
   );
   for reduce(lcm, for mofs, 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>

Tcl

Slavishly

Translation of: Python

<lang tcl>package require Tcl 8.5 package require struct::list

proc multOrder {a m} {

   assert {[gcd $a $m] == 1}
   set mofs [list]
   dict for {p e} [factor_num $m] {
       lappend mofs [multOrdr1 $a $p $e]
   }
   return [struct::list fold $mofs 1 lcm]

}

proc multOrdr1 {a p e} {

   set m [expr {$p ** $e}]
   set t [expr {($p - 1) * ($p ** ($e - 1))}]
   set qs [dict create 1 ""]
   
   dict for {f0 f1} [factor_num $t] {
       dict for {q -} $qs {
           foreach j [range [expr {1 + $f1}]] {
               dict set qs [expr {$q * $f0 ** $j}] ""
           }
       }
   }
   
   dict for {q -} $qs {
       if {pypow($a, $q, $m) == 1} break
   }
   return $q    

}

2. utility procs

proc assert {condition {message "Assertion failed!"}} {

   if { ! [uplevel 1 [list expr $condition]]} {
       return -code error $message
   }

}

proc gcd {a b} {

   while {$b != 0} {
       lassign [list $b [expr {$a % $b}]] a b
   }
   return $a

}

proc lcm {a b} {

   expr {$a * $b / [gcd $a $b]}

}

proc factor_num {num} {

   primes::restart
   set factors [dict create]
   for {set i [primes::get_next_prime]} {$i <= $num} {} {
       if {$num % $i == 0} {
           dict incr factors $i
           set num [expr {$num / $i}]
           continue
       } elseif {$i*$i > $num} {
           dict incr factors $num
           break
       } else {
           set i [primes::get_next_prime]
       }
   }
   return $factors

}

2. a range command akin to Python's

proc range args {

   foreach {start stop step} [switch -exact -- [llength $args] {
       1 {concat 0 $args 1}
       2 {concat   $args 1}
       3 {concat   $args  }
       default {error {wrong # of args: should be "range ?start? stop ?step?"}}
   }] break
   if {$step == 0} {error "cannot create a range when step == 0"}
   set range [list]
   while {$step > 0 ? $start < $stop : $stop < $start} {
       lappend range $start
       incr start $step
   }
   return $range

}

1. python's pow()

proc ::tcl::mathfunc::pypow {x y {z ""}} {

   expr {$z eq "" ? $x ** $y : ($x ** $y) % $z}

}

2. prime number generator
3. ref http://wiki.tcl.tk/5996

namespace eval primes {}

proc primes::reset {} {

   variable list [list]
   variable current_index end

}

namespace eval primes {reset}

proc primes::restart {} {

   variable list
   variable current_index
   if {[llength $list] > 0} {
       set current_index 0
   }

}

proc primes::is_prime {candidate} {

   variable list

   foreach prime $list {
       if {$candidate % $prime == 0} {
           return false
       }
       if {$prime * $prime > $candidate} {
           return true
       }
   }
   while true {
       set largest [get_next_prime]
       if {$largest * $largest >= $candidate} {
           return [is_prime $candidate]
       }
   }

}

proc primes::get_next_prime {} {

   variable list
   variable current_index
   
   if {$current_index ne "end"} {
       set p [lindex $list $current_index]
       if {[incr current_index] == [llength $list]} {
           set current_index end
       }
       return $p
   }
   
   switch -exact -- [llength $list] {
       0 {set candidate 2}
       1 {set candidate 3}
       default {
           set candidate [lindex $list end]
           while true {
               incr candidate 2
               if {[is_prime $candidate]} break
           }
       }
   }
   lappend list $candidate
   return $candidate

}

puts [multOrder 37 1000] ;# 100

set b [expr {10**20 - 1}] puts [multOrder 2 $b] ;# 3748806900 puts [multOrder 17 $b] ;# 1499522760

set a 54 set m 100001 puts [set n [multOrder $a $m]] ;# 9090 puts [expr {pypow($a, $n, $m)}] ;# 1

set lambda {{a n m} {expr {pypow($a, $n, $m) == 1}}} foreach r [lreverse [range 1 $n]] {

   if {[apply $lambda $a $r $m]} {
       error "Oops, $n is not the smallest:  {$a $r $m} satisfies $lambda"
   }
   if {$r % 1000 == 0} {puts "$r ..."}

} puts "OK, $n is the smallest n such that {$a $n $m} satisfies $lambda"</lang