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

Ada

Instead of assuming a library call to factorize the modulus, we assume the caller of our Find_Order function has already factorized it. The Multiplicative_Order package is specified as follows ("multiplicative_order.ads"). <lang Ada>package Multiplicative_Order is

  type Positive_Array is array (Positive range <>) of Positive;

  function Find_Order(Element, Modulus: Positive) return Positive;
  -- naive algorithm
  -- returns the smallest I such that (Element**I) mod Modulus = 1

  function Find_Order(Element: Positive;
                      Coprime_Factors: Positive_Array) return Positive;
  -- faster algorithm for the same task
  -- computes the order of all Coprime_Factors(I)
  -- and returns their least common multiple
  -- this gives the same result as Find_Order(Element, Modulus) 
  -- with Modulus being the product of all the Coprime_Factors(I)
  --
  -- preconditions: (1) 1 = GCD(Coprime_Factors(I), Coprime_Factors(J)) 
  --                    for all pairs I, I with I /= J
  --                (2) 1 < Coprime_Factors(I)   for all I

end Multiplicative_Order;</lang>

Here is the implementation ("multiplicative_order.adb"):

<lang Ada>package body Multiplicative_Order is

  function Find_Order(Element, Modulus: Positive) return Positive is

     function Power(Exp, Pow, M: Positive) return Positive is
        -- computes Exp**Pow mod M;
        -- note that Ada's native integer exponentiation "**" may overflow on
        -- computing Exp**Pow before ever computing the "mod M" part
        Result: Positive := 1;
        E: Positive := Exp;
        P: Natural := Pow;
     begin
        while P > 0 loop
           if P mod 2 = 1 then
              Result := (Result * E) mod M;
           end if;
           E := (E * E) mod M;
           P := P / 2;
        end loop;
        return Result;
     end Power;

  begin -- Find_Order(Element, Modulus)
     for I in 1 .. Modulus loop
        if Power(Element, I, Modulus) = 1 then
           return Positive(I);
        end if;
     end loop;
     raise Program_Error with
       Positive'Image(Element) &" is not coprime to" &Positive'Image(Modulus);
  end Find_Order;

  function Find_Order(Element: Positive;
                      Coprime_Factors: Positive_Array) return Positive is

        function GCD (A, B : Positive) return Integer is
           M : Natural := A;
           N : Natural := B;
           T : Natural;
        begin
           while N /= 0 loop
              T := M;
              M := N;
              N ;:= T mod N;
           end loop;
           return M;
        end GCD; -- from http://rosettacode.org/wiki/Least_common_multiple#Ada

        function LCM (A, B : Natural) return Integer is
        begin
           if A = 0 or B = 0 then
              return 0;
           end if;
           return abs (A * B) / Gcd (A, B);
        end LCM; -- from http://rosettacode.org/wiki/Least_common_multiple#Ada

        Result : Positive := 1;

  begin -- Find_Order(Element, Coprime_Factors)
     for I in Coprime_Factors'Range loop
        Result := LCM(Result, Find_Order(Element, Coprime_Factors(I)));
     end loop;
     return Result;
  end Find_Order;

end Multiplicative_Order;</lang>

This is a sample program using the Multiplicative_Order package: <lang Ada>with Ada.Text_IO, Multiplicative_Order;

procedure Main is

  package IIO is new Ada.Text_IO.Integer_IO(Integer);
  use Multiplicative_Order;

begin

  IIO.Put(Find_Order(3,10));
  IIO.Put(Find_Order(37,1000));
  IIO.Put(Find_Order(37,10_000));
  IIO.Put(Find_Order(37, 3343));
  IIO.Put(Find_Order(37, 3344));
  -- IIO.Put(Find_Order( 2,1000));
    --would raise Program_Error, because there is no I with 2**I=1 mod 1000
  Ada.Text_IO.New_Line;
  IIO.Put(Find_Order(3, (2,5)));           --  3 *   5 = 10
  IIO.Put(Find_Order(37, (8, 125)));       --  8 * 125 = 1000
  IIO.Put(Find_Order(37, (16, 625)));      -- 16 * 625 = 10_000
  IIO.Put(Find_Order(37, (1 => 3343)));    -- 1-element-array: 3343 is a prime
  IIO.Put(Find_Order(37, (11, 19, 16)));   -- 11 * 19 * 16 = 3344

  -- the following gives a wrong result, because 8 and 2 are not coprime
  IIO.Put(Find_Order(37, (11, 19, 8, 2)));

end Main;</lang>

The output from the sample program:

          4        100        500       1114         20
          4        100        500       1114         20         10

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 algol68>MODE LOOPINT = INT;

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

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

MODE GCDSTRUCT = LONG INT; PR READ "prelude/gcd.a68" PR;

PR READ "prelude/iterator.a68" PR;

PROC is prime = (LONG INT p)BOOL:

   ( p > 1 |#ANDF# ALL((YIELDBOOL yield)VOID: 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 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 factored = (LONG INT in a, PROC (LONG INT,LONG INT)VOID yield)VOID: (

   LONG INT a := in a;
 # FOR          p IN  # primes( # DO #
      (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
     )
 # ) OD #  );
   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 #, 
      (LONG INT f0,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
       )
 #   OD # );
   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

);

PROC 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 #

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

);

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

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

<lang haskell>primeFacsExp :: Integer -> [(Integer, Int)]</lang>

to calculate prime power factors, and another function

<lang haskell>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"</lang>

to efficiently calculate powers modulo some Integral, we get

<lang haskell>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</lang>

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.

<lang j>mo=: 4 : 0

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

)</lang>

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.

<lang j>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

)</lang>

For example:

<lang j> 37 mo 1000 100

  2 mo _1+10^80x

190174169488577769580266953193403101748804183400400</lang>

Mathematica

In Mathematica this is really easy, as this function is built-in: <lang Mathematica>MultiplicativeOrder[k,n]

   gives the multiplicative order of k modulo n, defined as the smallest integer m such that k^m == 1 mod n.

MultiplicativeOrder[k,n,{r1,r2,...}]

   gives the generalized multiplicative order of k modulo n, defined as the smallest integer m such that k^m==ri mod n for some i.</lang>

Examples: <lang Mathematica>MultiplicativeOrder[37, 1000] MultiplicativeOrder[10^100 + 1, 7919] (*10^3th prime number Prime[1000]*) MultiplicativeOrder[10^1000 + 1, 15485863] (*10^6th prime number*) MultiplicativeOrder[10^10000 - 1, 22801763489] (*10^9th prime number*) MultiplicativeOrder[13, 1 + 10^80] MultiplicativeOrder[11, 1 + 10^100]</lang> gives back: <lang Mathematica>100 3959 15485862 22801763488 109609547199756140150989321269669269476675495992554276140800 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800</lang>

PARI/GP

<lang>znorder(Mod(a,n))</lang>

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

Translation of: Python

Library: Tcllib (Package: struct::list)

<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