Multiplicative order

From Rosetta Code
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).


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 with regard 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.


Task

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

11l

Translation of: D
T PExp
   BigInt prime
   Int exp
   F (prime, exp)
      .prime = prime
      .exp = exp

F isqrt(self)
   V b = self
   L
      V a = b
      b = (self I/ a + a) I/ 2
      I b >= a
         R a

F factor(BigInt n)
   [PExp] pf
   V nn = n
   V b = 0
   L ((nn % 2) == 0)
      nn I/= 2
      b++

   I b > 0
      pf [+]= PExp(BigInt(2), b)

   V s = isqrt(nn)
   V d = BigInt(3)
   L nn > 1
      I d > s
         d = nn
      V e = 0
      L
         V (div, rem) = divmod(nn, d)
         I bits:length(rem) > 0
            L.break
         nn = div
         e++

      I e > 0
         pf [+]= PExp(d, e)
         s = isqrt(nn)

      d += 2

   R pf

F moBachShallit58(BigInt a, BigInt n; pf)
   V n1 = n - 1
   V mo = BigInt(1)
   L(pe) pf
      V y = n1 I/ pow(pe.prime, BigInt(pe.exp))
      V o = 0
      V x = pow(a, y, n)
      L x > 1
         x = pow(x, pe.prime, n)
         o++
      V o1 = pow(pe.prime, BigInt(o))
      o1 I/= gcd(mo, o1)
      mo *= o1
   R mo

F moTest(a, n)
   I bits:length(a) < 100
      print(‘ord(’a‘)’, end' ‘’)
   E
      print(‘ord([big])’, end' ‘’)
   print(‘ mod ’n‘ = ’moBachShallit58(a, n, factor(n - 1)))

moTest(37, 3343)

moTest(pow(BigInt(10), 100) + 1, 7919)
moTest(pow(BigInt(10), 1000) + 1, 15485863)
moTest(pow(BigInt(10), 10000) - 1, BigInt(22801763489))

moTest(1511678068, 7379191741)
moTest(BigInt(‘3047753288’), BigInt(‘2257683301’))
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

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

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;

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

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;

This is a sample program using the Multiplicative_Order package:

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;

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
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;
# 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
)

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.

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("The multiplicative order of %d related to %d is %lu \n", 37, 1000, multi_order(37, 1000));
	printf("The multiplicative order of %d related to %d is %lu \n", 54, 100001, multi_order(54, 100001));
	return 0;
}

C#

Translation of: Java
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Threading;

namespace MultiplicativeOrder {
    // Taken from https://stackoverflow.com/a/33918233
    public static class PrimeExtensions {
        // Random generator (thread safe)
        private static ThreadLocal<Random> s_Gen = new ThreadLocal<Random>(
          () => {
              return new Random();
          }
        );

        // Random generator (thread safe)
        private static Random Gen {
            get {
                return s_Gen.Value;
            }
        }

        public static bool IsProbablyPrime(this BigInteger value, int witnesses = 10) {
            if (value <= 1)
                return false;

            if (witnesses <= 0)
                witnesses = 10;

            BigInteger d = value - 1;
            int s = 0;

            while (d % 2 == 0) {
                d /= 2;
                s += 1;
            }

            byte[] bytes = new byte[value.ToByteArray().LongLength];
            BigInteger a;

            for (int i = 0; i < witnesses; i++) {
                do {
                    Gen.NextBytes(bytes);

                    a = new BigInteger(bytes);
                }
                while (a < 2 || a >= value - 2);

                BigInteger x = BigInteger.ModPow(a, d, value);
                if (x == 1 || x == value - 1)
                    continue;

                for (int r = 1; r < s; r++) {
                    x = BigInteger.ModPow(x, 2, value);

                    if (x == 1)
                        return false;
                    if (x == value - 1)
                        break;
                }

                if (x != value - 1)
                    return false;
            }

            return true;
        }
    }

    static class Helper {
        public static BigInteger Sqrt(this BigInteger self) {
            BigInteger b = self;
            while (true) {
                BigInteger a = b;
                b = self / a + a >> 1;
                if (b >= a) return a;
            }
        }

        public static long BitLength(this BigInteger self) {
            BigInteger bi = self;
            long bitlength = 0;
            while (bi != 0) {
                bitlength++;
                bi >>= 1;
            }
            return bitlength;
        }

        public static bool BitTest(this BigInteger self, int pos) {
            byte[] arr = self.ToByteArray();
            int idx = pos / 8;
            int mod = pos % 8;
            if (idx >= arr.Length) {
                return false;
            }
            return (arr[idx] & (1 << mod)) > 0;
        }
    }

    class PExp {
        public PExp(BigInteger prime, int exp) {
            Prime = prime;
            Exp = exp;
        }

        public BigInteger Prime { get; }

        public int Exp { get; }
    }

    class Program {
        static void MoTest(BigInteger a, BigInteger n) {
            if (!n.IsProbablyPrime(20)) {
                Console.WriteLine("Not computed. Modulus must be prime for this algorithm.");
                return;
            }
            if (a.BitLength() < 100) {
                Console.Write("ord({0})", a);
            } else {
                Console.Write("ord([big])");
            }
            if (n.BitLength() < 100) {
                Console.Write(" mod {0} ", n);
            } else {
                Console.Write(" mod [big] ");
            }
            BigInteger mob = MoBachShallit58(a, n, Factor(n - 1));
            Console.WriteLine("= {0}", mob);
        }

