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




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
