Multiplicative order

From Rosetta Code
Revision as of 05:34, 11 August 2016 by rosettacode>Purple24 (→‎{{header|Ruby}}: corrected require (deprecated 'mathn'), add output)
Task
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)4/(lg lg p)) bit operations.

Solution, page 337:

Let the prime factorization of p-1 be q1e1q2e2...qkek . We use the following observation: if x^((p-1)/qifi) = 1 (mod p) , and fi=ei or x^((p-1)/qifi+1) != 1 (mod p) , then qiei-fi||ordp 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 q1e1, q2e2, ... , qkek . This can be done using O((lg p)2) bit operations. Next, compute y1=(p-1)/q1e1, ... , yk=(p-1)/qkek . This can be done using O((lg p)2) bit operations. Now, using the binary method, compute x1=ay1(mod p), ... , xk=ayk(mod p) . This can be done using O(k(lg p)3) 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)3) steps. The total cost is dominated by O(k(lg p)3) , which is O((lg p)4/(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, J 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

  -- this violates the precondition, because 8 and 2 are not coprime
  -- it gives an incorrect result
  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.

C

Uses prime/factor functions from Factors of an integer#Prime factoring. This implementation is not robust because of integer overflows. To properly deal with even moderately large numbers, an arbitrary precision integer package is a must. <lang c>ulong mpow(ulong a, ulong p, ulong m) { ulong r = 1; while (p) { if ((1 & p)) r = r * a % m; a = a * a % m; p >>= 1; } return r; }

ulong ipow(ulong a, ulong p) { ulong r = 1; while (p) { if ((1 & p)) r = r * a; a *= a; p >>= 1; } return r; }

ulong gcd(ulong m, ulong n) { ulong t; while (m) { t = m; m = n % m; n = t; } return n; }

ulong lcm(ulong m, ulong n) { ulong g = gcd(m, n); return m / g * n; }

ulong multi_order_p(ulong a, ulong p, ulong e) { ulong fac[10000]; ulong m = ipow(p, e); ulong t = m / p * (p - 1); int i, len = get_factors(t, fac); for (i = 0; i < len; i++) if (mpow(a, fac[i], m) == 1) return fac[i]; return 0; }

ulong multi_order(ulong a, ulong m) { prime_factor pf[100]; int i, len = get_prime_factors(m, pf); ulong res = 1; for (i = 0; i < len; i++) res = lcm(res, multi_order_p(a, pf[i].p, pf[i].e)); return res; }

int main() { sieve(); printf("%lu\n", multi_order(37, 1000)); printf("%lu\n", multi_order(54, 100001)); return 0; }</lang>

EchoLisp

<lang scheme> (require 'bigint)

factor-exp returns a list ((p k) ..)
a = p1^k1 * p2^k2 ..

(define (factor-exp a) (map (lambda (g) (list (first g) (length g))) (group* (prime-factors a))))

copied from Ruby

(define (_mult_order a p k (x)) (define pk (expt p k)) (define t (* (1- p) (expt p (1- k)))) (define r 1) (for [((q e) (factor-exp t))] (set! x (powmod a (/ t (expt q e)) pk)) (while (!= x 1) (*= r q) (set! x (powmod x q pk)))) r)

(define (order a m) 

       "multiplicative order : (order a m) →  n : a^n = 1 (mod m)"

(assert (= 1 (gcd a m)) "a and m must be coprimes") (define mopks (for/list [((p k) (factor-exp m))] (_mult_order a p k))) (for/fold (n 1) ((mopk mopks)) (lcm n mopk)))

results

order 37 1000) 

  → 100

(order (+ (expt 10 100) 1) 7919) 

  → 3959

(order (+ (expt 10 1000) 1) 15485863) 

  → 15485862

</lang>

Go

<lang go>package main

import ( 

   "fmt"
   "math/big"

)

func main() { 

   moTest(big.NewInt(37), big.NewInt(3343))
   b := big.NewInt(100)
   moTest(b.Add(b.Exp(ten, b, nil), one), big.NewInt(7919))
   moTest(b.Add(b.Exp(ten, b.SetInt64(1000), nil), one), big.NewInt(15485863))
   moTest(b.Sub(b.Exp(ten, b.SetInt64(10000), nil), one),
       big.NewInt(22801763489))

   moTest(big.NewInt(1511678068), big.NewInt(7379191741))
   moTest(big.NewInt(3047753288), big.NewInt(2257683301))

}

func moTest(a, n *big.Int) { 

   if a.BitLen() < 100 {
       fmt.Printf("ord(%v)", a)
   } else {
       fmt.Print("ord([big])")
   }
   if n.BitLen() < 100 {
       fmt.Printf(" mod %v ", n)
   } else {
       fmt.Print(" mod [big] ")
   }
   if !n.ProbablyPrime(20) {
       fmt.Println("not computed.  modulus must be prime for this algorithm.")
       return
   }
   fmt.Println("=", moBachShallit58(a, n, factor(new(big.Int).Sub(n, one))))

}

var one = big.NewInt(1) var two = big.NewInt(2) var ten = big.NewInt(10)

func moBachShallit58(a, n *big.Int, pf []pExp) *big.Int { 

   n1 := new(big.Int).Sub(n, one)
   var x, y, o1, g big.Int
   mo := big.NewInt(1)
   for _, pe := range pf {
       y.Quo(n1, y.Exp(pe.prime, big.NewInt(pe.exp), nil))
       var o int64
       for x.Exp(a, &y, n); x.Cmp(one) > 0; o++ {
           x.Exp(&x, pe.prime, n)
       }
       o1.Exp(pe.prime, o1.SetInt64(o), nil)
       mo.Mul(mo, o1.Quo(&o1, g.GCD(nil, nil, mo, &o1)))
   }
   return mo

}

type pExp struct { 

   prime *big.Int
   exp   int64

}

func factor(n *big.Int) (pf []pExp) { 

   var e int64
   for ; n.Bit(int(e)) == 0; e++ {
   }
   if e > 0 {
       n.Rsh(n, uint(e))
       pf = []pExpTemplate:Big.NewInt(2), e
   }
   s := sqrt(n)
   q, r := new(big.Int), new(big.Int)
   for d := big.NewInt(3); n.Cmp(one) > 0; d.Add(d, two) {
       if d.Cmp(s) > 0 {
           d.Set(n)
       }
       for e = 0; ; e++ {
           q.QuoRem(n, d, r)
           if r.BitLen() > 0 {
               break
           }
           n.Set(q)
       }
       if e > 0 {
           pf = append(pf, pExp{new(big.Int).Set(d), e})
           s = sqrt(n)
       }
   }
   return

}

func sqrt(n *big.Int) *big.Int { 

   a := new(big.Int)
   for b := new(big.Int).Set(n); ; {
       a.Set(b)
       b.Rsh(b.Add(b.Quo(n, a), a), 1)
       if b.Cmp(a) >= 0 {
           return a
       }
   }
   return a.SetInt64(0)

}</lang>

Output:
ord(37) mod 3343 = 1114
ord([big]) mod 7919 = 3959
ord([big]) mod 15485863 = 15485862
ord([big]) mod 22801763489 = 22801763488
ord(1511678068) mod 7379191741 = 614932645
ord(3047753288) mod 2257683301 = 62713425

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>

Julia

(Uses the factors function from Factors of an integer#Julia.) <lang julia>function factors(n) 

   f = [one(n)]
   for (p,e) in factor(n)
       f = reduce(vcat, f, [f*p^j for j in 1:e])
   end
   return length(f) == 1 ? [one(n), n] : sort!(f)

end

function multorder(a, m) 

   gcd(a,m) == 1 || error("$a and $m are not coprime")
   res = one(m)
   for (p,e) in factor(m)
       m = p^e
       t = div(m, p) * (p-1)        
       for f in factors(t)
           if powermod(a, f, m) == 1
               res = lcm(res, f)
               break
           end
       end
   end
   res

end</lang>

Example output (using big to employ arbitrary-precision arithmetic where needed): 

julia> multorder(37, 1000)
100

julia> multorder(big(10)^100 + 1, 7919)
3959

julia> multorder(big(10)^1000 + 1, 15485863)
15485862

julia> multorder(big(10)^10000 - 1, 22801763489)
22801763488

Maple

<lang Maple>numtheory:-order( a, n )</lang> For example, <lang Maple>> numtheory:-order( 37, 1000 ); 

                      100</lang>

Mathematica

In Mathematica this is really easy, as this function is built-in: 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.
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: 

100
3959
15485862
22801763488
109609547199756140150989321269669269476675495992554276140800
2583496112724752500580158969425549088007844580826869433740066152289289764829816356800

Maxima

<lang maxima>zn_order(37, 1000); /* 100 */

zn_order(10^100 + 1, 7919); /* 3959 */

zn_order(10^1000 + 1, 15485863); /* 15485862 */

zn_order(10^10000 - 1, 22801763489); /* 22801763488 */

zn_order(13, 1 + 10^80); /* 109609547199756140150989321269669269476675495992554276140800 */

zn_order(11, 1 + 10^100); /* 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800 */</lang>

PARI/GP

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

Perl

Using modules:

Library: ntheory

<lang perl>use ntheory qw/znorder/; say znorder(54, 100001); use bigint; say znorder(11, 1 + 10**100);</lang> or <lang perl>use Math::Pari qw/znorder Mod/; say znorder(Mod(54, 100001)); say znorder(Mod(11, 1 + Math::Pari::PARI(10)**100));</lang>

Perl 6

<lang perl6>my @primes := 2, grep *.is-prime, (3,5,7...*);

sub factor($a is copy) { 

   gather {
       for @primes -> $p {
           my $j = 0;
           while $a %% $p {
               $a div= $p;
               $j++;
           }
           take $p => $j if $j > 0;
           last if $a < $p * $p;
       }
       
       take $a => 1 if $a > 1;
   }

}

sub mo-prime($a, $p, $e) { 

   my $m = $p ** $e;
   my $t = ($p - 1) * ($p ** ($e - 1)); #  = Phi($p**$e) where $p prime
   my @qs = 1;
   for factor($t) -> $f {
       @qs = @qs.map(-> $q { (0..$f.value).map(-> $j { $q * $f.key ** $j }) });
   }
   @qs.sort();

   @qs.first(-> $q { expmod( $a, $q, $m ) == 1});

}

sub mo($a, $m) { 

   $a gcd $m == 1 || die "$a and $m are not relatively prime";
   [lcm] 1, factor($m).map(-> $r { mo-prime($a, $r.key, $r.value) });

}

sub MAIN("test") { 

   use Test;
   
   for (10, 21, 25, 150, 1231, 123141, 34131) -> $n {
       # say factor($n).perl;
       # say factor($n).map(-> $pair { $pair.key ** $pair.value }).perl;
       is ([*] factor($n).map(-> $pair { $pair.key ** $pair.value })), $n, "$n factors correctly";
   }
   
   is mo(37, 1000), 100, 'mo(37,1000) == 100';
   my $b = 10**20-1;
   is mo(2, $b), 3748806900, 'mo(2,10**20-1) == 3748806900';
   is mo(17, $b), 1499522760, 'mo(17,10**20-1) == 1499522760'; 
   $b = 100001;
   is mo(54, $b), 9090, 'mo(54,100001) == 9090';

}</lang>

Output:
ok 1 - 10 factors correctly
ok 2 - 21 factors correctly
ok 3 - 25 factors correctly
ok 4 - 150 factors correctly
ok 5 - 1231 factors correctly
ok 6 - 123141 factors correctly
ok 7 - 34131 factors correctly
ok 8 - mo(37,1000) == 100
ok 9 - mo(2,10**20-1) == 3748806900
ok 10 - mo(17,10**20-1) == 1499522760
ok 11 - mo(54,100001) == 9090

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>

Racket

The Racket function unit-group-order from racket/math computes the multiplicative order of an element a in Zn. An implementation of the algorithm in the tast description is shown below.

<lang racket>

  1. lang racket

(require math)

(define (order a n) 

 (unless (coprime? a n) (error 'order "arguments must be coprime"))
 (for/fold ([o 1]) ([r (factorize n)])
   (lcm o (order1 a r))))

(define (order1 a p&e) 

 (match-define (list p e) p&e)
 (define m (expt p e))
 (define t (* (- p 1) (expt p (- e 1))))
 (define qs
   (for/fold ([qs '(1)]) ([f (factorize t)])
      (match f [(list f0 f1)
                (for*/list ([q qs] [j (in-range (+ 1 f1))])
                  (* q (expt f0 j)))])))
 (for/or ([q (sort qs <)] #:when (= (modular-expt a q m) 1)) q))


(order 37 1000) (order (+ (expt 10 100) 1) 7919) (order (+ (expt 10 1000) 1) 15485863) (order (- (expt 10 10000) 1) 22801763489) (order 13 (+ 1 (expt 10 80))) </lang> Output: <lang racket> 100 3959 15485862 22801763488 109609547199756140150989321269669269476675495992554276140800 </lang>

Ruby

<lang ruby>require 'prime'

def powerMod(b, p, m) 

 p.to_s(2).each_char.inject(1) do |result, bit|
   result = (result * result) % m
   bit=='1' ? (result * b) % m : result
 end

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) do |result, f|
   result.lcm(multOrder_(a, *f))
 end

end

puts multOrder(37, 1000) b = 10**20-1 puts multOrder(2, b) puts multOrder(17,b) 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>

Output:
100
3748806900
1499522760
9090
1
Everything checks.

Seed7

<lang seed7>$ include "seed7_05.s7i"; 

 include "bigint.s7i";

const type: oneFactor is new struct 

   var bigInteger: prime is 0_;
   var integer: exp is 0;
 end struct;

const func oneFactor: oneFactor (in bigInteger: prime, in integer: exp) is func 

 result
   var oneFactor: aFactor is oneFactor.value;
 begin
   aFactor.prime := prime;
   aFactor.exp := exp;
 end func;

const func array oneFactor: factor (in var bigInteger: n) is func 

 result
   var array oneFactor: pf is 0 times oneFactor.value;
 local
   var integer: e is 0;
   var bigInteger: d is 0_;
   var bigInteger: s is 0_;
 begin
   e := lowestSetBit(n);
   if e > 0 then
     n >>:= e;
     pf := [] (oneFactor(2_, e));
   end if;
   s := sqrt(n);
   d := 3_;
   while n > 1_ do
     if d > s then
       d := n;
     end if;
     e := 0;
     while n rem d = 0_ do
       n := n div d;
       incr(e);
     end while;
     if e > 0 then
       pf &:= oneFactor(d, e);
       s := sqrt(n);
     end if;
     d +:= 2_;
   end while;
 end func;

const func bigInteger: moBachShallit58(in bigInteger: a, in bigInteger: n, in array oneFactor: pf) is func 

 result
   var bigInteger: mo is 0_;
 local
   var bigInteger: n1 is 0_;
   var oneFactor: pe is oneFactor.value;
   var bigInteger: x is 0_;
   var bigInteger: y is 0_;
   var integer: o is 0;
   var bigInteger: o1 is 0_;
 begin
   n1 := n - 1_;
   mo := 1_;
   for pe range pf do
     y := n1 div pe.prime ** pe.exp;
     x := modPow(a, y, n);
     o := 0;
     while x > 1_ do
       x := modPow(x, pe.prime, n);
       incr(o);
     end while;
     o1 := pe.prime ** o;
     mo *:= o1 div gcd(mo, o1);
   end for;
 end func;

const func boolean: isProbablyPrime (in bigInteger: primeCandidate, in var integer: count) is func 

 result
   var boolean: isProbablyPrime is TRUE;
 local
   var bigInteger: aRandomNumber is 0_;
 begin
   while isProbablyPrime and count > 0 do
     aRandomNumber := rand(1_, pred(primeCandidate));
     isProbablyPrime := modPow(aRandomNumber, pred(primeCandidate), primeCandidate) = 1_;
     decr(count);
   end while;
   # writeln(count);
 end func;

const proc: moTest (in bigInteger: a, in bigInteger: n) is func 

 begin
   if bitLength(a) < 100 then
     write("ord(" <& a <& ")");
   else
     write("ord([big])");
   end if;
   if bitLength(n) < 100 then
     write(" mod " <& n <& " ");
   else
     write(" mod [big] ");
   end if;
   if not isProbablyPrime(n, 20) then
     writeln("not computed.  modulus must be prime for this algorithm.")
   else
     writeln("= " <& moBachShallit58(a, n, factor(n - 1_)));
   end if;
 end func;

const proc: main is func 

 local
   var bigInteger: b is 100_;
 begin
   moTest(37_, 3343_);
   moTest(10_ ** 100 + 1_, 7919_);
   moTest(10_ ** 1000 + 1_, 15485863_);
   moTest(10_ ** 10000 - 1_, 22801763489_);
   moTest(1511678068_, 7379191741_);
   moTest(3047753288_, 2257683301_);
 end func;</lang>
Output:
ord(37) mod 3343 = 1114
ord([big]) mod 7919 = 3959
ord([big]) mod 15485863 = 15485862
ord([big]) mod 22801763489 = 22801763488
ord(1511678068) mod 7379191741 = 614932645
ord(3047753288) mod 2257683301 = 62713425

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>

zkl

Translation of: Python

Using Extensible prime generator#zkl and the GMP library for lcm (least common multiple), pow and powm ((n^e)%m)

It would probably be nice to memoize the prime numbers but that isn't necessary for this task. <lang zkl>var BN =Import("zklBigNum"); var Sieve=Import("sieve"); 

   // factor n into powers of primes
   // eg 9090 == 2^1 * 3^2 * 5^1 * 101^1

fcn factor2PP(n){ // lazy factors using lazy primes --> (prime,power) ... 

  Utils.Generator(fcn(a){
     primes:=Utils.Generator(Sieve.postponed_sieve);
     foreach p in (primes){

e:=0; while(a%p == 0){ a /= p; e+=1; } if (e) vm.yield(p,e); if (a<p*p) break; 

     }
     if (a>1) vm.yield(a,1);
  },n)

}

fcn _multOrdr1(a,p,e){ 

  m:=p.pow(e);
  t:=m/p*(p - 1);
  qs:=L(BN(1));
  foreach p2,e2 in (factor2PP(t)){ 
     qs=[[(e,q); [0..e2]; qs; '{ q*BN(p2).pow(e) }]];
  }
  qs.filter1('wrap(q){ a.powm(q,m)==1 });

}

fcn multiOrder(a,m){ 

  if (m.gcd(a)!=1) throw(Exception.ValueError("Not co-prime"));
  res:=BN(1);
  foreach p,e in (factor2PP(m)){
     res = res.lcm(_multOrdr1(BN(a),BN(p),e));
  }
  return(res);

}</lang> <lang zkl>multiOrder(37,1000).println(); b:=BN(10).pow(20)-1; multiOrder(2,b).println(); multiOrder(17,b).println();

b=0d10_0001; [BN(1)..multiOrder(54,b)-1].filter1('wrap(r,b54){b54.powm(r,b)==1},BN(54)) : if (_) println("Exists a power r < 9090 where (54^r)%b)==1"); else println("Everything checks.");</lang>

Output:
100
3748806900
1499522760
Everything checks.
Retrieved from "https://rosettacode.org/wiki/Multiplicative_order?oldid=48310"