剰余の計算

 投稿者:GAIメール  投稿日:2008年12月24日(水)12時02分1秒
  どんな大きな(10桁ほどは欲しい)x,n,aに対しても、次の剰余が求まるプログラム
を作って頂きたいです。
x^n mod a
 

Re: 剰余の計算

 投稿者:山中和義  投稿日:2008年12月24日(水)12時52分41秒
  > No.190[元記事へ]

GAIさんへのお返事です。

> どんな大きな(10桁ほどは欲しい)x,n,aに対しても、次の剰余が求まるプログラム
> を作って頂きたいです。
> x^n mod a


・提供済みUBASIC関数のmodpow(1234,56,789)


・PRINT MOD(1234^56,789) !1234^56 mod 789
 多桁の整数なら有理数モードで実行します。
 x^nが大きいと時間がかかります。


nが負や実数なら検討しないといけません。
 

Re: 剰余の計算

 投稿者:荒田浩二  投稿日:2008年12月25日(木)00時05分19秒
  > No.190[元記事へ]

GAIさんへのお返事です。

> どんな大きな(10桁ほどは欲しい)x,n,aに対しても、次の剰余が求まるプログラム
> を作って頂きたいです。
> x^n mod a

  (x mod a) = r とすれば、
  (x^n mod a) = (r^n mod a) といえるので、
   x>a のときは数値が大きくなりすぎず効果があります。

  [注意] 十進BASICのヘルプによると、有理数モードでは
  「べき指数は-2147483647~2147483647の範囲の整数に限定」
  されるそうです。

! mod(x^n,a) 有利数モードで実行
LET x=1234
LET n=56
LET a=789
PRINT "x=";x
PRINT "n=";n
PRINT "a=";a
PRINT
PRINT "x^n=";x^n
PRINT "mod(x^n,a)=";MOD(x^n,a)
PRINT
LET r=MOD(x,a)
PRINT "r=mod(x,a)=";MOD(x,a)
PRINT "r^n=";r^n
PRINT "mod(r^n,a)=";MOD(r^n,a)
END
 

Re: 剰余の計算

 投稿者:荒田浩二  投稿日:2008年12月25日(木)02時09分38秒
  > No.190[元記事へ]

GAIさんへのお返事です。

> どんな大きな(10桁ほどは欲しい)x,n,aに対しても、次の剰余が求まるプログラム
> を作って頂きたいです。
> x^n mod a

  下のプログラムなら、計算中の最大値はr^2になります。
  rが8桁以下なら10進モードや2進モードで実行できます。
  rが8桁を超える場合は1000桁モードになります。

  x,n,aがともに10桁の下の例では、1000桁モードで約40分かかりました。

! mod(x^n,a)
LET x=9876543201
LET n=1234567890
LET a=6789012345

LET r=MOD(x,a)
LET k=r
FOR i=2 TO n
   LET k=MOD(k*r,a)
NEXT i
PRINT k
END
 

Re: 剰余の計算

 投稿者:SECOND  投稿日:2008年12月25日(木)04時29分36秒
  > No.190[元記事へ]

GAIさんへのお返事です。

山中氏の、EXTERNAL FUNCTION modpow(a,b,n) !a^b mod n ですが、
13桁まで、UBASIC の実行結果と一致します。(十進15桁defaultのまま)
14桁は、一致しない。
   10   print modpow(9999999999990,9999999999991,9999999999992)
   20   print modpow(9999999999991,9999999999992,9999999999993)
   30   print modpow(9999999999992,9999999999993,9999999999994)
   40   print modpow(9999999999993,9999999999994,9999999999995)
   50   print modpow(1999999999999,2999999999999,3999999999999)
   60   print modpow(2999999999999,3999999999999,4999999999999)
   70   print modpow(3999999999999,4999999999999,5999999999999)
   80   print modpow(4999999999999,5999999999999,6999999999999)
run
4328812006896
6296396057422
709006245058
6925339067979
2162020581010
1769850372477
5139725946770
6768106064755
OK
 