        static BigInteger MoBachShallit58(BigInteger a, BigInteger n, List<PExp> pf) {
            BigInteger n1 = n - 1;
            BigInteger mo = 1;
            foreach (PExp pe in pf) {
                BigInteger y = n1 / BigInteger.Pow(pe.Prime, pe.Exp);
                int o = 0;
                BigInteger x = BigInteger.ModPow(a, y, BigInteger.Abs(n));
                while (x > 1) {
                    x = BigInteger.ModPow(x, pe.Prime, BigInteger.Abs(n));
                    o++;
                }
                BigInteger o1 = BigInteger.Pow(pe.Prime, o);
                o1 = o1 / BigInteger.GreatestCommonDivisor(mo, o1);
                mo = mo * o1;
            }
            return mo;
        }

        static List<PExp> Factor(BigInteger n) {
            List<PExp> pf = new List<PExp>();
            BigInteger nn = n;
            int e = 0;
            while (!nn.BitTest(e)) e++;
            if (e > 0) {
                nn = nn >> e;
                pf.Add(new PExp(2, e));
            }
            BigInteger s = nn.Sqrt();
            BigInteger d = 3;
            while (nn > 1) {
                if (d > s) d = nn;
                e = 0;
                while (true) {
                    BigInteger div = BigInteger.DivRem(nn, d, out BigInteger rem);
                    if (rem.BitLength() > 0) break;
                    nn = div;
                    e++;
                }
                if (e > 0) {
                    pf.Add(new PExp(d, e));
                    s = nn.Sqrt();
                }
                d = d + 2;
            }

            return pf;
        }

        static void Main(string[] args) {
            MoTest(37, 3343);
            MoTest(BigInteger.Pow(10, 100) + 1, 7919);
            MoTest(BigInteger.Pow(10, 1000) + 1, 15485863);
            MoTest(BigInteger.Pow(10, 10000) - 1, 22801763489);
            MoTest(1511678068, 7379191741);
            MoTest(3047753288, 2257683301);
        }
    }
}
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

C++

Translation of: C
#include <algorithm>
#include <bitset>
#include <iostream>
#include <vector>

typedef unsigned long ulong;
std::vector<ulong> primes;

typedef struct {
    ulong p, e;
} prime_factor; /* prime, exponent */

void sieve() {
    /* 65536 = 2^16, so we can factor all 32 bit ints */
    constexpr int SIZE = 1 << 16;

    std::bitset<SIZE> bits;
    bits.flip(); // set all bits
    bits.reset(0);
    bits.reset(1);
    for (int i = 0; i < 256; i++) {
        if (bits.test(i)) {
            for (int j = i * i; j < SIZE; j += i) {
                bits.reset(j);
            }
        }
    }

    /* collect primes into a list. slightly faster this way if dealing with large numbers */
    for (int i = 0; i < SIZE; i++) {
        if (bits.test(i)) {
            primes.push_back(i);
        }
    }
}

auto get_prime_factors(ulong n) {
    std::vector<prime_factor> lst;
    ulong e, p;

    for (ulong i = 0; i < primes.size(); i++) {
        p = primes[i];
        if (p * p > n) break;
        for (e = 0; !(n % p); n /= p, e++);
        if (e) {
            lst.push_back({ p, e });
        }
    }

    if (n != 1) {
        lst.push_back({ n, 1 });
    }
    return lst;
}

auto get_factors(ulong n) {
    auto f = get_prime_factors(n);
    std::vector<ulong> lst{ 1 };

    size_t len2 = 1;
    /* L = (1); L = (L, L * p**(1 .. e)) forall((p, e)) */
    for (size_t i = 0; i < f.size(); i++, len2 = lst.size()) {
        for (ulong j = 0, p = f[i].p; j < f[i].e; j++, p *= f[i].p) {
            for (size_t k = 0; k < len2; k++) {
                lst.push_back(lst[k] * p);
            }
        }
    }

    std::sort(lst.begin(), lst.end());
    return lst;
}

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

ulong ipow(ulong a, ulong p) {
    ulong r = 1;
    while (p) {
        if (p & 1) 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 m = ipow(p, e);
    ulong t = m / p * (p - 1);
    auto fac = get_factors(t);
    for (size_t i = 0; i < fac.size(); i++) {
        if (mpow(a, fac[i], m) == 1) {
            return fac[i];
        }
    }
    return 0;
}

ulong multi_order(ulong a, ulong m) {
    auto pf = get_prime_factors(m);
    ulong res = 1;
    for (size_t i = 0; i < pf.size(); 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));   // expect 100
    printf("%lu\n", multi_order(54, 100001)); // expect 9090

    return 0;
}
Output:
100
9090

Clojure

Translation of Julie, then revised to be more clojure idiomatic. It references some external modules for factoring and integer exponentiation.

(defn gcd [a b]
   (if (zero? b)
      a
      (recur b (mod a b))))

(defn lcm [a b]
   (/ (* a b) (gcd a b)))

(def NaN  (Math/log -1))

