modpower

[Roger|Dyalog]

The modpower function, ⍺⍺|⍺*⍵, is known variously as modpow, powermod, modular exponentiation. It is useful in primality testing and public key cryptography.

Model

For small values, modpower can just be computed directly:

      1000|13*7
517
      13*7
62748517

The direct computation fails before too long, when the intermediate values exceed the precision of 64-bit floating point numbers:

      1000|13*21
0
      13*21
2.47065E23

      0⍕13*21
 2470645290734504________

The 0 answer for 1000|13*21 is obviously incorrect because it is saying that 1000 divides evenly into 13*21.

The problem can be overcome by using repeated squaring to do the exponentiation, applying residue every step of the way. Thus:

      ⍝
 z←a(m modpower)n;s
⍝ m|a*n
 z←1
 s←m|a
 :While 0≠n
   :If 1=2|n ⋄ z←m|z×s ⋄ :EndIf
   s←m|×⍨s
   n←⌊n÷2
 :EndWhile

      13 (1000 modpower) 21
413

      1000|13*7
517
      1000|517*3
413

The last two expressions serve as a check because 1000|13*21 ↔ 1000|(13*7)*3 ↔ 1000|(1000|13*7)*3.

Because repeated squaring is efficient, the computation is quick even for large exponents:

      16 (1000 modpower) 1e9
376

Rational modpower

For some applications an extended precision version of modpower is required. This can be done using the rational arithmetic facility, function qi for creating rational arrays and monadic operator Q for working with them.

      ⍝
qmodpower←{   ⍝ ⍺⍺|⍺*⍵
 m←⍺⍺                                  
 z←{⊃⌽⍵:m|Q×Q/⍺ ⋄ 0⌷⍺}                                 
 ((qi 1),m|Q ⍺) {⍬≡⍵:0⌷⍺ ⋄ ((⍺ z ⍵),m|Q×Q⍨1⌷⍺)∇ ¯1↓⍵} ⊤Q ⍵
}

      (qi 13) ((qi 1000) qmodpower) qi 7
┌───────────┐
│┌─────┬─┬─┐│
││5 1 7│1│1││
│└─────┴─┴─┘│
└───────────┘
      1000|13*7
517

      ⊢ t←(qi 2) ((qi 1e12) qmodpower) qi ⊂'1e8000'
┌───────────────────────────┐
│┌─────────────────────┬─┬─┐│
││8 1 7 8 7 1 0 9 3 7 6│1│1││
│└─────────────────────┴─┴─┘│
└───────────────────────────┘
      ⍕Q t
81787109376

(Misalignments in the display are due to defects in the APL Chat Forum software.)

The last example shows that the last 12 decimal digits of 2*1e8000 are 81787109376.

qmodpower uses to advantage ⊤Q ⍵, which produces the binary representation of ⍵. (This was the definition of ⊤⍵ in A Formal Description of System/360 in 1964, but was put aside from APL\360 due to lack of space.)