Solving a Single Congruence Equation

[Thomas Klemm]

This program for the HP-42S solves for x in the equation:

A * x = B mod N

Example:
5 * x = 3 mod 17

Solution:
5 ENTER 3 ENTER 17
XEQ "CONG"
4



Remark:
Since overflow is avoided it's possible to solve:
999,999,999,989 * x = 1 mod 999,999,999,999

Solution:
899,999,999,999

[Gerald H]

If you wanted to incorporate this version of BEZO in your programme you would probably see a reduction in the time required for the calculation.

This consideration is really only valid for the real calculator, the difference would not be noticeable on, eg, Free42.

0. { 75-Byte Prgm }
1. LBL “BEZO”
2. STO 08
3. 1
4. COMPLEX
5. X<>Y
6. STO 07
7. 0
8. COMPLEX
9. LBL 00
10. ENTER
11. COMPLEX
12. R↓
13. X=0?
14. GTO 01
15. RCL ST Z
16. X<>Y
17. /
18. COMPLEX
19. R↓
20. IP
21. RCL* ST Y
22. RCL- ST Z
23. +/-
24. GTO 00
25. LBL 01
26. RCL ST Z
27. COMPLEX
28. RCL ST Y
29. RCL ST Y
30. RCL* 08
31. –
32. RCL/ 07
33. STO 05
34. COMPLEX
35. RCL ST Y
36. SIGN
37. STO* 05
38. STO* ST Z
39. *
40. RCL 08
41. RCL 07
42. COMPLEX
43. END

[Ángel Martin]

(04-18-2014 05:24 PM)Thomas Klemm Wrote:  This program for the HP-42S solves for x in the equation:

A * x = B mod N

Example:
5 * x = 3 mod 17

Solution:
[font=Courier]5 ENTER 3 ENTER 17
XEQ "CONG"
4

Great example Thomas, thanks for sharing it.
I have used it to test the RCL Math functions from the Total_Rekall module, simply changing a few lines make it also applicable to the 41 platform now!

Cheers,
ÁM

[Thomas Klemm]

(04-22-2014 07:45 AM)Thomas Klemm Wrote:  

(04-22-2014 06:03 AM)Les Bell Wrote:  Now, about your "Solving a Single Congruence Equation" listing . . .

I know. That won't be easy without RCL-arithmetic.

(07-24-2015 12:23 PM)Ángel Martin Wrote:  I have used it to test the RCL Math functions from the Total_Rekall module, simply chaining a few lines make it also applicable to the 41 platform now!

Thanks for making it possible!

Cheers
Thomas

[Bunuel66]

I apologize for being too lazy for analyzing the algorithm given in the topic. Then please, forgive me if the following discussion is redundant.

