user
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)).
assistant
(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
```