剰余系における積と累乗,逆数

積

a×b mod mは，BASICでは，MOD(a*b, m)です。
a*bが数値を正確に表現できる上限（10進モードだと16桁強）を超えないように注意します。
十進モードだと，a,b＜1E8であれば大丈夫です。それを超えるときは有理数モードで計算してください。

累乗

an mod mをMOD(a^n,m)で計算すると，a^nが数値を正確に表現できる上限を超える可能性が高くなります。
たとえば，十進モードで 

10 PRINT MOD(9^88,89)
20 END

を実行すると41になりますが，正しくは1です。有理数モードだと正しい答えが得られます。

十進モードで累乗を計算する

法の数が1E8以下のときは，累乗を剰余系における乗算の繰り返しで計算すれば十進モードでも正しい値を求めることができます。
次のプログラムは十進モードでも正しい答えを計算します。

100 DECLARE FUNCTION PowerMod
110 PRINT PowerMod(9,88,89)
120 END
130 EXTERNAL FUNCTION PowerMod(a,n,m)
140 DECLARE NUMERIC x,i
150 LET x=1
160 FOR i=1 TO n
170    LET x=MOD(x*a,m)
180 NEXT i
190 LET PowerMod=x
200 END FUNCTION

しかし，繰り返し回数が多くなると時間がかかります。 p=99999989 として PowerMod(2,p-1,p)　を計算すると十進モードで40秒ほどかかりました。

累乗計算の高速化

a2, a4, a8, a16, …は，2乗することの反復で素早く計算できます。
これを利用して，たとえば，a37=a1+4+32=a1a4a32と変型すれば累乗計算が高速化できます。
この計算法は，特に，剰余系での累乗を計算するのに有効です。次のプログラムはほぼ一瞬で答えを出します。

100 DECLARE FUNCTION PowerMod
110 LET p=99999989 
120 PRINT PowerMod(2,p-1,p)
130 END
140 EXTERNAL FUNCTION PowerMod(a,n,m)
150 DECLARE NUMERIC x,i,y
160 LET x=1
170 LET y=a
180 DO WHILE n>0
190    IF MOD(n,2)=1 THEN LET x=MOD(x*y,m)
200    LET y=MOD(y^2,m)
210    LET n=INT(n/2)
220 LOOP
230 LET PowerMod=x
240 END FUNCTION

さらに桁数の大きい数を法としたいときは有理数モードで計算します。

100 OPTION ARITHMETIC RATIONAL
110 DECLARE FUNCTION PowerMod
120 LET a=2
130 LET p=2^32+1
140 PRINT PowerMod(a,p-1,p)
150 END
800 EXTERNAL FUNCTION PowerMod(a,n,m)
810 OPTION ARITHMETIC RATIONAL
820 DECLARE NUMERIC x,i,y
830 LET x=1
840 LET y=a
850 DO WHILE n>0
860    IF MOD(n,2)=1 THEN LET x=MOD(x*y,m)
870    LET y=MOD(y^2,m)
880    LET n=INT(n/2)
890 LOOP
900 LET PowerMod=x
910 END FUNCTION

p=232+1とすると，a=2のときap-1≡1 (mod p)となりますが，a=3だとap-1≡1 (mod p)とならないので，232+1は素数でないことがわかります（フェルマーの小定理による）。

剰余系における逆数

法mの剰余系におけるaの逆数とは，ax≡1 (mod m) となるxのことです。
これは，ax=1+kmとなる整数kの存在を意味しますから，ax-mk=1となる整数x,kから求まります。 ax-mk=1となるx,kは一次不定方程式ax+by=1の解でk=-yとおいて求まりますが，kを求める必要はないので，一次不定方程式ax+by=1を解いてxを求めればそれが剰余系における逆数です。
次のプログラムの800行で定義される関数InvMod(a,m)は，mを法とする剰余系におけるaの逆数を求めます。存在しないときは，EXTYPE 999の例外を生成します。
外部副プログラムEuclid(a,b,x,y,g)は，g=GCD(a,b)と,ax+by=gとなる整数x,yを求めます。

100 OPTION ARITHMETIC RATIONAL
110 DECLARE EXTERNAL FUNCTION InvMod
120 DECLARE NUMERIC a,p
130 INPUT a,p
140 PRINT InvMod(a,p)
150 END
1000 EXTERNAL FUNCTION InvMod(a,m)
1010 OPTION ARITHMETIC RATIONAL
1020 DECLARE EXTERNAL FUNCTION Euclid
1030 DECLARE NUMERIC x,y,g
1040 CALL Euclid(a,m,x,y,g)
1050 ! ax+my=g, g=GCD(a,m)
1060 IF g<>1 THEN CAUSE EXCEPTION 999 
1070 LET InvMod=MOD(x,m)
1080 END FUNCTION
2000 EXTERNAL SUB Euclid(a,b,x,y,g)
2010 OPTION ARITHMETIC RATIONAL
2020 DECLARE NUMERIC q,r,u,v
2030 LET q=INT(a/b)
2040 LET r=MOD(a,b) 
2050 IF r=0 THEN
2060    LET x=0
2070    LET y=1
2080    LET g=b
2090 ELSE 
2100    CALL Euclid(b,r,u,v,g)
2110    LET x=v
2120    LET y=u-q*v
2130 END IF
2140 END SUB
戻る