Strictly from an algorithmic standpoint (I haven't checked in terms of speed), it could be more straightforward to use Fermat's little theorem who says in a nutshell that if A is prime relatively to N then A^(N-1)=1 Mod N which is equivalent to say that A^(N-2) is an inverse of A Mod N.
Then, as the equation is unchanged if we consider A Mod N in case A>N, the solution is : x=B*A^(N-2) Mod N.
As A^(N-2) can be computed Mod N at every step the risk of overflow is reduced.
In addition, for reducing the number of times that the modulo is computed one can go through a test on the calculation of A^n and compute the modulo only when the product is near the accuracy limit of the machine.

Example: with the given equation 5*x=3 Mod 17 one gets x=3*5^(15) Mod 17
5^15 Mod 17=7 and 3*7 Mod 17=4

There is another approach using the extended Euclid Algorithm and the Bezout/Merignac theorem but it will be more complex to program. This alternate approach will be the faster at the price of an expense in memory and the use of a structure array like ;-(

My 2 cents...

[Thomas Klemm]

(07-25-2015 12:54 PM)Bunuel66 Wrote:  Strictly from an algorithmic standpoint (I haven't checked in terms of speed), it could be more straightforward to use Fermat's little theorem who says in a nutshell that if A is prime relatively to N then A^(N-1)=1 Mod N which is equivalent to say that A^(N-2) is an inverse of A Mod N.

Fermat's little theorem:
If p is a prime number and a is not divisible by p then: \(a^{p-1} \equiv 1 \pmod p.\)

Euler's theorem:
If n and a are coprime positive integers, then: \(a^{\varphi (n)} \equiv 1 \pmod{n}\)

Thus you'd first had to calculate \(\varphi (n)\).

For the given example this is: \(\varphi (999999999999)=461894400000\)

Quote:As A^(N-2) can be computed Mod N at every step the risk of overflow is reduced.

How do you intend to calculate 999999999989^461894399999 mod 999999999999 without overflow?

Kind regards
Thomas

[Gerald H]

Modulo powering on 42S:

a ^ p modulo m

Stack

Z: Integer to power, a
Y: Integer power, p
X: Integer modulus, m

Actuate MOD↑

& the answer is returned to the stack.

Code:


0.    { 42-Byte Prgm }
1.    LBL “SQM”
2.    STO ST Y
3.    1E6
4.    MOD
5.    STO ST Z
6.    –
7.    ENTER
8.    X^2
9.    RCL 00
10.    MOD
11.    X<>Y
12.    R↑
13.    STO* ST T
14.    *
15.    RCL 01
16.    MOD
17.    RCL- ST L
18.    RCL+ ST L
19.    RCL 01
20.    MOD
21.    +
22.    RCL 00
23.    MOD
24.    +
25.    RCL 01
26.    MOD
27.    END




0.    { 26-Byte Prgm }
1.    LBL “INVM”
2.    XEQ “BEZO”
3.    RCL ST Z
4.    DSE ST X
5.    GTO 00
6.    RCL 05
7.    RCL 08
8.    MOD
9.    RTN
10.    LBL 00
11.    CLX
12.    END


0.    { 77-Byte Prgm }
1.    LBL “M/”
2.    X<>Y
3.    STO 02
4.    R↓
5.    RCL 01
6.    XEQ “INVM”
7.    X=0?
8.    RTN
9.    RCL 02
10.    LBL “M*”
11.    RCL ST X
12.    1E6
13.    MOD
14.    X<>Y
15.    RCL- ST Y
16.    COMPLEX
17.    X<>Y
18.    1E6
19.    MOD
20.    RCL ST Z
21.    RCL- ST Y
22.    RCL* ST Z
23.    X<>Y
24.    RCL* ST Z
25.    COMPLEX
26.    RCL 00
27.    MOD
28.    +
29.    RCL 01
30.    MOD
31.    X<>Y
32.    COMPLEX
33.    RCL 00
34.    MOD
35.    X<>Y
36.    RCL 01
37.    MOD
38.    +
39.    RCL 00
40.    MOD
41.    +
42.    RCL 01
43.    MOD
44.    END


0.    { 60-Byte Prgm }
1.    LBL “MOD↑”
2.    STO 01
3.    +/-
4.    STO 00
5.    R↓
6.    LBL “M↑”
7.    STO 02
8.    R↓
9.    STO 03
10.    SIGN
11.    GTO 00
12.    LBL 01
13.    2
14.    MOD
15.    X≠0?
16.    GTO 02
17.    LASTX
18.    STO/ 02
19.    RCL 03
20.    XEQ “SQM”
21.    STO 03
22.    RCL 02
23.    GTO 01
24.    LBL 02
25.    STO- 02
26.    RCL 04
27.    RCL 03
28.    XEQ “M*”
29.    LBL 00
30.    STO 04
31.    RCL 02
32.    X≠0?
33.    GTO 01
34.    R↓
35.    END

[Bunuel66]

Quote:Fermat's little theorem:
If p is a prime number and a is not divisible by p then: \(a^{p-1} \equiv 1 \pmod p.\)

Yes you're right with the restriction you mention, but in that case it seems to be possible to use Bezout/Merignac for computing the inverse:

1/ if A and N are coprimes, then it exist u,v as: A.u+N.v=1 then u=\(A^{-1}\) Mod N
Ax=B Mod N => uAx=uB Mod N =>x=uB Mod N

2/ if A and N are not coprimes it exists C dividing A and N then let A=A/C, N=N/C and B=B/C

2.a if C doesn't divide B there is no solution

2.b if C divides B we are back to 1 as A, N are now coprimes.
x is solution of Ax=B mod N

Examples:
5x=3 Mod 17 Using Euclide extended one finds 7*5-2*17=1 then \(5^{-1}\)=7 Mod 17
x=3*7 Mod 17=4 Mod 17

12x=6 Mod 15 according to the above discussion: => 4x=2 Mod 5 and \(4^{-1}\)=4 Mod 5
then x=4*2 Mod 5 => x=3 Mod 5

Quote:Euler's theorem:
If n and a are coprime positive integers, then: \(a^{\varphi (n)} \equiv 1 \pmod{n}\)

Thus you'd first had to calculate \(\varphi (n)\).

For the given example this is: \(\varphi (999999999999)=461894400000\)

Well, it seems that \(\varphi (n)\) is not needed in the above approach. That said, I'm ready to admit that the extended Euclid Algorithm applied directly will probably lead to computations of that complexity. But it is certainly possible to conduct the computations in modular form for limiting the expansion of numbers.

Quote:How do you intend to calculate 999999999989^461894399999 mod 999999999999 without overflow?

It is mainly a question of what you have in hand. If you are able to compute at least the square of your number without overflow then the modulo operation will limit the expansion at every step
of the loop. It is not that difficult to write algorithms for the 4 elementary operations working on 24 figures if you have operators working on a lower number. It is a bit a brute force solution, but is should work.

There is certainly a more astute approach?

My 3 cents....

[Thomas Klemm]

There was a bug in my program. I had to replace the following lines:

Code:

41>LBL 03
42 2
43 ÷
44 ENTER
45 IP
46 STO 07
47 -
48 X=0?

These are the corrected lines:

Code:

41>LBL 03
42 2
43 MOD
44 STO- 07
45 LASTX
46 STO÷ 07
47 X<>Y
48 X=0?

The problem occurred only when the value in register 07 used 12 digits.
E.g. 937,499,999,999 ÷ 2 would yield 468,750,000,000 instead of 468,749,999,999.5 and thus got rounded up.
It appears that I've tested the example with Free42 where this doesn't happen.

Sorry for any inconveniences.
Thomas

[Thomas Klemm]

(07-25-2015 12:54 PM)Bunuel66 Wrote:  I apologize for being too lazy for analyzing the algorithm given in the topic.

No reason to apologize. I should have explained the algorithm or documented the program.
But sometimes I'm too lazy.

Quote:There is another approach using the extended Euclid Algorithm and the Bezout/Merignac theorem but it will be more complex to program. This alternate approach will be the faster at the price of an expense in memory and the use of a structure array like ;-(

Guess what: the extended Euclid Algorithm is used in my program.

(07-25-2015 09:04 PM)Bunuel66 Wrote:  It is mainly a question of what you have in hand. If you are able to compute at least the square of your number without overflow then the modulo operation will limit the expansion at every step
of the loop. It is not that difficult to write algorithms for the 4 elementary operations working on 24 figures if you have operators working on a lower number. It is a bit a brute force solution, but is should work.

For this approach you may have a look at Gerald's post. However I tried to avoid splitting numbers:

Code:

2.    STO ST Y
3.    1E6
4.    MOD
5.    STO ST Z
6.    –

Instead I implemented the multiplication based on addition and made sure that an overflow can't happen:

Code:

65>LBL 00      ; add ( y x -- y+x )
66 RCL 00      ; y x n
67 RCL ST Z    ; y x n y
68 RCL+ ST Z   ; y x n y+x
69 X<Y?        ; no overflow
70 RTN         ; y+x
71 Rv          ; y x n
72 -           ; y x-n
73 +           ; y+x-n
74 END         ; y+x MOD n

You might consider this approach a thought experiment.

If overflow isn't an issue you can resort to the solution in this post.
I assume it's faster than using Eddie's brute force approach but I haven't done any measurements.

Kind regards
Thomas