Problem Statement

The multiplicative order of a relative to m is the least positive integer n such that a^n is 1 (modulo m).

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.

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

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"). Here is the implementation ("multiplicative_order.adb"): This is a sample program using the Multiplicative_Order package: The output from the sample program: Output: 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. Translation of Julie, then revised to be more clojure idiomatic. It references some external modules for factoring and integer exponentiation. Assuming a function to efficiently calculate powers modulo some Integral, we get 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. 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. For example: (Uses the factors function from Factors of an integer#Julia.) Example output (using big to employ arbitrary-precision arithmetic where needed): For example, 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: gives back: Using modules: or 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. Output: (formerly Perl 6) Built-in: 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.

Let's start with the solution:

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