07-01-2018, 08:19 PM
Modular Arithmetic for the HP-48
We assume that n is the modulus: 0 < n < 1,000,000,000,000
Addition
When we add two integers a and b their sum a + b may overflow when it's bigger than 999,999,999,999.
However if we can assume that both values are between 0 and a modulus n we can calculate a - n + b instead.
This leads to the following program:
Example
This allows to calculate:
500,000,000,000 + 500,000,000,001 = 2 (mod 999,999,999,999)
Otherwise the last digit of the sum is cancelled and the result was 1.
Note:
The result is between 0 and n.
Multiplication
To multiply two numbers we can translate them to the base k = 1,000,000.
a = x * k + y
b = u * k + v
We can then calculate their product:
a * b
= (x * k + y) * (u * k + v)
= x * u * k2
+ x * v * k
+ y * u * k
+ y * v
Note:
Even though e.g. x * u * k2 is bigger than 1'000'000'000'000 the HP-48 calculates MOD n correctly.
We just have to make sure that the product x * u doesn't overflow.
Power
We will use the binary method for modular exponentiation.
Python
Note:
Python provides the function pow(a, b, n) out of the box.
UserRPL
Miller-Rabin Primality Test
If you only want to check if a number n is prime and don't need to find the factors there are more efficient algorithms to do so.
One of them is the Miller-Rabin Primality Test.
Here we use the deterministic variant.
Composite
Check with witness a if n is composite.
If the result is 1 we can be sure that n is composite.
If the result is 0 the number may still be composite.
In this case a is called a strong liar for n.
Example:
n = 221
n - 1 = 220 = 22 55
Thus 174 is a strong liar for 221.
But:
Thus we now can be sure that 221 = 13 * 17 is indeed composite.
Python
UserRPL
Prime?
According to Deterministic variants:
Python
UserRPL
Examples:
In fact: 999999999997 = 5507×181587071
It would be interesting to see how long this takes on a real HP-28.
We can use the division by 0 trick to raise an error since DOERR is missing.
Cheers
Thomas
We assume that n is the modulus: 0 < n < 1,000,000,000,000
Addition
Code:
PLUS ( a b n -- (a + b) % n )
When we add two integers a and b their sum a + b may overflow when it's bigger than 999,999,999,999.
However if we can assume that both values are between 0 and a modulus n we can calculate a - n + b instead.
This leads to the following program:
Code:
« → a b n
« a b +
IF DUP n <
THEN
ELSE DROP a n - b +
END
»
»
Example
This allows to calculate:
500,000,000,000 + 500,000,000,001 = 2 (mod 999,999,999,999)
Otherwise the last digit of the sum is cancelled and the result was 1.
Note:
The result is between 0 and n.
Multiplication
Code:
TIMES ( a b n -- a * b % n )
To multiply two numbers we can translate them to the base k = 1,000,000.
a = x * k + y
b = u * k + v
We can then calculate their product:
a * b
= (x * k + y) * (u * k + v)
= x * u * k2
+ x * v * k
+ y * u * k
+ y * v
Note:
Even though e.g. x * u * k2 is bigger than 1'000'000'000'000 the HP-48 calculates MOD n correctly.
We just have to make sure that the product x * u doesn't overflow.
Code:
« 1000000 → a b n k
« a k / IP a k MOD
b k / IP b k MOD
→ x y u v
« x u * k SQ * n MOD
x v * k * n MOD n PLUS
y u * k * n MOD n PLUS
y v * n MOD n PLUS
»
»
»
Power
Code:
POWER ( a b n -- a ^ b % n )
We will use the binary method for modular exponentiation.
Python
Code:
def power(a, b, n):
r = 1
while b > 0:
if b % 2 == 1:
r = r * a % n
a = a * a % n
b /= 2
return r
Note:
Python provides the function pow(a, b, n) out of the box.
UserRPL
Code:
« → n
« 1 ROT ROT
WHILE DUP
REPEAT → r a b
«
IF b 2 MOD
THEN a r n TIMES
ELSE r
END
a DUP n TIMES
b 2 / IP
»
END
DROP2
»
»
Miller-Rabin Primality Test
If you only want to check if a number n is prime and don't need to find the factors there are more efficient algorithms to do so.
One of them is the Miller-Rabin Primality Test.
Here we use the deterministic variant.
Composite
Check with witness a if n is composite.
If the result is 1 we can be sure that n is composite.
If the result is 0 the number may still be composite.
In this case a is called a strong liar for n.
Example:
n = 221
n - 1 = 220 = 22 55
Code:
174 2 55 221 COMPOSITE? → 0
Thus 174 is a strong liar for 221.
But:
Code:
137 2 55 221 COMPOSITE? → 1
Thus we now can be sure that 221 = 13 * 17 is indeed composite.
Code:
COMPOSITE? ( a s d n -- true if n is composite )
Python
Code:
def composite(a, s, d, n):
x = pow(a, d, n)
if x == 1:
return False
for i in range(s):
if x == n - 1:
return False
x = x * x % n
return True
UserRPL
Code:
« DUP 1 - → a s d n m
« a d n POWER
IF DUP 1 ==
THEN DROP 0
ELSE
IFERR s
WHILE DUP
REPEAT
SWAP
IF DUP m ==
THEN 0 DOERR
END
DUP n TIMES
SWAP 1 -
END
THEN DROP2 0
ELSE DROP2 1
END
END
»
»
Prime?
Code:
PRIME? ( n -- true if n is prime )
According to Deterministic variants:
Quote:if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
Python
Code:
def prime(n):
d, s = n - 1, 0
while not d % 2:
d, s = d >> 1, s + 1
return not any(composite(a, s, d, n) for a in (2, 13, 23, 1662803))
UserRPL
Code:
« → n
« 0 n 1 -
WHILE DUP 2 MOD NOT
REPEAT
SWAP 1 +
SWAP 2 /
END → s d
« { 2 13 23 1662803 }
IFERR
WHILE DUP SIZE
REPEAT
IF DUP HEAD s d n COMPOSITE?
THEN 0 DOERR
END
TAIL
END
THEN DROP 0
ELSE DROP 1
END
»
»
»
Examples:
Code:
999999999997 PRIM?
0
In fact: 999999999997 = 5507×181587071
Code:
999999999989 PRIM?
1
Code:
177777773 PRIM?
1
It would be interesting to see how long this takes on a real HP-28.
We can use the division by 0 trick to raise an error since DOERR is missing.
Cheers
Thomas