Re: 剰余の計算

 投稿者:白石 和夫  投稿日:2008年12月25日(木)09時11分50秒
  > No.190[元記事へ]

「基本のアルゴリズム」のページを見てください。

このプログラムは,a^n mod k を計算します。
kの値が8桁を超えるときは有理数モードで実行してください。

100 INPUT k
110 INPUT a,n
120 LET p=1
130 LET b=MOD(a,k)
140 DO UNTIL n=0
150    IF MOD(n,2)=1 THEN LET p=MOD(p*b,k)
160    LET b=MOD(b*b,k)
170    LET n=INT(n/2)
180 LOOP
190 PRINT p
200 END
 

Re: 剰余の計算

 投稿者:荒田浩二  投稿日:2008年12月25日(木)11時23分39秒
  > No.190[元記事へ]

GAIさんへのお返事です。

  x^nでべき指数をn=s*t*uと分解すれば、剰余計算の回数がs+t+uに減らせるのでnを素因数分解しました。
nが素数でなければ高速です。1000桁モードでも数秒です。

[投稿しようと掲示板を開いたら白石先生の投稿がありました。白石先生のプログラムではnが素数であるかに無関係で高速です。まあ、せっかく作ったので投稿します。]

追加編集[SECONDさんの投稿でさかのぼって確認しましたが、山中和義さんがすでに白石先生と同様の投稿をされていました。]

DECLARE EXTERNAL SUB prime_factor
PUBLIC NUMERIC prn(2000000),pf(100),f
DECLARE FUNCTION modn
LET x=9876543201
LET n=1234567890
LET a=6789012345
LET r=MOD(x,a)
CALL prime_factor(n) ! 素因数分解
LET k=r
FOR i=1 TO f
   LET k=modn(k,pf(i),a)
NEXT i
PRINT k
!
FUNCTION modn(r,n,a) ! mod(r^n,a)
   LET kk=r
   FOR ii=2 TO n
      LET kk=MOD(kk*r,a)
   NEXT ii
   LET modn=kk
END FUNCTION
END

REM ** 素因数分解 **
EXTERNAL SUB prime_factor(n)
DECLARE EXTERNAL SUB prime
LET f=0 ! 素因数の個数
LET q=n ! qは、商
LET m=1 ! 検証用変数
FOR i=1 TO 10
   READ prn(i)
NEXT i
DATA 2,3,5,7,11,13,17,19,23,29
LET i=1
DO WHILE MOD(q,prn(i))=0 ! 因数判定ルーチン
   LET f=f+1
   LET pf(f)=prn(i) ! pf(f)はnのf番目の素因数
   LET q=q/prn(i)
   LET m=m*prn(i)
LOOP
LET i=2
DO WHILE prn(i)<=SQR(q)
   DO WHILE MOD(q,prn(i))=0 ! 因数判定ルーチン
      LET f=f+1
      LET pf(f)=prn(i) ! pf(f)はnのf番目の素因数
      LET q=q/prn(i)
      LET m=m*prn(i)
   LOOP
   IF q=1 THEN EXIT DO
   LET i=i+1
   IF i>10 THEN CALL prime(i) ! i番目の素数
LOOP
IF q<>1 THEN
   LET f=f+1
   LET pf(f)=q
   LET m=m*q
END IF
IF m<>n THEN PRINT "error !!"
END SUB
!
REM ** 素数列生成(k番目の素数) **
EXTERNAL SUB prime(k)
LET m30=MOD(prn(k-1),30)
IF m30=1 OR m30=23 THEN LET a=prn(k-1)+6 ELSE LET a=prn(k-1)-2*MOD(m30,3)+6
DO
   LET sqra=SQR(a)
   FOR j=4 TO k-1
      IF MOD(a,prn(j))=0 THEN EXIT FOR
      IF prn(j)>=sqra THEN
         LET  prn(k)=a
         EXIT SUB
      END IF
   NEXT j
   LET m30=MOD(a,30)
   IF m30=1 OR m30=23 THEN LET a=a+6 ELSE LET a=a-2*MOD(m30,3)+6
LOOP
END SUB
 

戻る