(defn ord' [a [p e]]
   (let [m   (imath/expt p e)
         t   (* (quot m p) (dec p))]
            (loop [dv (factor/divisors t)]
               (let [d (first dv)]
                  (if (= (mmath/expm a d m) 1)
                     d
                     (recur (next dv)))))))

(defn ord [a n]
   (if (not= (gcd a n) 1)
      NaN
      (->>
         (factor/factorize n)
         (map (partial ord' a))
         (reduce lcm))))
Output:
user=> (ord 37 1000)
100

D

Translation of: Java
import std.bigint;
import std.random;
import std.stdio;

struct PExp {
    BigInt prime;
    int exp;
}

BigInt gcd(BigInt x, BigInt y) {
    if (y == 0) {
        return x;
    }
    return gcd(y, x % y);
}

/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
BigInt modPow(BigInt b, BigInt e, BigInt n) {
    if (n == 1) return BigInt(0);
    BigInt result = 1;
    b = b % n;
    while (e > 0) {
        if (e % 2 == 1) {
            result = (result * b) % n;
        }
        e >>= 1;
        b = (b*b) % n;
    }
    return result;
}

BigInt pow(long b, long e) {
    return pow(BigInt(b), BigInt(e));
}
BigInt pow(BigInt b, BigInt e) {
    if (e == 0) {
        return BigInt(1);
    }

    BigInt result = 1;
    while (e > 1) {
        if (e % 2 == 0) {
            b *= b;
            e /= 2;
        } else {
            result *= b;
            b *= b;
            e = (e - 1) / 2;
        }
    }

    return b * result;
}

BigInt sqrt(BigInt self) {
    BigInt b = self;
    while (true) {
        BigInt a = b;
        b = self / a + a >> 1;
        if (b >= a) return a;
    }
}

long bitLength(BigInt self) {
    BigInt bi = self;
    long length;
    while (bi != 0) {
        length++;
        bi >>= 1;
    }
    return length;
}

PExp[] factor(BigInt n) {
    PExp[] pf;
    BigInt nn = n;
    int b = 0;
    int e = 1;
    while ((nn & e) == 0) {
        e <<= 1;
        b++;
    }
    if (b > 0) {
        nn = nn >> b;
        pf ~= PExp(BigInt(2), b);
    }
    BigInt s = nn.sqrt();
    BigInt d = 3;
    while (nn > 1) {
        if (d > s) d = nn;
        e = 0;
        while (true) {
            BigInt div, rem;
            nn.divMod(d, div, rem);
            if (rem.bitLength > 0) break;
            nn = div;
            e++;
        }
        if (e > 0) {
            pf ~= PExp(d, e);
            s = nn.sqrt();
        }
        d += 2;
    }

    return pf;
}

BigInt moBachShallit58(BigInt a, BigInt n, PExp[] pf) {
    BigInt n1 = n - 1;
    BigInt mo = 1;
    foreach(pe; pf) {
        BigInt y = n1 / pe.prime.pow(BigInt(pe.exp));
        int o = 0;
        BigInt x = a.modPow(y, n);
        while (x > 1) {
            x = x.modPow(pe.prime, n);
            o++;
        }
        BigInt o1 = pe.prime.pow(BigInt(o));
        o1 = o1 / gcd(mo, o1);
        mo = mo * o1;
    }
    return mo;
}

void moTest(ulong a, ulong n) {
    moTest(BigInt(a), n);
}
void moTest(BigInt a, ulong n) {
    // Commented out because the implementations tried all failed for the -2 and -3 tests.
    // if (!n.isProbablePrime()) {
        // writeln("Not computed. Modulus must be prime for this algorithm.");
        // return;
    // }
    if (a.bitLength < 100) {
        write("ord(", a, ")");
    } else {
        write("ord([big])");
    }
    write(" mod ", n, " ");
    BigInt nn = n;
    BigInt mob = moBachShallit58(a, nn, factor(nn - 1));
    writeln("= ", mob);
}

void main() {
    moTest(37, 3343);

    moTest(pow(10, 100) + 1, 7919);
    moTest(pow(10, 1000) + 1, 15485863);
    moTest(pow(10, 10000) - 1, 22801763489);

    moTest(1511678068, 7379191741);
    moTest(3047753288, 2257683301);
}
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

EchoLisp

(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

Factor

Works with: Factor version 0.99 2020-01-23
USING: kernel math math.functions math.primes.factors sequences ;

: (ord) ( a pair -- n )
    first2 dupd ^ swap dupd [ /i ] keep 1 - * divisors
    [ swap ^mod 1 = ] 2with find nip ;
    

: ord ( a n -- m )
    2dup gcd nip 1 =
    [ group-factors [ (ord) ] with [ lcm ] map-reduce ]
    [ 2drop 0/0. ] if ;
Output:
IN: scratchpad 37 1000 ord .
100
IN: scratchpad 10 100 ^ 1 + 7919 ord .
3959

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 = []pExp{{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)
}
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

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

import Data.List (foldl1') --'

foldl1_ = foldl1' --'

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

Java

Translation of: Kotlin
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class MultiplicativeOrder {
    private static final BigInteger ONE = BigInteger.ONE;
    private static final BigInteger TWO = BigInteger.valueOf(2);
    private static final BigInteger THREE = BigInteger.valueOf(3);
    private static final BigInteger TEN = BigInteger.TEN;

    private static class PExp {
        BigInteger prime;
        long exp;

        PExp(BigInteger prime, long exp) {
            this.prime = prime;
            this.exp = exp;
        }
    }

    private static void moTest(BigInteger a, BigInteger n) {
        if (!n.isProbablePrime(20)) {
            System.out.println("Not computed. Modulus must be prime for this algorithm.");
            return;
        }
        if (a.bitLength() < 100) System.out.printf("ord(%s)", a);
        else System.out.print("ord([big])");
        if (n.bitLength() < 100) System.out.printf(" mod %s ", n);
        else System.out.print(" mod [big] ");
        BigInteger mob = moBachShallit58(a, n, factor(n.subtract(ONE)));
        System.out.println("= " + mob);
    }

    private static BigInteger moBachShallit58(BigInteger a, BigInteger n, List<PExp> pf) {
        BigInteger n1 = n.subtract(ONE);
        BigInteger mo = ONE;
        for (PExp pe : pf) {
            BigInteger y = n1.divide(pe.prime.pow((int) pe.exp));
            long o = 0;
            BigInteger x = a.modPow(y, n.abs());
            while (x.compareTo(ONE) > 0) {
                x = x.modPow(pe.prime, n.abs());
                o++;
            }
            BigInteger o1 = BigInteger.valueOf(o);
            o1 = pe.prime.pow(o1.intValue());
            o1 = o1.divide(mo.gcd(o1));
            mo = mo.multiply(o1);
        }
        return mo;
    }

    private static List<PExp> factor(BigInteger n) {
        List<PExp> pf = new ArrayList<>();
        BigInteger nn = n;
        Long e = 0L;
        while (!nn.testBit(e.intValue())) e++;
        if (e > 0L) {
            nn = nn.shiftRight(e.intValue());
            pf.add(new PExp(TWO, e));
        }
        BigInteger s = sqrt(nn);
        BigInteger d = THREE;
        while (nn.compareTo(ONE) > 0) {
            if (d.compareTo(s) > 0) d = nn;
            e = 0L;
            while (true) {
                BigInteger[] qr = nn.divideAndRemainder(d);
                if (qr[1].bitLength() > 0) break;
                nn = qr[0];
                e++;
            }
            if (e > 0L) {
                pf.add(new PExp(d, e));
                s = sqrt(nn);
            }
            d = d.add(TWO);
        }
        return pf;
    }

    private static BigInteger sqrt(BigInteger n) {
        BigInteger b = n;
        while (true) {
            BigInteger a = b;
            b = n.divide(a).add(a).shiftRight(1);
            if (b.compareTo(a) >= 0) return a;
        }
    }

    public static void main(String[] args) {
        moTest(BigInteger.valueOf(37), BigInteger.valueOf(3343));

        BigInteger b = TEN.pow(100).add(ONE);
        moTest(b, BigInteger.valueOf(7919));

        b = TEN.pow(1000).add(ONE);
        moTest(b, BigInteger.valueOf(15485863));

        b = TEN.pow(10000).subtract(ONE);
        moTest(b, BigInteger.valueOf(22801763489L));

        moTest(BigInteger.valueOf(1511678068), BigInteger.valueOf(7379191741L));
        moTest(BigInteger.valueOf(3047753288L), BigInteger.valueOf(2257683301L));
    }
}
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

Julia

(Uses the factors function from Factors of an integer#Julia.)

using Primes

function factors(n)
    f = [one(n)]
    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    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

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

Kotlin

Translation of: Go
// version 1.2.10

import java.math.BigInteger

val bigOne   = BigInteger.ONE
val bigTwo   = 2.toBigInteger()
val bigThree = 3.toBigInteger()
val bigTen   = BigInteger.TEN

class PExp(val prime: BigInteger, val exp: Long)

fun moTest(a: BigInteger, n: BigInteger) {
    if (!n.isProbablePrime(20)) {
        println("Not computed. Modulus must be prime for this algorithm.")
        return
    }
    if (a.bitLength() < 100) print("ord($a)") else print("ord([big])")
    if (n.bitLength() < 100) print(" mod $n ") else print(" mod [big] ")
    val mob = moBachShallit58(a, n, factor(n - bigOne))
    println("= $mob")
}

fun moBachShallit58(a: BigInteger, n: BigInteger, pf: List<PExp>): BigInteger {
    val n1 = n - bigOne
    var mo = bigOne
    for (pe in pf) {
        val y = n1 / pe.prime.pow(pe.exp.toInt())
        var o = 0L
        var x = a.modPow(y, n.abs())
        while (x > bigOne) {
            x = x.modPow(pe.prime, n.abs())
            o++
        }
        var o1 = o.toBigInteger()
        o1 = pe.prime.pow(o1.toInt())
        o1 /= mo.gcd(o1)
        mo *= o1
    }
    return mo
}

fun factor(n: BigInteger): List<PExp> {
    val pf = mutableListOf<PExp>()
    var nn = n
    var e = 0L
    while (!nn.testBit(e.toInt())) e++
    if (e > 0L) {
        nn = nn shr e.toInt()
        pf.add(PExp(bigTwo, e))
    }
    var s = bigSqrt(nn)
    var d = bigThree
    while (nn > bigOne) {
        if (d > s) d = nn
        e = 0L
        while (true) {
            val (q, r) = nn.divideAndRemainder(d)
            if (r.bitLength() > 0) break
            nn = q
            e++
        }
        if (e > 0L) {
            pf.add(PExp(d, e))
            s = bigSqrt(nn)
        }
        d += bigTwo
    }
    return pf
}

fun bigSqrt(n: BigInteger): BigInteger {
    var b = n
    while (true) {
        val a = b
        b = (n / a + a) shr 1
        if (b >= a) return a
    }
}

fun main(args: Array<String>) {
    moTest(37.toBigInteger(), 3343.toBigInteger())

    var b = bigTen.pow(100) + bigOne
    moTest(b, 7919.toBigInteger())

    b = bigTen.pow(1000) + bigOne
    moTest(b, BigInteger("15485863"))

    b = bigTen.pow(10000) - bigOne
    moTest(b, BigInteger("22801763489"))

    moTest(BigInteger("1511678068"), BigInteger("7379191741"))
    moTest(BigInteger("3047753288"), BigInteger("2257683301"))
}
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

Maple

numtheory:-order( a, n )

For example,

> numtheory:-order( 37, 1000 );
                       100

Mathematica/Wolfram Language

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:

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]

gives back:

100
3959
15485862
22801763488
109609547199756140150989321269669269476675495992554276140800
2583496112724752500580158969425549088007844580826869433740066152289289764829816356800

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 */

Nim

Translation of: Kotlin
Library: bignum
import strformat
import bignum

type PExp = tuple[prime: Int; exp: uint]

let
  one = newInt(1)
  two = newInt(2)
  ten = newInt(10)


func sqrt(n: Int): Int =
  var s = n
  while true:
    result = s
    s = (n div result + result) shr 1
    if s >= result: break


proc factor(n: Int): seq[PExp] =
  var n = n
  var e = 0u
  while n.bit(e) == 0: inc e
  if e != 0:
    n = n shr e
    result.add (two, e)
  var s = sqrt(n)
  var d = newInt(3)
  while n > one:
    if d > s: d = n
    e = 0u
    while true:
      let (q, r) = divMod(n, d)
      if not r.isZero: break
      n = q
      inc e
    if e != 0:
      result.add (d.clone, e)
      s = sqrt(n)
    inc d, two


proc moBachShallit58(a, n: Int; pf: seq[PExp]): Int =
  let n = abs(n)
  let n1 = n - one
  result = newInt(1)
  for pe in pf:
    let y = n1 div pe.prime.pow(pe.exp)
    var o = 0u
    var x = a.exp(y.toInt.uint, n)
    while x > one:
      x = x.exp(pe.prime.toInt.uint, n)
      inc o
    var o1 = pe.prime.pow(o)
    o1 = o1 div gcd(result, o1)
    result *= o1


proc moTest(a, n: Int) =
  if n.probablyPrime(25) == 0:
    echo "Not computed. Modulus must be prime for this algorithm."
    return

  stdout.write if a.bitLen < 100: &"ord({a})" else: "ord([big])"
  stdout.write if n.bitlen < 100: &" mod {n}" else: " mod [big]"
  let mob = moBachShallit58(a, n, factor(n - one))
  echo &" = {mob}"


when isMainModule:
  moTest(newInt(37), newInt(3343))

  var b = ten.pow(100) + one
  motest(b, newInt(7919))

  b = ten.pow(1000) + one
  moTest(b, newInt("15485863"))

  b = ten.pow(10000) - one
  moTest(b, newInt("22801763489"))

  moTest(newInt("1511678068"), newInt("7379191741"))

  moTest(newInt("3047753288"), newInt("2257683301"))
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

PARI/GP

znorder(Mod(a,n))

Perl

Using modules:

Library: ntheory
use ntheory qw/znorder/;
say znorder(54, 100001);
use bigint; say znorder(11, 1 + 10**100);

or

use Math::Pari qw/znorder Mod/;
say znorder(Mod(54, 100001));
say znorder(Mod(11, 1 + Math::Pari::PARI(10)**100));

Phix

Translation of: Ruby
Library: Phix/mpfr
with javascript_semantics
include mpfr.e
 
procedure multi_order(mpz res, a, sequence p_and_k)
    mpz {pk,t,x,q,pz} = mpz_inits(5)
    mpz_set_si(res,1)
    if length(p_and_k)=1 then
        string {ps} = p_and_k
        mpz_set_str(pk,ps)
        mpz_sub_ui(t,pk,1)
    else
        atom {p, k} = p_and_k
        mpz_set_d(pz,p)
        mpz_pow_ui(pk,pz,k)
        mpz_pow_ui(t,pz,k-1)
        mpz_sub_ui(pz,pz,1)
        mpz_mul(t,t,pz)
    end if
    sequence pf = mpz_prime_factors(t)
    for i=1 to length(pf) do
        if length(pf[i])=1 then
            string {fs} = pf[i]
            mpz_set_str(q,fs)
            mpz_set(x,q)
        else
            {integer qi, integer ei} = pf[i]
            mpz_set_si(q,qi)
            mpz_pow_ui(x,q,ei)
        end if
        mpz_fdiv_q(x, t, x)
        mpz_powm(x,a,x,pk)
        integer guard = 0
        while mpz_cmp_si(x,1)!=0 do
            mpz_mul(res,res,q)
            mpz_powm(x,x,q,pk)
            guard += 1
            if guard>100 then ?9/0 end if -- (increase if rqd)
        end while
    end for
    x = mpz_free(x)
end procedure
 
function multiplicative_order(mpz a, m)
    mpz res = mpz_init(1),
        ri = mpz_init()
    mpz_gcd(ri,a,m)
    if mpz_cmp_si(ri,1)!=0 then return "(a,m) not coprime" end if
    sequence pf = mpz_prime_factors(m,10000) -- (increase if rqd)
    for i=1 to length(pf) do
        multi_order(ri,a,pf[i])
        mpz_lcm(res,res,ri)
    end for
    return mpz_get_str(res)
end function
 
function shorta(mpz n)
    string res = mpz_get_str(n)
    integer lr = length(res)
    if lr>80 then
        res[6..-6] = "..."
        res &= sprintf(" (%d digits)",lr)
    end if
    return res
end function
 
procedure mo_test(mpz a, n)
    string res = multiplicative_order(a, n)
    printf(1,"ord(%s) mod %s = %s\n",{shorta(a),shorta(n),res})
end procedure
 
function i(atom i) return mpz_init(i) end function -- (ugh)
function p10(integer e,i)   -- init to 10^e+i
    mpz res = mpz_init()
    mpz_ui_pow_ui(res,10,e)
    mpz_add_si(res,res,i)
    return res
end function
 
atom t = time()
mo_test(i(3), i(10))
mo_test(i(37), i(1000))
mo_test(i(37), i(10000))
mo_test(i(37), i(3343))
mo_test(i(37), i(3344))
mo_test(i(2), i(1000))
mo_test(p10(100,+1), i(7919))
mo_test(p10(1000,+1), i(15485863))
mo_test(p10(10000,-1), i(22801763489))
mo_test(i(1511678068), i(7379191741))
mo_test(i(3047753288), i(2257683301))
?"==="
mpz b = p10(20,-1)
mo_test(i(2), b)
mo_test(i(17),b)
mo_test(i(54),i(100001))
string s9090 = multiplicative_order(mpz_init(54),mpz_init(100001))
if s9090!="9090" then ?9/0 end if
mpz m54 = mpz_init(54),
    m100001 = mpz_init(100001)
mpz_powm_ui(b,m54,9090,m100001)
printf(1,"%s\n",mpz_get_str(b))
bool error = false
for r=1 to 9090-1 do
    mpz_powm_ui(b,m54,r,m100001)
    if mpz_cmp_si(b,1)=0 then
        printf(1,"mpz_powm_ui(54,%d,100001) gives 1!\n",r)
        error = true
        exit
    end if
end for
if not error then
    printf(1,"Everything checks. (%s)\n",{elapsed(time()-t)})
end if
Output:
ord(3) mod 10 = 4
ord(37) mod 1000 = 100
ord(37) mod 10000 = 500
ord(37) mod 3343 = 1114
ord(37) mod 3344 = 20
ord(2) mod 1000 = (a,m) not coprime
ord(10000...00001 (101 digits)) mod 7919 = 3959
ord(10000...00001 (1001 digits)) mod 15485863 = 15485862
ord(99999...99999 (10000 digits)) mod 22801763489 = 22801763488
ord(1511678068) mod 7379191741 = 614932645
ord(3047753288) mod 2257683301 = 62713425
"==="
ord(2) mod 99999999999999999999 = 3748806900
ord(17) mod 99999999999999999999 = 1499522760
ord(54) mod 100001 = 9090
1
Everything checks. (0.2s)

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

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

Output:

100
3959
15485862
22801763488
109609547199756140150989321269669269476675495992554276140800

Raku

(formerly Perl 6)

use Prime::Factor;

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 prime-factors($t).Bag -> $f {
        @qs = flat @qs.map(-> $q { (0..$f.value).map(-> $j { $q * $f.key ** $j }) });
    }

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

sub mo($a, $m) {
    $a gcd $m == 1 or die "$a and $m are not relatively prime";
    [lcm] flat 1, prime-factors($m).Bag.map: { mo-prime($a, .key, .value) };
}

multi MAIN('test') {
    use Test;

    for (10, 21, 25, 150, 1231, 123141, 34131) -> $n {
        is ([*] prime-factors($n).Bag.map( { .key ** .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';
}
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

REXX

/*REXX pgm computes multiplicative order of a minimum integer  N  such that  a^n mod m≡1*/
wa= 0;  wm= 0     /*       ═a═   ══m══     */    /*maximum widths of the A and M values.*/
@.=.;                @.1=   3      10
                     @.2=  37    1000
                     @.3=  37   10000
                     @.4=  37    3343
                     @.5=  37    3344
                     @.6=   2    1000
pad= left('', 9)
d= 500                                           /*use 500 decimal digits for a starter.*/
     do w=1  for 2                               /*when W≡1, find max widths of A and M.*/
       do j=1  while @.j\==.;         parse var  @.j     a  .  1  r  m  ,  n
       if w==1  then do;  wa= max(wa, length(a) );     wm= max(wm, length(m) );    iterate
                     end
       if m//a==0  then n= ' [solution not possible]'     /*test co─prime for  A and B. */
       numeric digits d                          /*start with  100  decimal digits.     */
       if n==''  then do n= 2;    p= r * a       /*compute product──may have an exponent*/
                      parse  var  p  'E'  _      /*try to extract the exponent from  P. */
                      if _\==''   then do;  numeric digits _+d  /*bump the decimal digs.*/
                                            p=r*a               /*recalculate integer P.*/
                                       end
                      if p//m==1  then leave     /*now, perform the nitty─gritty modulo.*/
                      r= p                       /*assign product to R for next multiply*/
                      end   /*n*/                /* [↑]    //   is really  ÷  remainder.*/
       say pad  'a='  right(a,wa)  pad  "m=" right(m,wm)  pad  'multiplicative order:'   n
       end   /*j*/
     end     /*w*/                               /*stick a fork in it,  we're all done. */
output:
          a=  3           m=    10           multiplicative order: 4
          a= 37           m=  1000           multiplicative order: 100
          a= 37           m= 10000           multiplicative order: 500
          a= 37           m=  3343           multiplicative order: 1114
          a= 37           m=  3344           multiplicative order: 20
          a=  2           m=  1000           multiplicative order:  [solution not possible]

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
Output:
100
3748806900
1499522760
9090
1
Everything checks.

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;
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

Sidef

Built-in:

say 37.znorder(1000)     #=> 100
say 54.znorder(100001)   #=> 9090
Translation of: Raku
func mo_prime(a, p, e) {
    var m  = p**e
    var t  = (p-1)*(p**(e-1))
    var qs = [1]

    for f in (t.factor_exp) {
        qs.map! {|q|
            0..f[1] -> map {|j| q * f[0]**j }...
        }
    }

    qs.sort.first_by {|q| powmod(a, q, m) == 1 }
}

func mo(a, m) {
    gcd(a, m) == 1 || die "#{a} and #{m} are not relatively prime"
    Math.lcm(1, m.factor_exp.map {|r| mo_prime(a, r...) }...)
}

say mo(37, 1000)
say mo(54, 100001)

with (10**20 - 1) {|b|
    say mo(2, b)
    say mo(17, b)
}
Output:
100
9090
3748806900
1499522760

Tcl

Translation of: Python
Library: Tcllib (Package: struct::list)
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    
}

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

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

# python's pow()
proc ::tcl::mathfunc::pypow {x y {z ""}} {
    expr {$z eq "" ? $x ** $y : ($x ** $y) % $z}
}

####################################################
# prime number generator
# 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"

Visual Basic .NET

Translation of: C#
Imports System.Numerics
Imports System.Runtime.CompilerServices
Imports System.Threading

Module Module1
    Private s_gen As New ThreadLocal(Of Random)(Function() New Random())

    Private Function Gen()
        Return s_gen.Value
    End Function

    <Extension()>
    Public Function IsProbablyPrime(value As BigInteger, Optional witnesses As Integer = 10) As Boolean
        If value <= 1 Then
            Return False
        End If

        If witnesses <= 0 Then
            witnesses = 10
        End If

        Dim d = value - 1
        Dim s = 0

        While d Mod 2 = 0
            d /= 2
            s += 1
        End While

        Dim bytes(value.ToByteArray.LongLength - 1) As Byte
        Dim a As BigInteger

        For i = 1 To witnesses
            Do
                Gen.NextBytes(bytes)

                a = New BigInteger(bytes)
            Loop While a < 2 OrElse a >= value - 2

            Dim x = BigInteger.ModPow(a, d, value)
            If x = 1 OrElse x = value - 1 Then
                Continue For
            End If

            For r = 1 To s - 1
                x = BigInteger.ModPow(x, 2, value)

                If x = 1 Then
                    Return False
                End If
                If x = value - 1 Then
                    Exit For
                End If
            Next

            If x <> value - 1 Then
                Return False
            End If
        Next

        Return True
    End Function

    <Extension()>
    Function Sqrt(self As BigInteger) As BigInteger
        Dim b = self
        While True
            Dim a = b
            b = self / a + a >> 1
            If b >= a Then
                Return a
            End If
        End While
        Throw New Exception("Should not have happened")
    End Function

    <Extension()>
    Function BitLength(self As BigInteger) As Long
        Dim bi = self
        Dim len = 0L
        While bi <> 0
            len += 1
            bi >>= 1
        End While
        Return len
    End Function

    <Extension()>
    Function BitTest(self As BigInteger, pos As Integer) As Boolean
        Dim arr = self.ToByteArray
        Dim i = pos \ 8
        Dim m = pos Mod 8
        If i >= arr.Length Then
            Return False
        End If
        Return (arr(i) And (1 << m)) > 0
    End Function

    Class PExp
        Sub New(p As BigInteger, e As Integer)
            Prime = p
            Exp = e
        End Sub

        Public ReadOnly Property Prime As BigInteger
        Public ReadOnly Property Exp As Integer
    End Class

    Function MoBachShallit58(a As BigInteger, n As BigInteger, pf As List(Of PExp)) As BigInteger
        Dim n1 = n - 1
        Dim mo As BigInteger = 1
        For Each pe In pf
            Dim y = n1 / BigInteger.Pow(pe.Prime, pe.Exp)
            Dim o = 0
            Dim x = BigInteger.ModPow(a, y, BigInteger.Abs(n))
            While x > 1
                x = BigInteger.ModPow(x, pe.Prime, BigInteger.Abs(n))
                o += 1
            End While
            Dim o1 = BigInteger.Pow(pe.Prime, o)
            o1 /= BigInteger.GreatestCommonDivisor(mo, o1)
            mo *= o1
        Next
        Return mo
    End Function

    Function Factor(n As BigInteger) As List(Of PExp)
        Dim pf As New List(Of PExp)
        Dim nn = n
        Dim e = 0
        While Not nn.BitTest(e)
            e += 1
        End While
        If e > 0 Then
            nn >>= e
            pf.Add(New PExp(2, e))
        End If
        Dim s = nn.Sqrt
        Dim d As BigInteger = 3
        While nn > 1
            If d > s Then
                d = nn
            End If
            e = 0
            While True
                Dim remainder As New BigInteger
                Dim div = BigInteger.DivRem(nn, d, remainder)
                If remainder.BitLength > 0 Then
                    Exit While
                End If
                nn = div
                e += 1
            End While
            If e > 0 Then
                pf.Add(New PExp(d, e))
                s = nn.Sqrt
            End If
            d += 2
        End While
        Return pf
    End Function

    Sub MoTest(a As BigInteger, n As BigInteger)
        If Not n.IsProbablyPrime(20) Then
            Console.WriteLine("Not computed. Modulus must be prime for this algorithm.")
            Return
        End If
        If a.BitLength < 100 Then
            Console.Write("ord({0})", a)
        Else
            Console.Write("ord([big])")
        End If
        If n.BitLength < 100 Then
            Console.Write(" mod {0}", n)
        Else
            Console.Write(" mod [big]")
        End If
        Dim mob = MoBachShallit58(a, n, Factor(n - 1))
        Console.WriteLine(" = {0}", mob)
    End Sub

    Sub Main()
        MoTest(37, 3343)
        MoTest(BigInteger.Pow(10, 100) + 1, 7919)
        MoTest(BigInteger.Pow(10, 1000) + 1, 15485863)
        MoTest(BigInteger.Pow(10, 10000) - 1, 22801763489)
        MoTest(1511678068, 7379191741)
        MoTest(3047753288, 2257683301)
    End Sub

End Module
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

Wren

Translation of: Kotlin
Library: Wren-big
import "./big" for BigInt

class PExp {
    construct new(prime, exp) {
        _prime = prime
        _exp = exp
    }
    prime { _prime }
    exp   { _exp }
}

var moBachShallit58 = Fn.new { |a, n, pf|
    var n1 = n - BigInt.one
    var mo = BigInt.one
    for (pe in pf) {
        var y = n1 / pe.prime.pow(pe.exp)
        var o = 0
        var x = a.modPow(y, n.abs)
        while (x > BigInt.one) {
            x = x.modPow(pe.prime, n.abs)
            o = o + 1
        }
        var o1 = BigInt.new(o)
        o1 = pe.prime.pow(o1)
        o1 = o1 / BigInt.gcd(mo, o1)
        mo = mo * o1
    }
    return mo
}

var factor = Fn.new { |n|
    var pf = []
    var nn = n.copy()
    var e = 0
    while (!nn.testBit(e)) e = e + 1
    if (e > 0) {
        nn = nn >> e
        pf.add(PExp.new(BigInt.two, e))
    }
    var s = nn.isqrt
    var d = BigInt.three
    while (nn > BigInt.one) {
        if (d > s) d = nn
        e = 0
        while (true) {
            var dm = nn.divMod(d)
            if (dm[1].bitLength > 0) break
            nn = dm[0]
            e = e + 1
        }
        if (e > 0) {
            pf.add(PExp.new(d, e))
            s = nn.isqrt
        }
        d = d + BigInt.two
    }
    return pf
}

var moTest = Fn.new { |a, n|
    if (!n.isProbablePrime(10)) {
        System.print("Not computed. Modulus must be prime for this algorithm.")
        return
    }
    System.write((a.bitLength < 100) ? "ord(%(a))"  : "ord([big])")
    System.write((n.bitLength < 100) ? " mod %(n) " : "mod([big])")
    var mob = moBachShallit58.call(a, n, factor.call(n - BigInt.one))
    System.print("= %(mob)")
}

moTest.call(BigInt.new(37), BigInt.new(3343))

var b = BigInt.ten.pow(100) + BigInt.one
moTest.call(b, BigInt.new(7919))

b = BigInt.ten.pow(1000) + BigInt.one
moTest.call(b, BigInt.new(15485863))

b = BigInt.ten.pow(10000) - BigInt.one
moTest.call(b, BigInt.new(22801763489))

moTest.call(BigInt.new(1511678068), BigInt.new(7379191741))
moTest.call(BigInt.new(3047753288), BigInt.new(2257683301))
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

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.

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);
}
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.");
Output:
100
3748806900
1499522760
Everything checks.