Multiplicative order: Difference between revisions
Underscore (talk | contribs) m (→{{header|Haskell}}: Added lang tags.) |
m (Fixed lang tags.) |
||
Line 44: | Line 44: | ||
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}} |
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}} |
||
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput, also it generates a call to undefined C LONG externals }} --> |
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput, also it generates a call to undefined C LONG externals }} --> |
||
<lang |
<lang algol68>MODE LOOPINT = INT; |
||
MODE POWMODSTRUCT = LONG INT; |
MODE POWMODSTRUCT = LONG INT; |
||
Line 237: | Line 237: | ||
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 ''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. |
||
<lang j>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 |
The dyadic verb ''mopk'' expects a pair of prime and exponent |
||
Line 251: | Line 251: | ||
exponents ''q^d'' into a product. |
exponents ''q^d'' into a product. |
||
<lang j>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: |
For example: |
||
37 mo 1000 |
<lang j> 37 mo 1000 |
||
100 |
|||
2 mo _1+10^80x |
|||
190174169488577769580266953193403101748804183400400</lang> |
|||
=={{header|Mathematica}}== |
=={{header|Mathematica}}== |
||
In Mathematica this is really easy, as this function is built-in: |
In Mathematica this is really easy, as this function is built-in: |
||
<lang Mathematica> |
<lang Mathematica>MultiplicativeOrder[k,n] |
||
MultiplicativeOrder[k,n] |
|||
gives the multiplicative order of k modulo n, defined as the smallest integer m such that k^m == 1 mod 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,...}] |
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. |
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.</lang> |
||
⚫ | |||
Examples: |
Examples: |
||
<lang Mathematica> |
<lang Mathematica>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]</lang> |
|||
MultiplicativeOrder[11, 1 + 10^100] |
|||
⚫ | |||
gives back: |
gives back: |
||
<lang Mathematica> |
<lang Mathematica>100 |
||
100 |
|||
3959 |
3959 |
||
15485862 |
15485862 |
||
22801763488 |
22801763488 |
||
109609547199756140150989321269669269476675495992554276140800 |
109609547199756140150989321269669269476675495992554276140800 |
||
2583496112724752500580158969425549088007844580826869433740066152289289764829816356800 |
2583496112724752500580158969425549088007844580826869433740066152289289764829816356800</lang> |
||
</lang> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 12:14, 21 November 2009
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). For 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 wrt. 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(k(lg p)4/(lg lg p)) .
ALGOL 68
<lang algol68>MODE LOOPINT = INT;
MODE POWMODSTRUCT = LONG INT; PR READ "prelude/pow_mod.a68" PR;
MODE SORTSTRUCT = LONG INT; PR READ "prelude/sort.a68" PR;
PROC gcd = (LONG INT in a, in b)LONG INT: (
LONG INT a := in a, b := in b; WHILE b /= 0 DO LONG INT swap=a; a:=b; b:=swap %* a OD; a
);
PROC lcm = (LONG INT a, b)LONG INT : (a*b) OVER gcd(a, b);
MODE YIELDBOOL = PROC(BOOL)VOID; MODE FORBOOL = PROC(YIELDBOOL)VOID;
OP ANY = (FORBOOL sequence)BOOL: (
# FOR value IN sequence DO # sequence((BOOL value)VOID: IF value THEN return true FI ); else: FALSE EXIT return true: TRUE
);
OP ALL = (FORBOOL sequence)BOOL: (
# FOR value IN sequence DO # sequence((BOOL value)VOID: IF NOT value THEN return false FI ); else: TRUE EXIT return false: FALSE
);
PROC is prime = (LONG INT p)BOOL:
( p > 1 |#ANDF# ALL((YIELDBOOL yield)VOID: for 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 for 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 for factored = (LONG INT in a, PROC (LONG INT,LONG INT)VOID yield)VOID: (
LONG INT a := in a; # FOR p IN for primes DO # for primes( (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 ) ); 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 # for factored(t, (LONG INT f0,LONG INT 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 ) ); 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
);
MODE YIELDLONGINT = PROC(LONG INT)VOID; MODE FORLONGINT = PROC(YIELDLONGINT)VOID;
PROC for 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 #
iterator((LONG INT next)VOID: out := diadic(out, next) ); out
);
PROC mult order = (LONG INT a, LONG INT m)LONG INT: (
PROC for mofs = (YIELDLONGINT yield)VOID:( # FOR r IN factored(m) DO # for factored(m, (LONG INT p, LONG INT count)VOID: yield(mult0rdr1(a,p,count)) ) ); for reduce(lcm, for 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
)</lang> Output:
100 3748806900 1499522760 9090 1 Everything checks.
Haskell
Assuming a function
<lang haskell>primeFacsExp :: Integer -> [(Integer, Int)]</lang>
to calculate prime power factors, and another function
<lang haskell>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"</lang>
to efficiently calculate powers modulo some Integral
, we get
<lang haskell>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</lang>
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.
<lang j>mo=: 4 : 0
a=. x: x m=. x: y assert. 1=a+.m *./ a mopk"1 |: __ q: m
)</lang>
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.
<lang j>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
)</lang>
For example:
<lang j> 37 mo 1000 100
2 mo _1+10^80x
190174169488577769580266953193403101748804183400400</lang>
Mathematica
In Mathematica this is really easy, as this function is built-in: <lang Mathematica>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.</lang>
Examples: <lang Mathematica>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]</lang> gives back: <lang Mathematica>100 3959 15485862 22801763488 109609547199756140150989321269669269476675495992554276140800 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800</lang>
Python
<lang 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.'</lang>
Ruby
<lang ruby>require 'rational' # for lcm require 'mathn' # for prime_division
def powerMod(b, p, m)
result = 1 bits = p.to_s(2) for bit in bits.split() result = (result * result) % m if bit == '1' result = (result * b) % m end end result
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) {|result, f| result.lcm(multOrder_(a, *f)) }
end
puts multOrder(37, 1000) # 100 b = 10**20-1 puts multOrder(2, b) # 3748806900 puts multOrder(17,b) # 1499522760 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</lang>
Tcl
Slavishly
Uses struct::list
package from
<lang tcl>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"</lang