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




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"