Multiplicative order: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Sidef}}: Fix link: Perl 6 --> Raku) |
Alextretyak (talk | contribs) m (→{{header|11l}}: named tuples) |
||
(15 intermediate revisions by 9 users not shown) | |||
Line 45: | Line 45: | ||
The total cost is dominated by<tt> O(k(lg p)<sup>3</sup>)</tt> ,<tt> </tt>which is<tt> O((lg p)<sup>4</sup>/(lg lg p))</tt>. |
The total cost is dominated by<tt> O(k(lg p)<sup>3</sup>)</tt> ,<tt> </tt>which is<tt> O((lg p)<sup>4</sup>/(lg lg p))</tt>. |
||
<br><br> |
<br><br> |
||
=={{header|11l}}== |
|||
{{trans|D}} |
|||
<syntaxhighlight lang="11l">T PExp = (BigInt prime, Int exp) |
|||
F isqrt(self) |
|||
V b = self |
|||
L |
|||
V a = b |
|||
b = (self I/ a + a) I/ 2 |
|||
I b >= a |
|||
R a |
|||
F factor(BigInt n) |
|||
[PExp] pf |
|||
V nn = n |
|||
V b = 0 |
|||
L ((nn % 2) == 0) |
|||
nn I/= 2 |
|||
b++ |
|||
I b > 0 |
|||
pf [+]= PExp(BigInt(2), b) |
|||
V s = isqrt(nn) |
|||
V d = BigInt(3) |
|||
L nn > 1 |
|||
I d > s |
|||
d = nn |
|||
V e = 0 |
|||
L |
|||
V (div, rem) = divmod(nn, d) |
|||
I bits:length(rem) > 0 |
|||
L.break |
|||
nn = div |
|||
e++ |
|||
I e > 0 |
|||
pf [+]= PExp(d, e) |
|||
s = isqrt(nn) |
|||
d += 2 |
|||
R pf |
|||
F moBachShallit58(BigInt a, BigInt n; pf) |
|||
V n1 = n - 1 |
|||
V mo = BigInt(1) |
|||
L(pe) pf |
|||
V y = n1 I/ pow(pe.prime, BigInt(pe.exp)) |
|||
V o = 0 |
|||
V x = pow(a, y, n) |
|||
L x > 1 |
|||
x = pow(x, pe.prime, n) |
|||
o++ |
|||
V o1 = pow(pe.prime, BigInt(o)) |
|||
o1 I/= gcd(mo, o1) |
|||
mo *= o1 |
|||
R mo |
|||
F moTest(a, n) |
|||
I bits:length(a) < 100 |
|||
print(‘ord(’a‘)’, end' ‘’) |
|||
E |
|||
print(‘ord([big])’, end' ‘’) |
|||
print(‘ mod ’n‘ = ’moBachShallit58(a, n, factor(n - 1))) |
|||
moTest(37, 3343) |
|||
moTest(pow(BigInt(10), 100) + 1, 7919) |
|||
moTest(pow(BigInt(10), 1000) + 1, 15485863) |
|||
moTest(pow(BigInt(10), 10000) - 1, BigInt(22801763489)) |
|||
moTest(1511678068, 7379191741) |
|||
moTest(BigInt(‘3047753288’), BigInt(‘2257683301’))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
ord(37) mod 3343 = 1114 |
|||
ord([big]) mod 7919 = 3959 |
|||
ord([big]) mod 15485863 = 15485862 |
|||
ord([big]) mod 22801763489 = 22801763488 |
|||
ord(1511678068) mod 7379191741 = 614932645 |
|||
ord(3047753288) mod 2257683301 = 62713425 |
|||
</pre> |
|||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
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"). |
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"). |
||
< |
<syntaxhighlight lang="ada">package Multiplicative_Order is |
||
type Positive_Array is array (Positive range <>) of Positive; |
type Positive_Array is array (Positive range <>) of Positive; |
||
Line 68: | Line 154: | ||
-- (2) 1 < Coprime_Factors(I) for all I |
-- (2) 1 < Coprime_Factors(I) for all I |
||
end Multiplicative_Order;</ |
end Multiplicative_Order;</syntaxhighlight> |
||
Here is the implementation ("multiplicative_order.adb"): |
Here is the implementation ("multiplicative_order.adb"): |
||
< |
<syntaxhighlight lang="ada">package body Multiplicative_Order is |
||
function Find_Order(Element, Modulus: Positive) return Positive is |
function Find_Order(Element, Modulus: Positive) return Positive is |
||
Line 137: | Line 223: | ||
end Find_Order; |
end Find_Order; |
||
end Multiplicative_Order;</ |
end Multiplicative_Order;</syntaxhighlight> |
||
This is a sample program using the Multiplicative_Order package: |
This is a sample program using the Multiplicative_Order package: |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO, Multiplicative_Order; |
||
procedure Main is |
procedure Main is |
||
Line 163: | Line 249: | ||
-- it gives an incorrect result |
-- it gives an incorrect result |
||
IIO.Put(Find_Order(37, (11, 19, 8, 2))); |
IIO.Put(Find_Order(37, (11, 19, 8, 2))); |
||
end Main;</ |
end Main;</syntaxhighlight> |
||
The output from the sample program: |
The output from the sample program: |
||
Line 178: | Line 264: | ||
{{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 }} --> |
||
< |
<syntaxhighlight lang="algol68">MODE LOOPINT = INT; |
||
MODE POWMODSTRUCT = LONG INT; |
MODE POWMODSTRUCT = LONG INT; |
||
Line 297: | Line 383: | ||
printf(($g$, "Everything checks.", $l$)) |
printf(($g$, "Everything checks.", $l$)) |
||
FI |
FI |
||
)</ |
)</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 310: | Line 396: | ||
=={{header|C}}== |
=={{header|C}}== |
||
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. |
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. |
||
< |
<syntaxhighlight lang="c">ulong mpow(ulong a, ulong p, ulong m) |
||
{ |
{ |
||
ulong r = 1; |
ulong r = 1; |
||
Line 369: | Line 455: | ||
{ |
{ |
||
sieve(); |
sieve(); |
||
printf("%lu\n", multi_order(37, 1000)); |
printf("The multiplicative order of %d related to %d is %lu \n", 37, 1000, multi_order(37, 1000)); |
||
printf("%lu\n", multi_order(54, 100001)); |
printf("The multiplicative order of %d related to %d is %lu \n", 54, 100001, multi_order(54, 100001)); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using System.Collections.Generic; |
using System.Collections.Generic; |
||
using System.Numerics; |
using System.Numerics; |
||
Line 564: | Line 650: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>ord(37) mod 3343 = 1114 |
<pre>ord(37) mod 3343 = 1114 |
||
Line 575: | Line 661: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="cpp">#include <algorithm> |
||
#include <bitset> |
#include <bitset> |
||
#include <iostream> |
#include <iostream> |
||
Line 713: | Line 799: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>100 |
<pre>100 |
||
Line 720: | Line 806: | ||
=={{header|Clojure}}== |
=={{header|Clojure}}== |
||
Translation of Julie, then revised to be more clojure idiomatic. It references some external modules for factoring and integer exponentiation. |
Translation of Julie, then revised to be more clojure idiomatic. It references some external modules for factoring and integer exponentiation. |
||
< |
<syntaxhighlight lang="clojure">(defn gcd [a b] |
||
(if (zero? b) |
(if (zero? b) |
||
a |
a |
||
Line 746: | Line 832: | ||
(map (partial ord' a)) |
(map (partial ord' a)) |
||
(reduce lcm)))) |
(reduce lcm)))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 755: | Line 841: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="d">import std.bigint; |
||
import std.random; |
import std.random; |
||
import std.stdio; |
import std.stdio; |
||
Line 910: | Line 996: | ||
moTest(1511678068, 7379191741); |
moTest(1511678068, 7379191741); |
||
moTest(3047753288, 2257683301); |
moTest(3047753288, 2257683301); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>ord(37) mod 3343 = 1114 |
<pre>ord(37) mod 3343 = 1114 |
||
Line 920: | Line 1,006: | ||
=={{header|EchoLisp}}== |
=={{header|EchoLisp}}== |
||
< |
<syntaxhighlight lang="scheme"> |
||
(require 'bigint) |
(require 'bigint) |
||
Line 953: | Line 1,039: | ||
(order (+ (expt 10 1000) 1) 15485863) |
(order (+ (expt 10 1000) 1) 15485863) |
||
→ 15485862 |
→ 15485862 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{works with|Factor|0.99 2020-01-23}} |
{{works with|Factor|0.99 2020-01-23}} |
||
< |
<syntaxhighlight lang="factor">USING: kernel math math.functions math.primes.factors sequences ; |
||
: (ord) ( a pair -- n ) |
: (ord) ( a pair -- n ) |
||
Line 967: | Line 1,053: | ||
2dup gcd nip 1 = |
2dup gcd nip 1 = |
||
[ group-factors [ (ord) ] with [ lcm ] map-reduce ] |
[ group-factors [ (ord) ] with [ lcm ] map-reduce ] |
||
[ 2drop 0/0. ] if ;</ |
[ 2drop 0/0. ] if ;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 977: | Line 1,063: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 1,078: | Line 1,164: | ||
} |
} |
||
return a.SetInt64(0) |
return a.SetInt64(0) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,093: | Line 1,179: | ||
Assuming a function |
Assuming a function |
||
< |
<syntaxhighlight lang="haskell">powerMod |
||
:: (Integral a, Integral b) |
:: (Integral a, Integral b) |
||
=> a -> a -> b -> a |
=> a -> a -> b -> a |
||
Line 1,107: | Line 1,193: | ||
| even i = g (b * b `rem` m) (i `quot` 2) |
| even i = g (b * b `rem` m) (i `quot` 2) |
||
| otherwise = f b (i - 1) (b * y `rem` m) |
| otherwise = f b (i - 1) (b * y `rem` m) |
||
powerMod m _ _ = error "powerMod: negative exponent"</ |
powerMod m _ _ = error "powerMod: negative exponent"</syntaxhighlight> |
||
to efficiently calculate powers modulo some <code>Integral</code>, we get |
to efficiently calculate powers modulo some <code>Integral</code>, we get |
||
< |
<syntaxhighlight lang="haskell">import Data.List (foldl1') --' |
||
foldl1_ = foldl1' --' |
foldl1_ = foldl1' --' |
||
Line 1,127: | Line 1,213: | ||
where |
where |
||
x = powerMod pk a (t `div` (q ^ e)) |
x = powerMod pk a (t `div` (q ^ e)) |
||
d = length $ takeWhile (/= 1) $ iterate (\y -> powerMod pk y q) x</ |
d = length $ takeWhile (/= 1) $ iterate (\y -> powerMod pk y q) x</syntaxhighlight> |
||
=={{header|J}}== |
=={{header|J}}== |
||
Line 1,133: | Line 1,219: | ||
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. |
||
< |
<syntaxhighlight lang="j">mo=: 4 : 0 |
||
a=. x: x |
a=. x: x |
||
m=. x: y |
m=. x: y |
||
assert. 1=a+.m |
assert. 1=a+.m |
||
*./ a mopk"1 |: __ q: m |
*./ a mopk"1 |: __ q: m |
||
)</ |
)</syntaxhighlight> |
||
The dyadic verb ''mopk'' expects a pair of prime and exponent |
The dyadic verb ''mopk'' expects a pair of prime and exponent |
||
Line 1,147: | Line 1,233: | ||
exponents ''q^d'' into a product. |
exponents ''q^d'' into a product. |
||
< |
<syntaxhighlight lang="j">mopk=: 4 : 0 |
||
a=. x: x |
a=. x: x |
||
'p k'=. x: y |
'p k'=. x: y |
||
Line 1,156: | Line 1,242: | ||
d=. (1<x)+x (pm i. 1:)&> (e-1) */\@$&.> q |
d=. (1<x)+x (pm i. 1:)&> (e-1) */\@$&.> q |
||
*/q^d |
*/q^d |
||
)</ |
)</syntaxhighlight> |
||
For example: |
For example: |
||
< |
<syntaxhighlight lang="j"> 37 mo 1000 |
||
100 |
100 |
||
2 mo _1+10^80x |
2 mo _1+10^80x |
||
190174169488577769580266953193403101748804183400400</ |
190174169488577769580266953193403101748804183400400</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
< |
<syntaxhighlight lang="java">import java.math.BigInteger; |
||
import java.util.ArrayList; |
import java.util.ArrayList; |
||
import java.util.List; |
import java.util.List; |
||
Line 1,272: | Line 1,358: | ||
moTest(BigInteger.valueOf(3047753288L), BigInteger.valueOf(2257683301L)); |
moTest(BigInteger.valueOf(3047753288L), BigInteger.valueOf(2257683301L)); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>ord(37) mod 3343 = 1114 |
<pre>ord(37) mod 3343 = 1114 |
||
Line 1,280: | Line 1,366: | ||
ord(1511678068) mod 7379191741 = 614932645 |
ord(1511678068) mod 7379191741 = 614932645 |
||
ord(3047753288) mod 2257683301 = 62713425</pre> |
ord(3047753288) mod 2257683301 = 62713425</pre> |
||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
The Go implementation of jq supports unbounded-precision integer |
|||
arithmetic and so is suitable for the specified tasks. The program |
|||
given here also runs using the C implementation of jq but falters for |
|||
large integers. |
|||
<syntaxhighlight lang="jq"> |
|||
# Part 1: Library functions |
|||
### Counting and integer arithmetic |
|||
def count(s): reduce s as $x (0; .+1); |
|||
# If $j is 0, then an error condition is raised; |
|||
# otherwise, assuming infinite-precision integer arithmetic, |
|||
# if the input and $j are integers, then the result will be an integer. |
|||
def idivide($j): |
|||
(. - (. % $j)) / $j ; |
|||
def idivide($i; $j): |
|||
$i | idivide($j); |
|||
# Emit [dividend, mod] |
|||
def divmod($j): |
|||
(. % $j) as $mod |
|||
| [(. - $mod) / $j, $mod] ; |
|||
# input should be a non-negative integer for accuracy |
|||
# but may be any non-negative finite number |
|||
def isqrt: |
|||
def irt: |
|||
. as $x |
|||
| 1 | until(. > $x; . * 4) as $q |
|||
| {$q, $x, r: 0} |
|||
| until( .q <= 1; |
|||
.q |= idivide(4) |
|||
| .t = .x - .r - .q |
|||
| .r |= idivide(2) |
|||
| if .t >= 0 |
|||
then .x = .t |
|||
| .r += .q |
|||
else . |
|||
end) |
|||
| .r ; |
|||
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0 |
|||
then irt |
|||
else "isqrt requires a non-negative integer for accuracy" | error |
|||
end ; |
|||
# It is assumed that $n >= 0 |
|||
def power($n): |
|||
. as $in |
|||
| reduce range(0;$n) as $i (1; .* $in); |
|||
# For syntactic convenience |
|||
def power($in; $n): $in | power($n); |
|||
def gcd(a; b): |
|||
# subfunction expects [a,b] as input |
|||
# i.e. a ~ .[0] and b ~ .[1] |
|||
def rgcd: if .[1] == 0 then .[0] |
|||
else [.[1], .[0] % .[1]] | rgcd |
|||
end; |
|||
[a,b] | rgcd; |
|||
### Bit arrays and streams |
|||
def rightshift($n): |
|||
reduce range(0;$n) as $i (.; idivide(2)); |
|||
# Convert the input integer to a stream of 0s and 1s, least significant bit first |
|||
def bitwise: |
|||
recurse( if . >= 2 then idivide(2) else empty end) | . % 2; |
|||
def bitLength: count(bitwise); |
|||
def firstBit: |
|||
if . == 0 then empty |
|||
else first( foreach bitwise as $b (-1; .+1; if $b == 1 then . else empty end)) |
|||
end; |
|||
# Return true if the $i-th least-significant bit is 1, and false otherwise |
|||
def testBit($i): |
|||
(nth($i; bitwise) // 0) == 1; |
|||
# Part 2: "modulo" functions |
|||
# The multiplicative inverse of . modulo $n |
|||
def modInv($n): |
|||
. as $in |
|||
| { r: $n, |
|||
newR: length, # abs |
|||
t: 0, |
|||
newT: 1 } |
|||
| until (.newR != 0.; |
|||
idivide(.r; .newR) as $q |
|||
| .lastT = .t |
|||
| .lastR = .r |
|||
| .t = .newT |
|||
| .r = .newR |
|||
| .newT = .lastT - $q*.newT |
|||
| .newR = .lastR - $q*.newR ) |
|||
| if .r != 1 |
|||
then "\($in) and \($n) are not co-prime." | error |
|||
else if (.t < 0) then .t += $n end |
|||
| if ($in < 0) then - .t else .t end |
|||
end; |
|||
# Return . to the power $exp modulo $mod |
|||
def modPow($exp; $mod): |
|||
def isOdd: . % 2 == 1; |
|||
if $mod == 0 then "Cannot take modPow with modulus 0." | error |
|||
else {r: 1, base: (. % $mod), $exp} |
|||
| if .exp < 0 |
|||
then .exp *= -1 |
|||
| .base |= modInv($mod) |
|||
end |
|||
| until ((.exp == 0) or .emit; |
|||
if .base == 0 then .emit = 0 |
|||
else if (.exp | isOdd) then .r = (.r * .base) % $mod end |
|||
| .exp |= idivide(2) |
|||
| .base |= (.*.) % $mod |
|||
end ) |
|||
| (.emit // .r) |
|||
end ; |
|||
# Part 3: Multiplicative order |
|||
def moBachShallit58($a; $n; $pf): |
|||
{n1: ($n - 1), |
|||
mo: 1 } |
|||
| reduce $pf[] as $pe (.; |
|||
(.n1 | idivide($pe.prime | power($pe.exp))) as $y |
|||
| .o = 0 |
|||
| .x = ($a | modPow($y; ($n|length))) |
|||
| until (.x <= 1; |
|||
.x |= modPow($pe.prime; ($n|length) ) |
|||
| .o += 1 ) |
|||
| .o1 = .o |
|||
| .o1 = power($pe.prime;.o1) |
|||
| .o1 = idivide(.o1; gcd(.mo; .o1) ) |
|||
| .mo = .mo * .o1 ) |
|||
| .mo ; |
|||
def factor($n): |
|||
{ pf: [], |
|||
nn: $n, |
|||
e: ($n | firstBit)} |
|||
| if .e > 0 |
|||
then .e as $e |
|||
| .nn |= rightshift($e) |
|||
| .pf = [{prime: 2, exp: .e}] |
|||
end |
|||
| (.nn | isqrt) as $s |
|||
| .d = 3 |
|||
| until (.nn <= 1; |
|||
if .d > $s then .d = .nn end |
|||
| .e = 0 |
|||
| .done = null |
|||
| until( .done; |
|||
.d as $d |
|||
| (.nn | divmod($d)) as $dm |
|||
| if $dm[1] > 0 |
|||
then .done = true |
|||
else .nn = $dm[0] |
|||
| .e += 1 |
|||
end ) |
|||
| if .e > 0 |
|||
then .pf += [{prime: .d, exp: .e}] |
|||
|.s = (.nn|isqrt) |
|||
end |
|||
| .d += 2 |
|||
) |
|||
| .pf ; |
|||
# $n should be prime |
|||
def moTest($a; $n): |
|||
if ($a|bitLength) < 100 then "ord(\($a)) " else "ord([big]) " end + |
|||
if ($n|bitLength) < 100 then "mod \($n) " else "mod [big] " end + |
|||
"= \(moBachShallit58($a; $n; factor($n - 1)))" ; |
|||
moTest(37; 3343), |
|||
moTest(1 + power(10;100); 7919), |
|||
moTest(1 + power(10;100); 15485863), |
|||
moTest(power(10;10000) - 1; 22801763489), |
|||
moTest(1511678068; 7379191741), |
|||
moTest(3047753288; 2257683301) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
ord(37) mod 3343 = 1114 |
|||
ord([big]) mod 7919 = 3959 |
|||
ord([big]) mod 15485863 = 15485862 |
|||
ord([big]) mod 22801763489 = 22801763488 |
|||
ord(1511678068) mod 7379191741 = 614932645 |
|||
ord(3047753288) mod 2257683301 = 62713425 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
(Uses the <code>factors</code> function from [[Factors of an integer#Julia]].) |
(Uses the <code>factors</code> function from [[Factors of an integer#Julia]].) |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
function factors(n) |
function factors(n) |
||
Line 1,307: | Line 1,596: | ||
end |
end |
||
res |
res |
||
end</ |
end</syntaxhighlight> |
||
Example output (using <code>big</code> to employ arbitrary-precision arithmetic where needed): |
Example output (using <code>big</code> to employ arbitrary-precision arithmetic where needed): |
||
Line 1,326: | Line 1,615: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="scala">// version 1.2.10 |
||
import java.math.BigInteger |
import java.math.BigInteger |
||
Line 1,419: | Line 1,708: | ||
moTest(BigInteger("1511678068"), BigInteger("7379191741")) |
moTest(BigInteger("1511678068"), BigInteger("7379191741")) |
||
moTest(BigInteger("3047753288"), BigInteger("2257683301")) |
moTest(BigInteger("3047753288"), BigInteger("2257683301")) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,432: | Line 1,721: | ||
=={{header|Maple}}== |
=={{header|Maple}}== |
||
< |
<syntaxhighlight lang="maple">numtheory:-order( a, n )</syntaxhighlight> |
||
For example, |
For example, |
||
< |
<syntaxhighlight lang="maple">> numtheory:-order( 37, 1000 ); |
||
100</ |
100</syntaxhighlight> |
||
=={{header|Mathematica}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
In Mathematica this is really easy, as this function is built-in: |
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.<br> |
MultiplicativeOrder[k,n] gives the multiplicative order of k modulo n, defined as the smallest integer m such that k^m == 1 mod n.<br> |
||
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.<BR> |
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.<BR> |
||
Examples: |
Examples: |
||
< |
<syntaxhighlight lang="mathematica">MultiplicativeOrder[37, 1000] |
||
MultiplicativeOrder[10^100 + 1, 7919] (*10^3th prime number Prime[1000]*) |
MultiplicativeOrder[10^100 + 1, 7919] (*10^3th prime number Prime[1000]*) |
||
MultiplicativeOrder[10^1000 + 1, 15485863] (*10^6th prime number*) |
MultiplicativeOrder[10^1000 + 1, 15485863] (*10^6th prime number*) |
||
MultiplicativeOrder[10^10000 - 1, 22801763489] (*10^9th prime number*) |
MultiplicativeOrder[10^10000 - 1, 22801763489] (*10^9th prime number*) |
||
MultiplicativeOrder[13, 1 + 10^80] |
MultiplicativeOrder[13, 1 + 10^80] |
||
MultiplicativeOrder[11, 1 + 10^100]</ |
MultiplicativeOrder[11, 1 + 10^100]</syntaxhighlight> |
||
gives back: |
gives back: |
||
<pre>100 |
<pre>100 |
||
Line 1,457: | Line 1,746: | ||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<syntaxhighlight lang="maxima">zn_order(37, 1000); |
||
/* 100 */ |
/* 100 */ |
||
Line 1,473: | Line 1,762: | ||
zn_order(11, 1 + 10^100); |
zn_order(11, 1 + 10^100); |
||
/* 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800 */</ |
/* 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800 */</syntaxhighlight> |
||
=={{header|Nim}}== |
|||
{{trans|Kotlin}} |
|||
{{libheader|bignum}} |
|||
<syntaxhighlight lang="nim">import strformat |
|||
import bignum |
|||
type PExp = tuple[prime: Int; exp: uint] |
|||
let |
|||
one = newInt(1) |
|||
two = newInt(2) |
|||
ten = newInt(10) |
|||
func sqrt(n: Int): Int = |
|||
var s = n |
|||
while true: |
|||
result = s |
|||
s = (n div result + result) shr 1 |
|||
if s >= result: break |
|||
proc factor(n: Int): seq[PExp] = |
|||
var n = n |
|||
var e = 0u |
|||
while n.bit(e) == 0: inc e |
|||
if e != 0: |
|||
n = n shr e |
|||
result.add (two, e) |
|||
var s = sqrt(n) |
|||
var d = newInt(3) |
|||
while n > one: |
|||
if d > s: d = n |
|||
e = 0u |
|||
while true: |
|||
let (q, r) = divMod(n, d) |
|||
if not r.isZero: break |
|||
n = q |
|||
inc e |
|||
if e != 0: |
|||
result.add (d.clone, e) |
|||
s = sqrt(n) |
|||
inc d, two |
|||
proc moBachShallit58(a, n: Int; pf: seq[PExp]): Int = |
|||
let n = abs(n) |
|||
let n1 = n - one |
|||
result = newInt(1) |
|||
for pe in pf: |
|||
let y = n1 div pe.prime.pow(pe.exp) |
|||
var o = 0u |
|||
var x = a.exp(y.toInt.uint, n) |
|||
while x > one: |
|||
x = x.exp(pe.prime.toInt.uint, n) |
|||
inc o |
|||
var o1 = pe.prime.pow(o) |
|||
o1 = o1 div gcd(result, o1) |
|||
result *= o1 |
|||
proc moTest(a, n: Int) = |
|||
if n.probablyPrime(25) == 0: |
|||
echo "Not computed. Modulus must be prime for this algorithm." |
|||
return |
|||
stdout.write if a.bitLen < 100: &"ord({a})" else: "ord([big])" |
|||
stdout.write if n.bitlen < 100: &" mod {n}" else: " mod [big]" |
|||
let mob = moBachShallit58(a, n, factor(n - one)) |
|||
echo &" = {mob}" |
|||
when isMainModule: |
|||
moTest(newInt(37), newInt(3343)) |
|||
var b = ten.pow(100) + one |
|||
motest(b, newInt(7919)) |
|||
b = ten.pow(1000) + one |
|||
moTest(b, newInt("15485863")) |
|||
b = ten.pow(10000) - one |
|||
moTest(b, newInt("22801763489")) |
|||
moTest(newInt("1511678068"), newInt("7379191741")) |
|||
moTest(newInt("3047753288"), newInt("2257683301"))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>ord(37) mod 3343 = 1114 |
|||
ord([big]) mod 7919 = 3959 |
|||
ord([big]) mod 15485863 = 15485862 |
|||
ord([big]) mod 22801763489 = 22801763488 |
|||
ord(1511678068) mod 7379191741 = 614932645 |
|||
ord(3047753288) mod 2257683301 = 62713425</pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
<lang |
<syntaxhighlight lang="parigp">znorder(Mod(a,n))</syntaxhighlight> |
||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
Using modules: |
Using modules: |
||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
< |
<syntaxhighlight lang="perl">use ntheory qw/znorder/; |
||
say znorder(54, 100001); |
say znorder(54, 100001); |
||
use bigint; say znorder(11, 1 + 10**100);</ |
use bigint; say znorder(11, 1 + 10**100);</syntaxhighlight> |
||
or |
or |
||
< |
<syntaxhighlight lang="perl">use Math::Pari qw/znorder Mod/; |
||
say znorder(Mod(54, 100001)); |
say znorder(Mod(54, 100001)); |
||
say znorder(Mod(11, 1 + Math::Pari::PARI(10)**100));</ |
say znorder(Mod(11, 1 + Math::Pari::PARI(10)**100));</syntaxhighlight> |
||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Ruby}} |
{{trans|Ruby}} |
||
{{libheader|mpfr}} |
{{libheader|Phix/mpfr}} |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>include mpfr.e |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
procedure multi_order(mpz res, a, sequence p_and_k) |
|||
mpz pk = mpz_init(), |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">multi_order</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">res</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">p_and_k</span><span style="color: #0000FF;">)</span> |
|||
t = mpz_init(), |
|||
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span> |
|||
x = mpz_init(), |
|||
<span style="color: #7060A8;">mpz_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
q = mpz_init() |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p_and_k</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
mpz_set_si(res,1) |
|||
<span style="color: #004080;">string</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">ps</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p_and_k</span> |
|||
if length(p_and_k)=1 then |
|||
<span style="color: #7060A8;">mpz_set_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ps</span><span style="color: #0000FF;">)</span> |
|||
string {ps} = p_and_k |
|||
<span style="color: #7060A8;">mpz_sub_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
mpz_set_str(pk,ps) |
|||
<span style="color: #008080;">else</span> |
|||
mpz_sub_ui(t,pk,1) |
|||
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p_and_k</span> |
|||
else |
|||
<span style="color: #7060A8;">mpz_set_d</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
atom {p, k} = p_and_k |
|||
<span style="color: #7060A8;">mpz_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
mpz_ui_pow_ui(pk,p,k) |
|||
<span style="color: #7060A8;">mpz_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
mpz_ui_pow_ui(t,p,k-1) |
|||
<span style="color: #7060A8;">mpz_sub_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
mpz_mul_si(t,t,p-1) |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pz</span><span style="color: #0000FF;">)</span> |
|||
end if |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
sequence pf = mpz_prime_factors(t) |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pf</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_prime_factors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
for i=1 to length(pf) do |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pf</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
if length(pf[i])=1 then |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pf</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
string {fs} = pf[i] |
|||
<span style="color: #004080;">string</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pf</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
mpz_set_str(q,fs) |
|||
<span style="color: #7060A8;">mpz_set_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">)</span> |
|||
mpz_set(x,q) |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
else |
|||
<span style="color: #008080;">else</span> |
|||
{integer qi, integer ei} = pf[i] |
|||
<span style="color: #0000FF;">{</span><span style="color: #004080;">integer</span> <span style="color: #000000;">qi</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">ei</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pf</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
mpz_set_si(q,qi) |
|||
<span style="color: #7060A8;">mpz_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">qi</span><span style="color: #0000FF;">)</span> |
|||
mpz_pow_ui(x,q,ei) |
|||
<span style="color: #7060A8;">mpz_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ei</span><span style="color: #0000FF;">)</span> |
|||
end if |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
mpz_fdiv_q(x, t, x) |
|||
<span style="color: #7060A8;">mpz_fdiv_q</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
mpz_powm(x,a,x,pk) |
|||
<span style="color: #7060A8;">mpz_powm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">)</span> |
|||
integer guard = 0 |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">guard</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
while mpz_cmp_si(x,1)!=0 do |
|||
<span style="color: #008080;">while</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span> |
|||
mpz_mul(res,res,q) |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
mpz_powm(x,x,q,pk) |
|||
<span style="color: #7060A8;">mpz_powm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pk</span><span style="color: #0000FF;">)</span> |
|||
guard += 1 |
|||
<span style="color: #000000;">guard</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
if guard>100 then ?9/0 end if -- (increase if rqd) |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">guard</span><span style="color: #0000FF;">></span><span style="color: #000000;">100</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- (increase if rqd)</span> |
|||
end while |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
x = mpz_free(x) |
|||
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_free</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
end procedure |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
function multiplicative_order(mpz a, m) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">multiplicative_order</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
mpz res = mpz_init(1), |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> |
|||
ri = mpz_init() |
|||
<span style="color: #000000;">ri</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span> |
|||
mpz_gcd(ri,a,m) |
|||
<span style="color: #7060A8;">mpz_gcd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ri</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
if mpz_cmp_si(ri,1)!=0 then return "(a,m) not coprime" end if |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ri</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #008000;">"(a,m) not coprime"</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
sequence pf = mpz_prime_factors(m,10000) -- (increase if rqd) |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pf</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_prime_factors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10000</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (increase if rqd)</span> |
|||
for i=1 to length(pf) do |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pf</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
multi_order(ri,a,pf[i]) |
|||
<span style="color: #000000;">multi_order</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ri</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pf</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span> |
|||
mpz_lcm(res,res,ri) |
|||
<span style="color: #7060A8;">mpz_lcm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ri</span><span style="color: #0000FF;">)</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
return mpz_get_str(res) |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
function shorta(mpz n) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">shorta</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
string res = mpz_get_str(n) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
integer lr = length(res) |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">lr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span> |
|||
if lr>80 then |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">lr</span><span style="color: #0000FF;">></span><span style="color: #000000;">80</span> <span style="color: #008080;">then</span> |
|||
res[6..-6] = "..." |
|||
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">6</span><span style="color: #0000FF;">..-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"..."</span> |
|||
res &= sprintf(" (%d digits)",lr) |
|||
<span style="color: #000000;">res</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" (%d digits)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">)</span> |
|||
end if |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
return res |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
procedure mo_test(mpz a, n) |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
string res = multiplicative_order(a, n) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">multiplicative_order</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
printf(1,"ord(%s) mod %s = %s\n",{shorta(a),shorta(n),res}) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"ord(%s) mod %s = %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">shorta</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span><span style="color: #000000;">shorta</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">res</span><span style="color: #0000FF;">})</span> |
|||
end procedure |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
function i(atom i) return mpz_init(i) end function -- (ugh) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> <span style="color: #000080;font-style:italic;">-- (ugh)</span> |
|||
function p10(integer e,i) -- init to 10^e+i |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">p10</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- init to 10^e+i</span> |
|||
mpz res = mpz_init() |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span> |
|||
mpz_ui_pow_ui(res,10,e) |
|||
<span style="color: #7060A8;">mpz_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e</span><span style="color: #0000FF;">)</span> |
|||
mpz_add_si(res,res,i) |
|||
<span style="color: #7060A8;">mpz_add_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
return res |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
atom t = time() |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
mo_test(i(3), i(10)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(37), i(1000)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">37</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(37), i(10000)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">37</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10000</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(37), i(3343)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">37</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3343</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(37), i(3344)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">37</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3344</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(2), i(1000)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">))</span> |
|||
mo_test(p10(100,+1), i(7919)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100</span><span style="color: #0000FF;">,+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">7919</span><span style="color: #0000FF;">))</span> |
|||
mo_test(p10(1000,+1), i(15485863)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">,+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">15485863</span><span style="color: #0000FF;">))</span> |
|||
mo_test(p10(10000,-1), i(22801763489)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10000</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">22801763489</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(1511678068), i(7379191741)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1511678068</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">7379191741</span><span style="color: #0000FF;">))</span> |
|||
mo_test(i(3047753288), i(2257683301)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3047753288</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2257683301</span><span style="color: #0000FF;">))</span> |
|||
?"===" |
|||
<span style="color: #0000FF;">?</span><span style="color: #008000;">"==="</span> |
|||
mpz b = p10(20,-1) |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">20</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
mo_test(i(2), b) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
mo_test(i(17),b) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">17</span><span style="color: #0000FF;">),</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
mo_test(i(54),i(100001)) |
|||
<span style="color: #000000;">mo_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">54</span><span style="color: #0000FF;">),</span><span style="color: #000000;">i</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100001</span><span style="color: #0000FF;">))</span> |
|||
string s9090 = multiplicative_order(mpz_init(54),mpz_init(100001)) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">s9090</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">multiplicative_order</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">54</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100001</span><span style="color: #0000FF;">))</span> |
|||
if s9090!="9090" then ?9/0 end if |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">s9090</span><span style="color: #0000FF;">!=</span><span style="color: #008000;">"9090"</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
mpz m54 = mpz_init(54), |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">m54</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">54</span><span style="color: #0000FF;">),</span> |
|||
m100001 = mpz_init(100001) |
|||
<span style="color: #000000;">m100001</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100001</span><span style="color: #0000FF;">)</span> |
|||
mpz_powm_ui(b,m54,9090,m100001) |
|||
<span style="color: #7060A8;">mpz_powm_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m54</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9090</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m100001</span><span style="color: #0000FF;">)</span> |
|||
printf(1,"%s\n",mpz_get_str(b)) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">))</span> |
|||
bool error = false |
|||
<span style="color: #004080;">bool</span> <span style="color: #000000;">error</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">false</span> |
|||
for r=1 to 9090-1 do |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9090</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
|||
mpz_powm_ui(b,m54,r,m100001) |
|||
<span style="color: #7060A8;">mpz_powm_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m54</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m100001</span><span style="color: #0000FF;">)</span> |
|||
if mpz_cmp_si(b,1)=0 then |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> |
|||
printf(1,"mpz_powm_ui(54,%d,100001) gives 1!\n",r) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"mpz_powm_ui(54,%d,100001) gives 1!\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
error = true |
|||
<span style="color: #000000;">error</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">true</span> |
|||
exit |
|||
<span style="color: #008080;">exit</span> |
|||
end if |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
if not error then |
|||
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #000000;">error</span> <span style="color: #008080;">then</span> |
|||
printf(1,"Everything checks. (%s)\n",{elapsed(time()-t)}) |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Everything checks. (%s)\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)})</span> |
|||
end if</lang> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,628: | Line 2,015: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
< |
<syntaxhighlight lang="python">def gcd(a, b): |
||
while b != 0: |
while b != 0: |
||
a, b = b, a % b |
a, b = b, a % b |
||
Line 1,693: | Line 2,080: | ||
print 'Exists a power r < 9090 where pow(54,r,b)==1' |
print 'Exists a power r < 9090 where pow(54,r,b)==1' |
||
else: |
else: |
||
print 'Everything checks.'</ |
print 'Everything checks.'</syntaxhighlight> |
||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
Line 1,700: | Line 2,087: | ||
shown below. |
shown below. |
||
< |
<syntaxhighlight lang="racket"> |
||
#lang racket |
#lang racket |
||
(require math) |
(require math) |
||
Line 1,726: | Line 2,113: | ||
(order (- (expt 10 10000) 1) 22801763489) |
(order (- (expt 10 10000) 1) 22801763489) |
||
(order 13 (+ 1 (expt 10 80))) |
(order 13 (+ 1 (expt 10 80))) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Output: |
Output: |
||
< |
<syntaxhighlight lang="racket"> |
||
100 |
100 |
||
3959 |
3959 |
||
Line 1,734: | Line 2,121: | ||
22801763488 |
22801763488 |
||
109609547199756140150989321269669269476675495992554276140800 |
109609547199756140150989321269669269476675495992554276140800 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
(formerly Perl 6) |
(formerly Perl 6) |
||
<syntaxhighlight lang="raku" line>use Prime::Factor; |
|||
<lang perl6>my @primes := 2, |grep *.is-prime, (3,5,7...*); |
|||
sub factor($a is copy) { |
|||
gather { |
|||
for @primes -> $p { |
|||
my $j = 0; |
|||
while $a %% $p { |
|||
$a div= $p; |
|||
$j++; |
|||
} |
|||
take $p => $j if $j > 0; |
|||
last if $a < $p * $p; |
|||
} |
|||
take $a => 1 if $a > 1; |
|||
} |
|||
} |
|||
sub mo-prime($a, $p, $e) { |
sub mo-prime($a, $p, $e) { |
||
my $m = $p ** $e; |
my $m = $p ** $e; |
||
my $t = ($p - 1) * ($p ** ($e - 1)); # = Phi($p**$e) where $p prime |
my $t = ($p - 1) * ($p ** ($e - 1)); # = Phi($p**$e) where $p prime |
||
my @qs = 1; |
my @qs = 1; |
||
for |
for prime-factors($t).Bag -> $f { |
||
@qs = flat @qs.map(-> $q { (0..$f.value).map(-> $j { $q * $f.key ** $j }) }); |
@qs = flat @qs.map(-> $q { (0..$f.value).map(-> $j { $q * $f.key ** $j }) }); |
||
} |
} |
||
@qs.sort.first |
@qs.sort.first: -> $q { expmod( $a, $q, $m ) == 1 }; |
||
} |
} |
||
sub mo($a, $m) { |
sub mo($a, $m) { |
||
$a gcd $m == 1 |
$a gcd $m == 1 or die "$a and $m are not relatively prime"; |
||
[lcm] flat 1, |
[lcm] flat 1, prime-factors($m).Bag.map: { mo-prime($a, .key, .value) }; |
||
} |
} |
||
multi MAIN('test') { |
|||
use Test; |
use Test; |
||
for (10, 21, 25, 150, 1231, 123141, 34131) -> $n { |
for (10, 21, 25, 150, 1231, 123141, 34131) -> $n { |
||
is ([*] |
is ([*] prime-factors($n).Bag.map( { .key ** .value } )), $n, "$n factors correctly"; |
||
} |
} |
||
is mo(37, 1000), 100, 'mo(37,1000) == 100'; |
is mo(37, 1000), 100, 'mo(37,1000) == 100'; |
||
my $b = 10**20-1; |
my $b = 10**20-1; |
||
is mo(2, $b), 3748806900, 'mo(2,10**20-1) == 3748806900'; |
is mo(2, $b), 3748806900, 'mo(2,10**20-1) == 3748806900'; |
||
is mo(17, $b), 1499522760, 'mo(17,10**20-1) == 1499522760'; |
is mo(17, $b), 1499522760, 'mo(17,10**20-1) == 1499522760'; |
||
$b = 100001; |
$b = 100001; |
||
is mo(54, $b), 9090, 'mo(54,100001) == 9090'; |
is mo(54, $b), 9090, 'mo(54,100001) == 9090'; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>ok 1 - 10 factors correctly |
<pre>ok 1 - 10 factors correctly |
||
Line 1,800: | Line 2,171: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
< |
<syntaxhighlight lang="rexx">/*REXX pgm computes multiplicative order of a minimum integer N such that a^n mod m≡1*/ |
||
wa= 0; wm= 0 /* ═a═ ══m══ */ /*maximum widths of the A and M values.*/ |
wa= 0; wm= 0 /* ═a═ ══m══ */ /*maximum widths of the A and M values.*/ |
||
@.=.; @.1= 3 10 |
@.=.; @.1= 3 10 |
||
Line 1,826: | Line 2,197: | ||
say pad 'a=' right(a,wa) pad "m=" right(m,wm) pad 'multiplicative order:' n |
say pad 'a=' right(a,wa) pad "m=" right(m,wm) pad 'multiplicative order:' n |
||
end /*j*/ |
end /*j*/ |
||
end /*w*/ /*stick a fork in it, we're all done. */</ |
end /*w*/ /*stick a fork in it, we're all done. */</syntaxhighlight> |
||
{{out|output|.}} |
{{out|output|.}} |
||
<pre> |
<pre> |
||
Line 1,839: | Line 2,210: | ||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
< |
<syntaxhighlight lang="ruby">require 'prime' |
||
def powerMod(b, p, m) |
def powerMod(b, p, m) |
||
Line 1,879: | Line 2,250: | ||
else |
else |
||
puts 'Everything checks.' |
puts 'Everything checks.' |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,892: | Line 2,263: | ||
=={{header|Seed7}}== |
=={{header|Seed7}}== |
||
< |
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
||
include "bigint.s7i"; |
include "bigint.s7i"; |
||
Line 2,009: | Line 2,380: | ||
moTest(1511678068_, 7379191741_); |
moTest(1511678068_, 7379191741_); |
||
moTest(3047753288_, 2257683301_); |
moTest(3047753288_, 2257683301_); |
||
end func;</ |
end func;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,024: | Line 2,395: | ||
Built-in: |
Built-in: |
||
< |
<syntaxhighlight lang="ruby">say 37.znorder(1000) #=> 100 |
||
say 54.znorder(100001) #=> 9090</ |
say 54.znorder(100001) #=> 9090</syntaxhighlight> |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="ruby">func mo_prime(a, p, e) { |
||
var m = p**e |
var m = p**e |
||
var t = (p-1)*(p**(e-1)) |
var t = (p-1)*(p**(e-1)) |
||
Line 2,053: | Line 2,424: | ||
say mo(2, b) |
say mo(2, b) |
||
say mo(17, b) |
say mo(17, b) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,065: | Line 2,436: | ||
{{trans|Python}} |
{{trans|Python}} |
||
{{tcllib|struct::list}} |
{{tcllib|struct::list}} |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
package require struct::list |
package require struct::list |
||
Line 2,242: | Line 2,613: | ||
if {$r % 1000 == 0} {puts "$r ..."} |
if {$r % 1000 == 0} {puts "$r ..."} |
||
} |
} |
||
puts "OK, $n is the smallest n such that {$a $n $m} satisfies $lambda"</ |
puts "OK, $n is the smallest n such that {$a $n $m} satisfies $lambda"</syntaxhighlight> |
||
=={{header|Visual Basic .NET}}== |
|||
{{trans|C#}} |
|||
<syntaxhighlight lang="vbnet">Imports System.Numerics |
|||
Imports System.Runtime.CompilerServices |
|||
Imports System.Threading |
|||
Module Module1 |
|||
Private s_gen As New ThreadLocal(Of Random)(Function() New Random()) |
|||
Private Function Gen() |
|||
Return s_gen.Value |
|||
End Function |
|||
<Extension()> |
|||
Public Function IsProbablyPrime(value As BigInteger, Optional witnesses As Integer = 10) As Boolean |
|||
If value <= 1 Then |
|||
Return False |
|||
End If |
|||
If witnesses <= 0 Then |
|||
witnesses = 10 |
|||
End If |
|||
Dim d = value - 1 |
|||
Dim s = 0 |
|||
While d Mod 2 = 0 |
|||
d /= 2 |
|||
s += 1 |
|||
End While |
|||
Dim bytes(value.ToByteArray.LongLength - 1) As Byte |
|||
Dim a As BigInteger |
|||
For i = 1 To witnesses |
|||
Do |
|||
Gen.NextBytes(bytes) |
|||
a = New BigInteger(bytes) |
|||
Loop While a < 2 OrElse a >= value - 2 |
|||
Dim x = BigInteger.ModPow(a, d, value) |
|||
If x = 1 OrElse x = value - 1 Then |
|||
Continue For |
|||
End If |
|||
For r = 1 To s - 1 |
|||
x = BigInteger.ModPow(x, 2, value) |
|||
If x = 1 Then |
|||
Return False |
|||
End If |
|||
If x = value - 1 Then |
|||
Exit For |
|||
End If |
|||
Next |
|||
If x <> value - 1 Then |
|||
Return False |
|||
End If |
|||
Next |
|||
Return True |
|||
End Function |
|||
<Extension()> |
|||
Function Sqrt(self As BigInteger) As BigInteger |
|||
Dim b = self |
|||
While True |
|||
Dim a = b |
|||
b = self / a + a >> 1 |
|||
If b >= a Then |
|||
Return a |
|||
End If |
|||
End While |
|||
Throw New Exception("Should not have happened") |
|||
End Function |
|||
<Extension()> |
|||
Function BitLength(self As BigInteger) As Long |
|||
Dim bi = self |
|||
Dim len = 0L |
|||
While bi <> 0 |
|||
len += 1 |
|||
bi >>= 1 |
|||
End While |
|||
Return len |
|||
End Function |
|||
<Extension()> |
|||
Function BitTest(self As BigInteger, pos As Integer) As Boolean |
|||
Dim arr = self.ToByteArray |
|||
Dim i = pos \ 8 |
|||
Dim m = pos Mod 8 |
|||
If i >= arr.Length Then |
|||
Return False |
|||
End If |
|||
Return (arr(i) And (1 << m)) > 0 |
|||
End Function |
|||
Class PExp |
|||
Sub New(p As BigInteger, e As Integer) |
|||
Prime = p |
|||
Exp = e |
|||
End Sub |
|||
Public ReadOnly Property Prime As BigInteger |
|||
Public ReadOnly Property Exp As Integer |
|||
End Class |
|||
Function MoBachShallit58(a As BigInteger, n As BigInteger, pf As List(Of PExp)) As BigInteger |
|||
Dim n1 = n - 1 |
|||
Dim mo As BigInteger = 1 |
|||
For Each pe In pf |
|||
Dim y = n1 / BigInteger.Pow(pe.Prime, pe.Exp) |
|||
Dim o = 0 |
|||
Dim x = BigInteger.ModPow(a, y, BigInteger.Abs(n)) |
|||
While x > 1 |
|||
x = BigInteger.ModPow(x, pe.Prime, BigInteger.Abs(n)) |
|||
o += 1 |
|||
End While |
|||
Dim o1 = BigInteger.Pow(pe.Prime, o) |
|||
o1 /= BigInteger.GreatestCommonDivisor(mo, o1) |
|||
mo *= o1 |
|||
Next |
|||
Return mo |
|||
End Function |
|||
Function Factor(n As BigInteger) As List(Of PExp) |
|||
Dim pf As New List(Of PExp) |
|||
Dim nn = n |
|||
Dim e = 0 |
|||
While Not nn.BitTest(e) |
|||
e += 1 |
|||
End While |
|||
If e > 0 Then |
|||
nn >>= e |
|||
pf.Add(New PExp(2, e)) |
|||
End If |
|||
Dim s = nn.Sqrt |
|||
Dim d As BigInteger = 3 |
|||
While nn > 1 |
|||
If d > s Then |
|||
d = nn |
|||
End If |
|||
e = 0 |
|||
While True |
|||
Dim remainder As New BigInteger |
|||
Dim div = BigInteger.DivRem(nn, d, remainder) |
|||
If remainder.BitLength > 0 Then |
|||
Exit While |
|||
End If |
|||
nn = div |
|||
e += 1 |
|||
End While |
|||
If e > 0 Then |
|||
pf.Add(New PExp(d, e)) |
|||
s = nn.Sqrt |
|||
End If |
|||
d += 2 |
|||
End While |
|||
Return pf |
|||
End Function |
|||
Sub MoTest(a As BigInteger, n As BigInteger) |
|||
If Not n.IsProbablyPrime(20) Then |
|||
Console.WriteLine("Not computed. Modulus must be prime for this algorithm.") |
|||
Return |
|||
End If |
|||
If a.BitLength < 100 Then |
|||
Console.Write("ord({0})", a) |
|||
Else |
|||
Console.Write("ord([big])") |
|||
End If |
|||
If n.BitLength < 100 Then |
|||
Console.Write(" mod {0}", n) |
|||
Else |
|||
Console.Write(" mod [big]") |
|||
End If |
|||
Dim mob = MoBachShallit58(a, n, Factor(n - 1)) |
|||
Console.WriteLine(" = {0}", mob) |
|||
End Sub |
|||
Sub Main() |
|||
MoTest(37, 3343) |
|||
MoTest(BigInteger.Pow(10, 100) + 1, 7919) |
|||
MoTest(BigInteger.Pow(10, 1000) + 1, 15485863) |
|||
MoTest(BigInteger.Pow(10, 10000) - 1, 22801763489) |
|||
MoTest(1511678068, 7379191741) |
|||
MoTest(3047753288, 2257683301) |
|||
End Sub |
|||
End Module</syntaxhighlight> |
|||
{{out}} |
|||
<pre>ord(37) mod 3343 = 1114 |
|||
ord([big]) mod 7919 = 3959 |
|||
ord([big]) mod 15485863 = 15485862 |
|||
ord([big]) mod 22801763489 = 22801763488 |
|||
ord(1511678068) mod 7379191741 = 614932645 |
|||
ord(3047753288) mod 2257683301 = 62713425</pre> |
|||
=={{header|Wren}}== |
|||
{{trans|Kotlin}} |
|||
{{libheader|Wren-big}} |
|||
<syntaxhighlight lang="wren">import "./big" for BigInt |
|||
class PExp { |
|||
construct new(prime, exp) { |
|||
_prime = prime |
|||
_exp = exp |
|||
} |
|||
prime { _prime } |
|||
exp { _exp } |
|||
} |
|||
var moBachShallit58 = Fn.new { |a, n, pf| |
|||
var n1 = n - BigInt.one |
|||
var mo = BigInt.one |
|||
for (pe in pf) { |
|||
var y = n1 / pe.prime.pow(pe.exp) |
|||
var o = 0 |
|||
var x = a.modPow(y, n.abs) |
|||
while (x > BigInt.one) { |
|||
x = x.modPow(pe.prime, n.abs) |
|||
o = o + 1 |
|||
} |
|||
var o1 = BigInt.new(o) |
|||
o1 = pe.prime.pow(o1) |
|||
o1 = o1 / BigInt.gcd(mo, o1) |
|||
mo = mo * o1 |
|||
} |
|||
return mo |
|||
} |
|||
var factor = Fn.new { |n| |
|||
var pf = [] |
|||
var nn = n.copy() |
|||
var e = 0 |
|||
while (!nn.testBit(e)) e = e + 1 |
|||
if (e > 0) { |
|||
nn = nn >> e |
|||
pf.add(PExp.new(BigInt.two, e)) |
|||
} |
|||
var s = nn.isqrt |
|||
var d = BigInt.three |
|||
while (nn > BigInt.one) { |
|||
if (d > s) d = nn |
|||
e = 0 |
|||
while (true) { |
|||
var dm = nn.divMod(d) |
|||
if (dm[1].bitLength > 0) break |
|||
nn = dm[0] |
|||
e = e + 1 |
|||
} |
|||
if (e > 0) { |
|||
pf.add(PExp.new(d, e)) |
|||
s = nn.isqrt |
|||
} |
|||
d = d + BigInt.two |
|||
} |
|||
return pf |
|||
} |
|||
var moTest = Fn.new { |a, n| |
|||
if (!n.isProbablePrime(10)) { |
|||
System.print("Not computed. Modulus must be prime for this algorithm.") |
|||
return |
|||
} |
|||
System.write((a.bitLength < 100) ? "ord(%(a))" : "ord([big])") |
|||
System.write((n.bitLength < 100) ? " mod %(n) " : "mod([big])") |
|||
var mob = moBachShallit58.call(a, n, factor.call(n - BigInt.one)) |
|||
System.print("= %(mob)") |
|||
} |
|||
moTest.call(BigInt.new(37), BigInt.new(3343)) |
|||
var b = BigInt.ten.pow(100) + BigInt.one |
|||
moTest.call(b, BigInt.new(7919)) |
|||
b = BigInt.ten.pow(1000) + BigInt.one |
|||
moTest.call(b, BigInt.new(15485863)) |
|||
b = BigInt.ten.pow(10000) - BigInt.one |
|||
moTest.call(b, BigInt.new(22801763489)) |
|||
moTest.call(BigInt.new(1511678068), BigInt.new(7379191741)) |
|||
moTest.call(BigInt.new(3047753288), BigInt.new(2257683301))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
ord(37) mod 3343 = 1114 |
|||
ord([big]) mod 7919 = 3959 |
|||
ord([big]) mod 15485863 = 15485862 |
|||
ord([big]) mod 22801763489 = 22801763488 |
|||
ord(1511678068) mod 7379191741 = 614932645 |
|||
ord(3047753288) mod 2257683301 = 62713425 |
|||
</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
Line 2,249: | Line 2,918: | ||
It would probably be nice to memoize the prime numbers but that isn't necessary for this task. |
It would probably be nice to memoize the prime numbers but that isn't necessary for this task. |
||
< |
<syntaxhighlight lang="zkl">var BN =Import("zklBigNum"); |
||
var Sieve=Import("sieve"); |
var Sieve=Import("sieve"); |
||
Line 2,283: | Line 2,952: | ||
} |
} |
||
return(res); |
return(res); |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="zkl">multiOrder(37,1000).println(); |
||
b:=BN(10).pow(20)-1; |
b:=BN(10).pow(20)-1; |
||
multiOrder(2,b).println(); |
multiOrder(2,b).println(); |
||
Line 2,292: | Line 2,961: | ||
[BN(1)..multiOrder(54,b)-1].filter1('wrap(r,b54){b54.powm(r,b)==1},BN(54)) : |
[BN(1)..multiOrder(54,b)-1].filter1('wrap(r,b54){b54.powm(r,b)==1},BN(54)) : |
||
if (_) println("Exists a power r < 9090 where (54^r)%b)==1"); |
if (_) println("Exists a power r < 9090 where (54^r)%b)==1"); |
||
else println("Everything checks.");</ |
else println("Everything checks.");</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 11:28, 18 May 2024
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).
- 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)).
11l
T PExp = (BigInt prime, Int exp)
F isqrt(self)
V b = self
L
V a = b
b = (self I/ a + a) I/ 2
I b >= a
R a
F factor(BigInt n)
[PExp] pf
V nn = n
V b = 0
L ((nn % 2) == 0)
nn I/= 2
b++
I b > 0
pf [+]= PExp(BigInt(2), b)
V s = isqrt(nn)
V d = BigInt(3)
L nn > 1
I d > s
d = nn
V e = 0
L
V (div, rem) = divmod(nn, d)
I bits:length(rem) > 0
L.break
nn = div
e++
I e > 0
pf [+]= PExp(d, e)
s = isqrt(nn)
d += 2
R pf
F moBachShallit58(BigInt a, BigInt n; pf)
V n1 = n - 1
V mo = BigInt(1)
L(pe) pf
V y = n1 I/ pow(pe.prime, BigInt(pe.exp))
V o = 0
V x = pow(a, y, n)
L x > 1
x = pow(x, pe.prime, n)
o++
V o1 = pow(pe.prime, BigInt(o))
o1 I/= gcd(mo, o1)
mo *= o1
R mo
F moTest(a, n)
I bits:length(a) < 100
print(‘ord(’a‘)’, end' ‘’)
E
print(‘ord([big])’, end' ‘’)
print(‘ mod ’n‘ = ’moBachShallit58(a, n, factor(n - 1)))
moTest(37, 3343)
moTest(pow(BigInt(10), 100) + 1, 7919)
moTest(pow(BigInt(10), 1000) + 1, 15485863)
moTest(pow(BigInt(10), 10000) - 1, BigInt(22801763489))
moTest(1511678068, 7379191741)
moTest(BigInt(‘3047753288’), BigInt(‘2257683301’))
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Ada
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").
package Multiplicative_Order is
type Positive_Array is array (Positive range <>) of Positive;
function Find_Order(Element, Modulus: Positive) return Positive;
-- naive algorithm
-- returns the smallest I such that (Element**I) mod Modulus = 1
function Find_Order(Element: Positive;
Coprime_Factors: Positive_Array) return Positive;
-- faster algorithm for the same task
-- computes the order of all Coprime_Factors(I)
-- and returns their least common multiple
-- this gives the same result as Find_Order(Element, Modulus)
-- with Modulus being the product of all the Coprime_Factors(I)
--
-- preconditions: (1) 1 = GCD(Coprime_Factors(I), Coprime_Factors(J))
-- for all pairs I, J with I /= J
-- (2) 1 < Coprime_Factors(I) for all I
end Multiplicative_Order;
Here is the implementation ("multiplicative_order.adb"):
package body Multiplicative_Order is
function Find_Order(Element, Modulus: Positive) return Positive is
function Power(Exp, Pow, M: Positive) return Positive is
-- computes Exp**Pow mod M;
-- note that Ada's native integer exponentiation "**" may overflow on
-- computing Exp**Pow before ever computing the "mod M" part
Result: Positive := 1;
E: Positive := Exp;
P: Natural := Pow;
begin
while P > 0 loop
if P mod 2 = 1 then
Result := (Result * E) mod M;
end if;
E := (E * E) mod M;
P := P / 2;
end loop;
return Result;
end Power;
begin -- Find_Order(Element, Modulus)
for I in 1 .. Modulus loop
if Power(Element, I, Modulus) = 1 then
return Positive(I);
end if;
end loop;
raise Program_Error with
Positive'Image(Element) &" is not coprime to" &Positive'Image(Modulus);
end Find_Order;
function Find_Order(Element: Positive;
Coprime_Factors: Positive_Array) return Positive is
function GCD (A, B : Positive) return Integer is
M : Natural := A;
N : Natural := B;
T : Natural;
begin
while N /= 0 loop
T := M;
M := N;
N ;:= T mod N;
end loop;
return M;
end GCD; -- from http://rosettacode.org/wiki/Least_common_multiple#Ada
function LCM (A, B : Natural) return Integer is
begin
if A = 0 or B = 0 then
return 0;
end if;
return abs (A * B) / Gcd (A, B);
end LCM; -- from http://rosettacode.org/wiki/Least_common_multiple#Ada
Result : Positive := 1;
begin -- Find_Order(Element, Coprime_Factors)
for I in Coprime_Factors'Range loop
Result := LCM(Result, Find_Order(Element, Coprime_Factors(I)));
end loop;
return Result;
end Find_Order;
end Multiplicative_Order;
This is a sample program using the Multiplicative_Order package:
with Ada.Text_IO, Multiplicative_Order;
procedure Main is
package IIO is new Ada.Text_IO.Integer_IO(Integer);
use Multiplicative_Order;
begin
IIO.Put(Find_Order(3,10));
IIO.Put(Find_Order(37,1000));
IIO.Put(Find_Order(37,10_000));
IIO.Put(Find_Order(37, 3343));
IIO.Put(Find_Order(37, 3344));
-- IIO.Put(Find_Order( 2,1000));
--would raise Program_Error, because there is no I with 2**I=1 mod 1000
Ada.Text_IO.New_Line;
IIO.Put(Find_Order(3, (2,5))); -- 3 * 5 = 10
IIO.Put(Find_Order(37, (8, 125))); -- 8 * 125 = 1000
IIO.Put(Find_Order(37, (16, 625))); -- 16 * 625 = 10_000
IIO.Put(Find_Order(37, (1 => 3343))); -- 1-element-array: 3343 is a prime
IIO.Put(Find_Order(37, (11, 19, 16))); -- 11 * 19 * 16 = 3344
-- this violates the precondition, because 8 and 2 are not coprime
-- it gives an incorrect result
IIO.Put(Find_Order(37, (11, 19, 8, 2)));
end Main;
The output from the sample program:
4 100 500 1114 20 4 100 500 1114 20 10
ALGOL 68
MODE LOOPINT = INT;
MODE POWMODSTRUCT = LONG INT;
PR READ "prelude/pow_mod.a68" PR;
MODE SORTSTRUCT = LONG INT;
PR READ "prelude/sort.a68" PR;
MODE GCDSTRUCT = LONG INT;
PR READ "prelude/gcd.a68" PR;
PR READ "prelude/iterator.a68" PR;
PROC is prime = (LONG INT p)BOOL:
( p > 1 |#ANDF# ALL((YIELDBOOL yield)VOID: 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 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 factored = (LONG INT in a, PROC (LONG INT,LONG INT)VOID yield)VOID: (
LONG INT a := in a;
# FOR p IN # primes( # DO #
(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
)
# ) OD # );
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 #,
(LONG INT f0,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
)
# OD # );
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
);
PROC 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 #
(LONG INT next)VOID:
out := diadic(out, next)
# OD # );
out
);
PROC mult order = (LONG INT a, LONG INT m)LONG INT: (
PROC mofs = (YIELDLONGINT yield)VOID:(
# FOR p, count IN # factored(m, # DO #
(LONG INT p, LONG INT count)VOID:
yield(mult0rdr1(a,p,count))
)
# OD # );
reduce(lcm, 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
)
Output:
100 3748806900 1499522760 9090 1 Everything checks.
C
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.
ulong mpow(ulong a, ulong p, ulong m)
{
ulong r = 1;
while (p) {
if ((1 & p)) r = r * a % m;
a = a * a % m;
p >>= 1;
}
return r;
}
ulong ipow(ulong a, ulong p) {
ulong r = 1;
while (p) {
if ((1 & p)) r = r * a;
a *= a;
p >>= 1;
}
return r;
}
ulong gcd(ulong m, ulong n)
{
ulong t;
while (m) { t = m; m = n % m; n = t; }
return n;
}
ulong lcm(ulong m, ulong n)
{
ulong g = gcd(m, n);
return m / g * n;
}
ulong multi_order_p(ulong a, ulong p, ulong e)
{
ulong fac[10000];
ulong m = ipow(p, e);
ulong t = m / p * (p - 1);
int i, len = get_factors(t, fac);
for (i = 0; i < len; i++)
if (mpow(a, fac[i], m) == 1)
return fac[i];
return 0;
}
ulong multi_order(ulong a, ulong m)
{
prime_factor pf[100];
int i, len = get_prime_factors(m, pf);
ulong res = 1;
for (i = 0; i < len; i++)
res = lcm(res, multi_order_p(a, pf[i].p, pf[i].e));
return res;
}
int main()
{
sieve();
printf("The multiplicative order of %d related to %d is %lu \n", 37, 1000, multi_order(37, 1000));
printf("The multiplicative order of %d related to %d is %lu \n", 54, 100001, multi_order(54, 100001));
return 0;
}
C#
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Threading;
namespace MultiplicativeOrder {
// Taken from https://stackoverflow.com/a/33918233
public static class PrimeExtensions {
// Random generator (thread safe)
private static ThreadLocal<Random> s_Gen = new ThreadLocal<Random>(
() => {
return new Random();
}
);
// Random generator (thread safe)
private static Random Gen {
get {
return s_Gen.Value;
}
}
public static bool IsProbablyPrime(this BigInteger value, int witnesses = 10) {
if (value <= 1)
return false;
if (witnesses <= 0)
witnesses = 10;
BigInteger d = value - 1;
int s = 0;
while (d % 2 == 0) {
d /= 2;
s += 1;
}
byte[] bytes = new byte[value.ToByteArray().LongLength];
BigInteger a;
for (int i = 0; i < witnesses; i++) {
do {
Gen.NextBytes(bytes);
a = new BigInteger(bytes);
}
while (a < 2 || a >= value - 2);
BigInteger x = BigInteger.ModPow(a, d, value);
if (x == 1 || x == value - 1)
continue;
for (int r = 1; r < s; r++) {
x = BigInteger.ModPow(x, 2, value);
if (x == 1)
return false;
if (x == value - 1)
break;
}
if (x != value - 1)
return false;
}
return true;
}
}
static class Helper {
public static BigInteger Sqrt(this BigInteger self) {
BigInteger b = self;
while (true) {
BigInteger a = b;
b = self / a + a >> 1;
if (b >= a) return a;
}
}
public static long BitLength(this BigInteger self) {
BigInteger bi = self;
long bitlength = 0;
while (bi != 0) {
bitlength++;
bi >>= 1;
}
return bitlength;
}
public static bool BitTest(this BigInteger self, int pos) {
byte[] arr = self.ToByteArray();
int idx = pos / 8;
int mod = pos % 8;
if (idx >= arr.Length) {
return false;
}
return (arr[idx] & (1 << mod)) > 0;
}
}
class PExp {
public PExp(BigInteger prime, int exp) {
Prime = prime;
Exp = exp;
}
public BigInteger Prime { get; }
public int Exp { get; }
}
class Program {
static void MoTest(BigInteger a, BigInteger n) {
if (!n.IsProbablyPrime(20)) {
Console.WriteLine("Not computed. Modulus must be prime for this algorithm.");
return;
}
if (a.BitLength() < 100) {
Console.Write("ord({0})", a);
} else {
Console.Write("ord([big])");
}
if (n.BitLength() < 100) {
Console.Write(" mod {0} ", n);
} else {
Console.Write(" mod [big] ");
}
BigInteger mob = MoBachShallit58(a, n, Factor(n - 1));
Console.WriteLine("= {0}", mob);
}
static BigInteger MoBachShallit58(BigInteger a, BigInteger n, List<PExp> pf) {
BigInteger n1 = n - 1;
BigInteger mo = 1;
foreach (PExp pe in pf) {
BigInteger y = n1 / BigInteger.Pow(pe.Prime, pe.Exp);
int o = 0;
BigInteger x = BigInteger.ModPow(a, y, BigInteger.Abs(n));
while (x > 1) {
x = BigInteger.ModPow(x, pe.Prime, BigInteger.Abs(n));
o++;
}
BigInteger o1 = BigInteger.Pow(pe.Prime, o);
o1 = o1 / BigInteger.GreatestCommonDivisor(mo, o1);
mo = mo * o1;
}
return mo;
}
static List<PExp> Factor(BigInteger n) {
List<PExp> pf = new List<PExp>();
BigInteger nn = n;
int e = 0;
while (!nn.BitTest(e)) e++;
if (e > 0) {
nn = nn >> e;
pf.Add(new PExp(2, e));
}
BigInteger s = nn.Sqrt();
BigInteger d = 3;
while (nn > 1) {
if (d > s) d = nn;
e = 0;
while (true) {
BigInteger div = BigInteger.DivRem(nn, d, out BigInteger rem);
if (rem.BitLength() > 0) break;
nn = div;
e++;
}
if (e > 0) {
pf.Add(new PExp(d, e));
s = nn.Sqrt();
}
d = d + 2;
}
return pf;
}
static void Main(string[] args) {
MoTest(37, 3343);
MoTest(BigInteger.Pow(10, 100) + 1, 7919);
MoTest(BigInteger.Pow(10, 1000) + 1, 15485863);
MoTest(BigInteger.Pow(10, 10000) - 1, 22801763489);
MoTest(1511678068, 7379191741);
MoTest(3047753288, 2257683301);
}
}
}
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
C++
#include <algorithm>
#include <bitset>
#include <iostream>
#include <vector>
typedef unsigned long ulong;
std::vector<ulong> primes;
typedef struct {
ulong p, e;
} prime_factor; /* prime, exponent */
void sieve() {
/* 65536 = 2^16, so we can factor all 32 bit ints */
constexpr int SIZE = 1 << 16;
std::bitset<SIZE> bits;
bits.flip(); // set all bits
bits.reset(0);
bits.reset(1);
for (int i = 0; i < 256; i++) {
if (bits.test(i)) {
for (int j = i * i; j < SIZE; j += i) {
bits.reset(j);
}
}
}
/* collect primes into a list. slightly faster this way if dealing with large numbers */
for (int i = 0; i < SIZE; i++) {
if (bits.test(i)) {
primes.push_back(i);
}
}
}
auto get_prime_factors(ulong n) {
std::vector<prime_factor> lst;
ulong e, p;
for (ulong i = 0; i < primes.size(); i++) {
p = primes[i];
if (p * p > n) break;
for (e = 0; !(n % p); n /= p, e++);
if (e) {
lst.push_back({ p, e });
}
}
if (n != 1) {
lst.push_back({ n, 1 });
}
return lst;
}
auto get_factors(ulong n) {
auto f = get_prime_factors(n);
std::vector<ulong> lst{ 1 };
size_t len2 = 1;
/* L = (1); L = (L, L * p**(1 .. e)) forall((p, e)) */
for (size_t i = 0; i < f.size(); i++, len2 = lst.size()) {
for (ulong j = 0, p = f[i].p; j < f[i].e; j++, p *= f[i].p) {
for (size_t k = 0; k < len2; k++) {
lst.push_back(lst[k] * p);
}
}
}
std::sort(lst.begin(), lst.end());
return lst;
}
ulong mpow(ulong a, ulong p, ulong m) {
ulong r = 1;
while (p) {
if (p & 1) {
r = r * a % m;
}
a = a * a % m;
p >>= 1;
}
return r;
}
ulong ipow(ulong a, ulong p) {
ulong r = 1;
while (p) {
if (p & 1) r *= a;
a *= a;
p >>= 1;
}
return r;
}
ulong gcd(ulong m, ulong n) {
ulong t;
while (m) {
t = m;
m = n % m;
n = t;
}
return n;
}
ulong lcm(ulong m, ulong n) {
ulong g = gcd(m, n);
return m / g * n;
}
ulong multi_order_p(ulong a, ulong p, ulong e) {
ulong m = ipow(p, e);
ulong t = m / p * (p - 1);
auto fac = get_factors(t);
for (size_t i = 0; i < fac.size(); i++) {
if (mpow(a, fac[i], m) == 1) {
return fac[i];
}
}
return 0;
}
ulong multi_order(ulong a, ulong m) {
auto pf = get_prime_factors(m);
ulong res = 1;
for (size_t i = 0; i < pf.size(); i++) {
res = lcm(res, multi_order_p(a, pf[i].p, pf[i].e));
}
return res;
}
int main() {
sieve();
printf("%lu\n", multi_order(37, 1000)); // expect 100
printf("%lu\n", multi_order(54, 100001)); // expect 9090
return 0;
}
- Output:
100 9090
Clojure
Translation of Julie, then revised to be more clojure idiomatic. It references some external modules for factoring and integer exponentiation.
(defn gcd [a b]
(if (zero? b)
a
(recur b (mod a b))))
(defn lcm [a b]
(/ (* a b) (gcd a b)))
(def NaN (Math/log -1))
(defn ord' [a [p e]]
(let [m (imath/expt p e)
t (* (quot m p) (dec p))]
(loop [dv (factor/divisors t)]
(let [d (first dv)]
(if (= (mmath/expm a d m) 1)
d
(recur (next dv)))))))
(defn ord [a n]
(if (not= (gcd a n) 1)
NaN
(->>
(factor/factorize n)
(map (partial ord' a))
(reduce lcm))))
- Output:
user=> (ord 37 1000) 100
D
import std.bigint;
import std.random;
import std.stdio;
struct PExp {
BigInt prime;
int exp;
}
BigInt gcd(BigInt x, BigInt y) {
if (y == 0) {
return x;
}
return gcd(y, x % y);
}
/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
BigInt modPow(BigInt b, BigInt e, BigInt n) {
if (n == 1) return BigInt(0);
BigInt result = 1;
b = b % n;
while (e > 0) {
if (e % 2 == 1) {
result = (result * b) % n;
}
e >>= 1;
b = (b*b) % n;
}
return result;
}
BigInt pow(long b, long e) {
return pow(BigInt(b), BigInt(e));
}
BigInt pow(BigInt b, BigInt e) {
if (e == 0) {
return BigInt(1);
}
BigInt result = 1;
while (e > 1) {
if (e % 2 == 0) {
b *= b;
e /= 2;
} else {
result *= b;
b *= b;
e = (e - 1) / 2;
}
}
return b * result;
}
BigInt sqrt(BigInt self) {
BigInt b = self;
while (true) {
BigInt a = b;
b = self / a + a >> 1;
if (b >= a) return a;
}
}
long bitLength(BigInt self) {
BigInt bi = self;
long length;
while (bi != 0) {
length++;
bi >>= 1;
}
return length;
}
PExp[] factor(BigInt n) {
PExp[] pf;
BigInt nn = n;
int b = 0;
int e = 1;
while ((nn & e) == 0) {
e <<= 1;
b++;
}
if (b > 0) {
nn = nn >> b;
pf ~= PExp(BigInt(2), b);
}
BigInt s = nn.sqrt();
BigInt d = 3;
while (nn > 1) {
if (d > s) d = nn;
e = 0;
while (true) {
BigInt div, rem;
nn.divMod(d, div, rem);
if (rem.bitLength > 0) break;
nn = div;
e++;
}
if (e > 0) {
pf ~= PExp(d, e);
s = nn.sqrt();
}
d += 2;
}
return pf;
}
BigInt moBachShallit58(BigInt a, BigInt n, PExp[] pf) {
BigInt n1 = n - 1;
BigInt mo = 1;
foreach(pe; pf) {
BigInt y = n1 / pe.prime.pow(BigInt(pe.exp));
int o = 0;
BigInt x = a.modPow(y, n);
while (x > 1) {
x = x.modPow(pe.prime, n);
o++;
}
BigInt o1 = pe.prime.pow(BigInt(o));
o1 = o1 / gcd(mo, o1);
mo = mo * o1;
}
return mo;
}
void moTest(ulong a, ulong n) {
moTest(BigInt(a), n);
}
void moTest(BigInt a, ulong n) {
// Commented out because the implementations tried all failed for the -2 and -3 tests.
// if (!n.isProbablePrime()) {
// writeln("Not computed. Modulus must be prime for this algorithm.");
// return;
// }
if (a.bitLength < 100) {
write("ord(", a, ")");
} else {
write("ord([big])");
}
write(" mod ", n, " ");
BigInt nn = n;
BigInt mob = moBachShallit58(a, nn, factor(nn - 1));
writeln("= ", mob);
}
void main() {
moTest(37, 3343);
moTest(pow(10, 100) + 1, 7919);
moTest(pow(10, 1000) + 1, 15485863);
moTest(pow(10, 10000) - 1, 22801763489);
moTest(1511678068, 7379191741);
moTest(3047753288, 2257683301);
}
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
EchoLisp
(require 'bigint)
;; factor-exp returns a list ((p k) ..) : a = p1^k1 * p2^k2 ..
(define (factor-exp a)
(map (lambda (g) (list (first g) (length g)))
(group* (prime-factors a))))
;; copied from Ruby
(define (_mult_order a p k (x))
(define pk (expt p k))
(define t (* (1- p) (expt p (1- k))))
(define r 1)
(for [((q e) (factor-exp t))]
(set! x (powmod a (/ t (expt q e)) pk))
(while (!= x 1)
(*= r q)
(set! x (powmod x q pk))))
r)
(define (order a m)
"multiplicative order : (order a m) → n : a^n = 1 (mod m)"
(assert (= 1 (gcd a m)) "a and m must be coprimes")
(define mopks (for/list [((p k) (factor-exp m))] (_mult_order a p k)))
(for/fold (n 1) ((mopk mopks)) (lcm n mopk)))
;; results
order 37 1000)
→ 100
(order (+ (expt 10 100) 1) 7919)
→ 3959
(order (+ (expt 10 1000) 1) 15485863)
→ 15485862
Factor
USING: kernel math math.functions math.primes.factors sequences ;
: (ord) ( a pair -- n )
first2 dupd ^ swap dupd [ /i ] keep 1 - * divisors
[ swap ^mod 1 = ] 2with find nip ;
: ord ( a n -- m )
2dup gcd nip 1 =
[ group-factors [ (ord) ] with [ lcm ] map-reduce ]
[ 2drop 0/0. ] if ;
- Output:
IN: scratchpad 37 1000 ord . 100 IN: scratchpad 10 100 ^ 1 + 7919 ord . 3959
Go
package main
import (
"fmt"
"math/big"
)
func main() {
moTest(big.NewInt(37), big.NewInt(3343))
b := big.NewInt(100)
moTest(b.Add(b.Exp(ten, b, nil), one), big.NewInt(7919))
moTest(b.Add(b.Exp(ten, b.SetInt64(1000), nil), one), big.NewInt(15485863))
moTest(b.Sub(b.Exp(ten, b.SetInt64(10000), nil), one),
big.NewInt(22801763489))
moTest(big.NewInt(1511678068), big.NewInt(7379191741))
moTest(big.NewInt(3047753288), big.NewInt(2257683301))
}
func moTest(a, n *big.Int) {
if a.BitLen() < 100 {
fmt.Printf("ord(%v)", a)
} else {
fmt.Print("ord([big])")
}
if n.BitLen() < 100 {
fmt.Printf(" mod %v ", n)
} else {
fmt.Print(" mod [big] ")
}
if !n.ProbablyPrime(20) {
fmt.Println("not computed. modulus must be prime for this algorithm.")
return
}
fmt.Println("=", moBachShallit58(a, n, factor(new(big.Int).Sub(n, one))))
}
var one = big.NewInt(1)
var two = big.NewInt(2)
var ten = big.NewInt(10)
func moBachShallit58(a, n *big.Int, pf []pExp) *big.Int {
n1 := new(big.Int).Sub(n, one)
var x, y, o1, g big.Int
mo := big.NewInt(1)
for _, pe := range pf {
y.Quo(n1, y.Exp(pe.prime, big.NewInt(pe.exp), nil))
var o int64
for x.Exp(a, &y, n); x.Cmp(one) > 0; o++ {
x.Exp(&x, pe.prime, n)
}
o1.Exp(pe.prime, o1.SetInt64(o), nil)
mo.Mul(mo, o1.Quo(&o1, g.GCD(nil, nil, mo, &o1)))
}
return mo
}
type pExp struct {
prime *big.Int
exp int64
}
func factor(n *big.Int) (pf []pExp) {
var e int64
for ; n.Bit(int(e)) == 0; e++ {
}
if e > 0 {
n.Rsh(n, uint(e))
pf = []pExp{{big.NewInt(2), e}}
}
s := sqrt(n)
q, r := new(big.Int), new(big.Int)
for d := big.NewInt(3); n.Cmp(one) > 0; d.Add(d, two) {
if d.Cmp(s) > 0 {
d.Set(n)
}
for e = 0; ; e++ {
q.QuoRem(n, d, r)
if r.BitLen() > 0 {
break
}
n.Set(q)
}
if e > 0 {
pf = append(pf, pExp{new(big.Int).Set(d), e})
s = sqrt(n)
}
}
return
}
func sqrt(n *big.Int) *big.Int {
a := new(big.Int)
for b := new(big.Int).Set(n); ; {
a.Set(b)
b.Rsh(b.Add(b.Quo(n, a), a), 1)
if b.Cmp(a) >= 0 {
return a
}
}
return a.SetInt64(0)
}
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Haskell
Assuming a function
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"
to efficiently calculate powers modulo some Integral
, we get
import Data.List (foldl1') --'
foldl1_ = foldl1' --'
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
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.
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 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.
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:
37 mo 1000
100
2 mo _1+10^80x
190174169488577769580266953193403101748804183400400
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
public class MultiplicativeOrder {
private static final BigInteger ONE = BigInteger.ONE;
private static final BigInteger TWO = BigInteger.valueOf(2);
private static final BigInteger THREE = BigInteger.valueOf(3);
private static final BigInteger TEN = BigInteger.TEN;
private static class PExp {
BigInteger prime;
long exp;
PExp(BigInteger prime, long exp) {
this.prime = prime;
this.exp = exp;
}
}
private static void moTest(BigInteger a, BigInteger n) {
if (!n.isProbablePrime(20)) {
System.out.println("Not computed. Modulus must be prime for this algorithm.");
return;
}
if (a.bitLength() < 100) System.out.printf("ord(%s)", a);
else System.out.print("ord([big])");
if (n.bitLength() < 100) System.out.printf(" mod %s ", n);
else System.out.print(" mod [big] ");
BigInteger mob = moBachShallit58(a, n, factor(n.subtract(ONE)));
System.out.println("= " + mob);
}
private static BigInteger moBachShallit58(BigInteger a, BigInteger n, List<PExp> pf) {
BigInteger n1 = n.subtract(ONE);
BigInteger mo = ONE;
for (PExp pe : pf) {
BigInteger y = n1.divide(pe.prime.pow((int) pe.exp));
long o = 0;
BigInteger x = a.modPow(y, n.abs());
while (x.compareTo(ONE) > 0) {
x = x.modPow(pe.prime, n.abs());
o++;
}
BigInteger o1 = BigInteger.valueOf(o);
o1 = pe.prime.pow(o1.intValue());
o1 = o1.divide(mo.gcd(o1));
mo = mo.multiply(o1);
}
return mo;
}
private static List<PExp> factor(BigInteger n) {
List<PExp> pf = new ArrayList<>();
BigInteger nn = n;
Long e = 0L;
while (!nn.testBit(e.intValue())) e++;
if (e > 0L) {
nn = nn.shiftRight(e.intValue());
pf.add(new PExp(TWO, e));
}
BigInteger s = sqrt(nn);
BigInteger d = THREE;
while (nn.compareTo(ONE) > 0) {
if (d.compareTo(s) > 0) d = nn;
e = 0L;
while (true) {
BigInteger[] qr = nn.divideAndRemainder(d);
if (qr[1].bitLength() > 0) break;
nn = qr[0];
e++;
}
if (e > 0L) {
pf.add(new PExp(d, e));
s = sqrt(nn);
}
d = d.add(TWO);
}
return pf;
}
private static BigInteger sqrt(BigInteger n) {
BigInteger b = n;
while (true) {
BigInteger a = b;
b = n.divide(a).add(a).shiftRight(1);
if (b.compareTo(a) >= 0) return a;
}
}
public static void main(String[] args) {
moTest(BigInteger.valueOf(37), BigInteger.valueOf(3343));
BigInteger b = TEN.pow(100).add(ONE);
moTest(b, BigInteger.valueOf(7919));
b = TEN.pow(1000).add(ONE);
moTest(b, BigInteger.valueOf(15485863));
b = TEN.pow(10000).subtract(ONE);
moTest(b, BigInteger.valueOf(22801763489L));
moTest(BigInteger.valueOf(1511678068), BigInteger.valueOf(7379191741L));
moTest(BigInteger.valueOf(3047753288L), BigInteger.valueOf(2257683301L));
}
}
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
jq
Adapted from Wren
Works with gojq, the Go implementation of jq
The Go implementation of jq supports unbounded-precision integer arithmetic and so is suitable for the specified tasks. The program given here also runs using the C implementation of jq but falters for large integers.
# Part 1: Library functions
### Counting and integer arithmetic
def count(s): reduce s as $x (0; .+1);
# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($j):
(. - (. % $j)) / $j ;
def idivide($i; $j):
$i | idivide($j);
# Emit [dividend, mod]
def divmod($j):
(. % $j) as $mod
| [(. - $mod) / $j, $mod] ;
# input should be a non-negative integer for accuracy
# but may be any non-negative finite number
def isqrt:
def irt:
. as $x
| 1 | until(. > $x; . * 4) as $q
| {$q, $x, r: 0}
| until( .q <= 1;
.q |= idivide(4)
| .t = .x - .r - .q
| .r |= idivide(2)
| if .t >= 0
then .x = .t
| .r += .q
else .
end)
| .r ;
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0
then irt
else "isqrt requires a non-negative integer for accuracy" | error
end ;
# It is assumed that $n >= 0
def power($n):
. as $in
| reduce range(0;$n) as $i (1; .* $in);
# For syntactic convenience
def power($in; $n): $in | power($n);
def gcd(a; b):
# subfunction expects [a,b] as input
# i.e. a ~ .[0] and b ~ .[1]
def rgcd: if .[1] == 0 then .[0]
else [.[1], .[0] % .[1]] | rgcd
end;
[a,b] | rgcd;
### Bit arrays and streams
def rightshift($n):
reduce range(0;$n) as $i (.; idivide(2));
# Convert the input integer to a stream of 0s and 1s, least significant bit first
def bitwise:
recurse( if . >= 2 then idivide(2) else empty end) | . % 2;
def bitLength: count(bitwise);
def firstBit:
if . == 0 then empty
else first( foreach bitwise as $b (-1; .+1; if $b == 1 then . else empty end))
end;
# Return true if the $i-th least-significant bit is 1, and false otherwise
def testBit($i):
(nth($i; bitwise) // 0) == 1;
# Part 2: "modulo" functions
# The multiplicative inverse of . modulo $n
def modInv($n):
. as $in
| { r: $n,
newR: length, # abs
t: 0,
newT: 1 }
| until (.newR != 0.;
idivide(.r; .newR) as $q
| .lastT = .t
| .lastR = .r
| .t = .newT
| .r = .newR
| .newT = .lastT - $q*.newT
| .newR = .lastR - $q*.newR )
| if .r != 1
then "\($in) and \($n) are not co-prime." | error
else if (.t < 0) then .t += $n end
| if ($in < 0) then - .t else .t end
end;
# Return . to the power $exp modulo $mod
def modPow($exp; $mod):
def isOdd: . % 2 == 1;
if $mod == 0 then "Cannot take modPow with modulus 0." | error
else {r: 1, base: (. % $mod), $exp}
| if .exp < 0
then .exp *= -1
| .base |= modInv($mod)
end
| until ((.exp == 0) or .emit;
if .base == 0 then .emit = 0
else if (.exp | isOdd) then .r = (.r * .base) % $mod end
| .exp |= idivide(2)
| .base |= (.*.) % $mod
end )
| (.emit // .r)
end ;
# Part 3: Multiplicative order
def moBachShallit58($a; $n; $pf):
{n1: ($n - 1),
mo: 1 }
| reduce $pf[] as $pe (.;
(.n1 | idivide($pe.prime | power($pe.exp))) as $y
| .o = 0
| .x = ($a | modPow($y; ($n|length)))
| until (.x <= 1;
.x |= modPow($pe.prime; ($n|length) )
| .o += 1 )
| .o1 = .o
| .o1 = power($pe.prime;.o1)
| .o1 = idivide(.o1; gcd(.mo; .o1) )
| .mo = .mo * .o1 )
| .mo ;
def factor($n):
{ pf: [],
nn: $n,
e: ($n | firstBit)}
| if .e > 0
then .e as $e
| .nn |= rightshift($e)
| .pf = [{prime: 2, exp: .e}]
end
| (.nn | isqrt) as $s
| .d = 3
| until (.nn <= 1;
if .d > $s then .d = .nn end
| .e = 0
| .done = null
| until( .done;
.d as $d
| (.nn | divmod($d)) as $dm
| if $dm[1] > 0
then .done = true
else .nn = $dm[0]
| .e += 1
end )
| if .e > 0
then .pf += [{prime: .d, exp: .e}]
|.s = (.nn|isqrt)
end
| .d += 2
)
| .pf ;
# $n should be prime
def moTest($a; $n):
if ($a|bitLength) < 100 then "ord(\($a)) " else "ord([big]) " end +
if ($n|bitLength) < 100 then "mod \($n) " else "mod [big] " end +
"= \(moBachShallit58($a; $n; factor($n - 1)))" ;
moTest(37; 3343),
moTest(1 + power(10;100); 7919),
moTest(1 + power(10;100); 15485863),
moTest(power(10;10000) - 1; 22801763489),
moTest(1511678068; 7379191741),
moTest(3047753288; 2257683301)
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Julia
(Uses the factors
function from Factors of an integer#Julia.)
using Primes
function factors(n)
f = [one(n)]
for (p,e) in factor(n)
f = reduce(vcat, [f*p^j for j in 1:e], init=f)
end
return length(f) == 1 ? [one(n), n] : sort!(f)
end
function multorder(a, m)
gcd(a,m) == 1 || error("$a and $m are not coprime")
res = one(m)
for (p,e) in factor(m)
m = p^e
t = div(m, p) * (p-1)
for f in factors(t)
if powermod(a, f, m) == 1
res = lcm(res, f)
break
end
end
end
res
end
Example output (using big
to employ arbitrary-precision arithmetic where needed):
julia> multorder(37, 1000) 100 julia> multorder(big(10)^100 + 1, 7919) 3959 julia> multorder(big(10)^1000 + 1, 15485863) 15485862 julia> multorder(big(10)^10000 - 1, 22801763489) 22801763488
Kotlin
// version 1.2.10
import java.math.BigInteger
val bigOne = BigInteger.ONE
val bigTwo = 2.toBigInteger()
val bigThree = 3.toBigInteger()
val bigTen = BigInteger.TEN
class PExp(val prime: BigInteger, val exp: Long)
fun moTest(a: BigInteger, n: BigInteger) {
if (!n.isProbablePrime(20)) {
println("Not computed. Modulus must be prime for this algorithm.")
return
}
if (a.bitLength() < 100) print("ord($a)") else print("ord([big])")
if (n.bitLength() < 100) print(" mod $n ") else print(" mod [big] ")
val mob = moBachShallit58(a, n, factor(n - bigOne))
println("= $mob")
}
fun moBachShallit58(a: BigInteger, n: BigInteger, pf: List<PExp>): BigInteger {
val n1 = n - bigOne
var mo = bigOne
for (pe in pf) {
val y = n1 / pe.prime.pow(pe.exp.toInt())
var o = 0L
var x = a.modPow(y, n.abs())
while (x > bigOne) {
x = x.modPow(pe.prime, n.abs())
o++
}
var o1 = o.toBigInteger()
o1 = pe.prime.pow(o1.toInt())
o1 /= mo.gcd(o1)
mo *= o1
}
return mo
}
fun factor(n: BigInteger): List<PExp> {
val pf = mutableListOf<PExp>()
var nn = n
var e = 0L
while (!nn.testBit(e.toInt())) e++
if (e > 0L) {
nn = nn shr e.toInt()
pf.add(PExp(bigTwo, e))
}
var s = bigSqrt(nn)
var d = bigThree
while (nn > bigOne) {
if (d > s) d = nn
e = 0L
while (true) {
val (q, r) = nn.divideAndRemainder(d)
if (r.bitLength() > 0) break
nn = q
e++
}
if (e > 0L) {
pf.add(PExp(d, e))
s = bigSqrt(nn)
}
d += bigTwo
}
return pf
}
fun bigSqrt(n: BigInteger): BigInteger {
var b = n
while (true) {
val a = b
b = (n / a + a) shr 1
if (b >= a) return a
}
}
fun main(args: Array<String>) {
moTest(37.toBigInteger(), 3343.toBigInteger())
var b = bigTen.pow(100) + bigOne
moTest(b, 7919.toBigInteger())
b = bigTen.pow(1000) + bigOne
moTest(b, BigInteger("15485863"))
b = bigTen.pow(10000) - bigOne
moTest(b, BigInteger("22801763489"))
moTest(BigInteger("1511678068"), BigInteger("7379191741"))
moTest(BigInteger("3047753288"), BigInteger("2257683301"))
}
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Maple
numtheory:-order( a, n )
For example,
> numtheory:-order( 37, 1000 );
100
Mathematica/Wolfram Language
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:
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]
gives back:
100 3959 15485862 22801763488 109609547199756140150989321269669269476675495992554276140800 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800
Maxima
zn_order(37, 1000);
/* 100 */
zn_order(10^100 + 1, 7919);
/* 3959 */
zn_order(10^1000 + 1, 15485863);
/* 15485862 */
zn_order(10^10000 - 1, 22801763489);
/* 22801763488 */
zn_order(13, 1 + 10^80);
/* 109609547199756140150989321269669269476675495992554276140800 */
zn_order(11, 1 + 10^100);
/* 2583496112724752500580158969425549088007844580826869433740066152289289764829816356800 */
Nim
import strformat
import bignum
type PExp = tuple[prime: Int; exp: uint]
let
one = newInt(1)
two = newInt(2)
ten = newInt(10)
func sqrt(n: Int): Int =
var s = n
while true:
result = s
s = (n div result + result) shr 1
if s >= result: break
proc factor(n: Int): seq[PExp] =
var n = n
var e = 0u
while n.bit(e) == 0: inc e
if e != 0:
n = n shr e
result.add (two, e)
var s = sqrt(n)
var d = newInt(3)
while n > one:
if d > s: d = n
e = 0u
while true:
let (q, r) = divMod(n, d)
if not r.isZero: break
n = q
inc e
if e != 0:
result.add (d.clone, e)
s = sqrt(n)
inc d, two
proc moBachShallit58(a, n: Int; pf: seq[PExp]): Int =
let n = abs(n)
let n1 = n - one
result = newInt(1)
for pe in pf:
let y = n1 div pe.prime.pow(pe.exp)
var o = 0u
var x = a.exp(y.toInt.uint, n)
while x > one:
x = x.exp(pe.prime.toInt.uint, n)
inc o
var o1 = pe.prime.pow(o)
o1 = o1 div gcd(result, o1)
result *= o1
proc moTest(a, n: Int) =
if n.probablyPrime(25) == 0:
echo "Not computed. Modulus must be prime for this algorithm."
return
stdout.write if a.bitLen < 100: &"ord({a})" else: "ord([big])"
stdout.write if n.bitlen < 100: &" mod {n}" else: " mod [big]"
let mob = moBachShallit58(a, n, factor(n - one))
echo &" = {mob}"
when isMainModule:
moTest(newInt(37), newInt(3343))
var b = ten.pow(100) + one
motest(b, newInt(7919))
b = ten.pow(1000) + one
moTest(b, newInt("15485863"))
b = ten.pow(10000) - one
moTest(b, newInt("22801763489"))
moTest(newInt("1511678068"), newInt("7379191741"))
moTest(newInt("3047753288"), newInt("2257683301"))
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
PARI/GP
znorder(Mod(a,n))
Perl
Using modules:
use ntheory qw/znorder/;
say znorder(54, 100001);
use bigint; say znorder(11, 1 + 10**100);
or
use Math::Pari qw/znorder Mod/;
say znorder(Mod(54, 100001));
say znorder(Mod(11, 1 + Math::Pari::PARI(10)**100));
Phix
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
- Output:
ord(3) mod 10 = 4 ord(37) mod 1000 = 100 ord(37) mod 10000 = 500 ord(37) mod 3343 = 1114 ord(37) mod 3344 = 20 ord(2) mod 1000 = (a,m) not coprime ord(10000...00001 (101 digits)) mod 7919 = 3959 ord(10000...00001 (1001 digits)) mod 15485863 = 15485862 ord(99999...99999 (10000 digits)) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425 "===" ord(2) mod 99999999999999999999 = 3748806900 ord(17) mod 99999999999999999999 = 1499522760 ord(54) mod 100001 = 9090 1 Everything checks. (0.2s)
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.'
Racket
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.
#lang racket
(require math)
(define (order a n)
(unless (coprime? a n) (error 'order "arguments must be coprime"))
(for/fold ([o 1]) ([r (factorize n)])
(lcm o (order1 a r))))
(define (order1 a p&e)
(match-define (list p e) p&e)
(define m (expt p e))
(define t (* (- p 1) (expt p (- e 1))))
(define qs
(for/fold ([qs '(1)]) ([f (factorize t)])
(match f [(list f0 f1)
(for*/list ([q qs] [j (in-range (+ 1 f1))])
(* q (expt f0 j)))])))
(for/or ([q (sort qs <)] #:when (= (modular-expt a q m) 1)) q))
(order 37 1000)
(order (+ (expt 10 100) 1) 7919)
(order (+ (expt 10 1000) 1) 15485863)
(order (- (expt 10 10000) 1) 22801763489)
(order 13 (+ 1 (expt 10 80)))
Output:
100
3959
15485862
22801763488
109609547199756140150989321269669269476675495992554276140800
Raku
(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
REXX
/*REXX pgm computes multiplicative order of a minimum integer N such that a^n mod m≡1*/
wa= 0; wm= 0 /* ═a═ ══m══ */ /*maximum widths of the A and M values.*/
@.=.; @.1= 3 10
@.2= 37 1000
@.3= 37 10000
@.4= 37 3343
@.5= 37 3344
@.6= 2 1000
pad= left('', 9)
d= 500 /*use 500 decimal digits for a starter.*/
do w=1 for 2 /*when W≡1, find max widths of A and M.*/
do j=1 while @.j\==.; parse var @.j a . 1 r m , n
if w==1 then do; wa= max(wa, length(a) ); wm= max(wm, length(m) ); iterate
end
if m//a==0 then n= ' [solution not possible]' /*test co─prime for A and B. */
numeric digits d /*start with 100 decimal digits. */
if n=='' then do n= 2; p= r * a /*compute product──may have an exponent*/
parse var p 'E' _ /*try to extract the exponent from P. */
if _\=='' then do; numeric digits _+d /*bump the decimal digs.*/
p=r*a /*recalculate integer P.*/
end
if p//m==1 then leave /*now, perform the nitty─gritty modulo.*/
r= p /*assign product to R for next multiply*/
end /*n*/ /* [↑] // is really ÷ remainder.*/
say pad 'a=' right(a,wa) pad "m=" right(m,wm) pad 'multiplicative order:' n
end /*j*/
end /*w*/ /*stick a fork in it, we're all done. */
- output:
a= 3 m= 10 multiplicative order: 4 a= 37 m= 1000 multiplicative order: 100 a= 37 m= 10000 multiplicative order: 500 a= 37 m= 3343 multiplicative order: 1114 a= 37 m= 3344 multiplicative order: 20 a= 2 m= 1000 multiplicative order: [solution not possible]
Ruby
require 'prime'
def powerMod(b, p, m)
p.to_s(2).each_char.inject(1) do |result, bit|
result = (result * result) % m
bit=='1' ? (result * b) % m : result
end
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) do |result, f|
result.lcm(multOrder_(a, *f))
end
end
puts multOrder(37, 1000)
b = 10**20-1
puts multOrder(2, b)
puts multOrder(17,b)
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
- Output:
100 3748806900 1499522760 9090 1 Everything checks.
Seed7
$ 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;
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Sidef
Built-in:
say 37.znorder(1000) #=> 100
say 54.znorder(100001) #=> 9090
func mo_prime(a, p, e) {
var m = p**e
var t = (p-1)*(p**(e-1))
var qs = [1]
for f in (t.factor_exp) {
qs.map! {|q|
0..f[1] -> map {|j| q * f[0]**j }...
}
}
qs.sort.first_by {|q| powmod(a, q, m) == 1 }
}
func mo(a, m) {
gcd(a, m) == 1 || die "#{a} and #{m} are not relatively prime"
Math.lcm(1, m.factor_exp.map {|r| mo_prime(a, r...) }...)
}
say mo(37, 1000)
say mo(54, 100001)
with (10**20 - 1) {|b|
say mo(2, b)
say mo(17, b)
}
- Output:
100 9090 3748806900 1499522760
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"
Visual Basic .NET
Imports System.Numerics
Imports System.Runtime.CompilerServices
Imports System.Threading
Module Module1
Private s_gen As New ThreadLocal(Of Random)(Function() New Random())
Private Function Gen()
Return s_gen.Value
End Function
<Extension()>
Public Function IsProbablyPrime(value As BigInteger, Optional witnesses As Integer = 10) As Boolean
If value <= 1 Then
Return False
End If
If witnesses <= 0 Then
witnesses = 10
End If
Dim d = value - 1
Dim s = 0
While d Mod 2 = 0
d /= 2
s += 1
End While
Dim bytes(value.ToByteArray.LongLength - 1) As Byte
Dim a As BigInteger
For i = 1 To witnesses
Do
Gen.NextBytes(bytes)
a = New BigInteger(bytes)
Loop While a < 2 OrElse a >= value - 2
Dim x = BigInteger.ModPow(a, d, value)
If x = 1 OrElse x = value - 1 Then
Continue For
End If
For r = 1 To s - 1
x = BigInteger.ModPow(x, 2, value)
If x = 1 Then
Return False
End If
If x = value - 1 Then
Exit For
End If
Next
If x <> value - 1 Then
Return False
End If
Next
Return True
End Function
<Extension()>
Function Sqrt(self As BigInteger) As BigInteger
Dim b = self
While True
Dim a = b
b = self / a + a >> 1
If b >= a Then
Return a
End If
End While
Throw New Exception("Should not have happened")
End Function
<Extension()>
Function BitLength(self As BigInteger) As Long
Dim bi = self
Dim len = 0L
While bi <> 0
len += 1
bi >>= 1
End While
Return len
End Function
<Extension()>
Function BitTest(self As BigInteger, pos As Integer) As Boolean
Dim arr = self.ToByteArray
Dim i = pos \ 8
Dim m = pos Mod 8
If i >= arr.Length Then
Return False
End If
Return (arr(i) And (1 << m)) > 0
End Function
Class PExp
Sub New(p As BigInteger, e As Integer)
Prime = p
Exp = e
End Sub
Public ReadOnly Property Prime As BigInteger
Public ReadOnly Property Exp As Integer
End Class
Function MoBachShallit58(a As BigInteger, n As BigInteger, pf As List(Of PExp)) As BigInteger
Dim n1 = n - 1
Dim mo As BigInteger = 1
For Each pe In pf
Dim y = n1 / BigInteger.Pow(pe.Prime, pe.Exp)
Dim o = 0
Dim x = BigInteger.ModPow(a, y, BigInteger.Abs(n))
While x > 1
x = BigInteger.ModPow(x, pe.Prime, BigInteger.Abs(n))
o += 1
End While
Dim o1 = BigInteger.Pow(pe.Prime, o)
o1 /= BigInteger.GreatestCommonDivisor(mo, o1)
mo *= o1
Next
Return mo
End Function
Function Factor(n As BigInteger) As List(Of PExp)
Dim pf As New List(Of PExp)
Dim nn = n
Dim e = 0
While Not nn.BitTest(e)
e += 1
End While
If e > 0 Then
nn >>= e
pf.Add(New PExp(2, e))
End If
Dim s = nn.Sqrt
Dim d As BigInteger = 3
While nn > 1
If d > s Then
d = nn
End If
e = 0
While True
Dim remainder As New BigInteger
Dim div = BigInteger.DivRem(nn, d, remainder)
If remainder.BitLength > 0 Then
Exit While
End If
nn = div
e += 1
End While
If e > 0 Then
pf.Add(New PExp(d, e))
s = nn.Sqrt
End If
d += 2
End While
Return pf
End Function
Sub MoTest(a As BigInteger, n As BigInteger)
If Not n.IsProbablyPrime(20) Then
Console.WriteLine("Not computed. Modulus must be prime for this algorithm.")
Return
End If
If a.BitLength < 100 Then
Console.Write("ord({0})", a)
Else
Console.Write("ord([big])")
End If
If n.BitLength < 100 Then
Console.Write(" mod {0}", n)
Else
Console.Write(" mod [big]")
End If
Dim mob = MoBachShallit58(a, n, Factor(n - 1))
Console.WriteLine(" = {0}", mob)
End Sub
Sub Main()
MoTest(37, 3343)
MoTest(BigInteger.Pow(10, 100) + 1, 7919)
MoTest(BigInteger.Pow(10, 1000) + 1, 15485863)
MoTest(BigInteger.Pow(10, 10000) - 1, 22801763489)
MoTest(1511678068, 7379191741)
MoTest(3047753288, 2257683301)
End Sub
End Module
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
Wren
import "./big" for BigInt
class PExp {
construct new(prime, exp) {
_prime = prime
_exp = exp
}
prime { _prime }
exp { _exp }
}
var moBachShallit58 = Fn.new { |a, n, pf|
var n1 = n - BigInt.one
var mo = BigInt.one
for (pe in pf) {
var y = n1 / pe.prime.pow(pe.exp)
var o = 0
var x = a.modPow(y, n.abs)
while (x > BigInt.one) {
x = x.modPow(pe.prime, n.abs)
o = o + 1
}
var o1 = BigInt.new(o)
o1 = pe.prime.pow(o1)
o1 = o1 / BigInt.gcd(mo, o1)
mo = mo * o1
}
return mo
}
var factor = Fn.new { |n|
var pf = []
var nn = n.copy()
var e = 0
while (!nn.testBit(e)) e = e + 1
if (e > 0) {
nn = nn >> e
pf.add(PExp.new(BigInt.two, e))
}
var s = nn.isqrt
var d = BigInt.three
while (nn > BigInt.one) {
if (d > s) d = nn
e = 0
while (true) {
var dm = nn.divMod(d)
if (dm[1].bitLength > 0) break
nn = dm[0]
e = e + 1
}
if (e > 0) {
pf.add(PExp.new(d, e))
s = nn.isqrt
}
d = d + BigInt.two
}
return pf
}
var moTest = Fn.new { |a, n|
if (!n.isProbablePrime(10)) {
System.print("Not computed. Modulus must be prime for this algorithm.")
return
}
System.write((a.bitLength < 100) ? "ord(%(a))" : "ord([big])")
System.write((n.bitLength < 100) ? " mod %(n) " : "mod([big])")
var mob = moBachShallit58.call(a, n, factor.call(n - BigInt.one))
System.print("= %(mob)")
}
moTest.call(BigInt.new(37), BigInt.new(3343))
var b = BigInt.ten.pow(100) + BigInt.one
moTest.call(b, BigInt.new(7919))
b = BigInt.ten.pow(1000) + BigInt.one
moTest.call(b, BigInt.new(15485863))
b = BigInt.ten.pow(10000) - BigInt.one
moTest.call(b, BigInt.new(22801763489))
moTest.call(BigInt.new(1511678068), BigInt.new(7379191741))
moTest.call(BigInt.new(3047753288), BigInt.new(2257683301))
- Output:
ord(37) mod 3343 = 1114 ord([big]) mod 7919 = 3959 ord([big]) mod 15485863 = 15485862 ord([big]) mod 22801763489 = 22801763488 ord(1511678068) mod 7379191741 = 614932645 ord(3047753288) mod 2257683301 = 62713425
zkl
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.
var BN =Import("zklBigNum");
var Sieve=Import("sieve");
// factor n into powers of primes
// eg 9090 == 2^1 * 3^2 * 5^1 * 101^1
fcn factor2PP(n){ // lazy factors using lazy primes --> (prime,power) ...
Utils.Generator(fcn(a){
primes:=Utils.Generator(Sieve.postponed_sieve);
foreach p in (primes){
e:=0; while(a%p == 0){ a /= p; e+=1; }
if (e) vm.yield(p,e);
if (a<p*p) break;
}
if (a>1) vm.yield(a,1);
},n)
}
fcn _multOrdr1(a,p,e){
m:=p.pow(e);
t:=m/p*(p - 1);
qs:=L(BN(1));
foreach p2,e2 in (factor2PP(t)){
qs=[[(e,q); [0..e2]; qs; '{ q*BN(p2).pow(e) }]];
}
qs.filter1('wrap(q){ a.powm(q,m)==1 });
}
fcn multiOrder(a,m){
if (m.gcd(a)!=1) throw(Exception.ValueError("Not co-prime"));
res:=BN(1);
foreach p,e in (factor2PP(m)){
res = res.lcm(_multOrdr1(BN(a),BN(p),e));
}
return(res);
}
multiOrder(37,1000).println();
b:=BN(10).pow(20)-1;
multiOrder(2,b).println();
multiOrder(17,b).println();
b=0d10_0001;
[BN(1)..multiOrder(54,b)-1].filter1('wrap(r,b54){b54.powm(r,b)==1},BN(54)) :
if (_) println("Exists a power r < 9090 where (54^r)%b)==1");
else println("Everything checks.");
- Output:
100 3748806900 1499522760 Everything checks.