    1MODE 12
    2FOR d = 1 TO 2
    3PRINT'' "bbc basic v/64 prog ""calc.factors.LgeFactor.powrMOd2"",  24.05.98."
    4PRINT "n.b.  * basic64.    ds mcdonald, wellington, 20/5/99.  14:25 pm."CHR$7
    5IF d=1 THEN INPUT ''"spool to *ram. moddata,   1 = yes, 0 = no.";sp
    6IF sp = 1 THEN
    7PRINT "yes   CONT."GET$
    8*ram
    9*sp. moddata
   10ENDIF
   11NEXT d
   12INPUT ' "enter base [retn 2], a = ";a$      "   reduced printout/  index  "P%
   13n = P%   :  NOW = TIME
   14IF a$ = "" THEN a =2   ELSE  a = EVAL a$
   15
   16PRINT  a$  "  =  "; a
   17ON ERROR REPORT : PRINT "   at erl.  "ERL :  IF ERR=17 THEN PRINT"close #0"CHR$7 : CLOSE#0  : STOP
   18PRINT  CHR$7
   19REM INPUT LINE ' "  enter index, n = [ retn +/ next prime  q quit,  r run.] ";n$
   20
   21REM IF n$ = "" OR INSTR("+/",n$) THEN
   22REPEAT
   23   n +=1
   24   IF n MOD 2 = 0 THEN n +=1
   25   FOR i = 3 TO SQR n STEP 2
   26   IF n MOD i  THEN NEXT i  :  i = n
   27   UNTIL i = n
   28   PRINT '' "next prime = "; n;
   29REM    ELSE
   30
   31REM  IF INSTR("qQ", n$) THEN  PRINT P% : END
   32REM  IF INSTR("rR", n$) THEN RUN
   33REM  n = EVAL n$
   34REM  ENDIF
   35
   36REM  PRINT  n$ "  =  "; n
   37REM  INPUTLINE "enter test factors,"'"  numeric = sgle modulus, p = "'"b = 2 ~ 200."'"r retn = range ";p$
   38REM      CASE p$ OF
   39REM     REM IF (INSTR("bB", p$) TH.
   40REM  
   41REM  WHEN "b", "B"
   42REM   FOR p = 2 TO 200
   43REM           test = FNPower(a, n, p)
   44REM             PRINT "test, p = "; p "  reslt = "; test
   45REM             IF test=1 THEN PRINT  "Yes."GET$
   46REM             NEXT p
   47REM  
   48REM  WHEN "r", "R" ,  ""
   49REM   INPUT "start [retn 1]= ", s
   50REM   IF s<1 THEN s=1
   51REM   INPUT "step =  (  [retn  k] = 2kp+1 ) ", st$
   52REM   IF st$="" OR INSTR("kK", st$) THEN st = 2*n ELSE st = EVAL st$
   53REM   PRINT st
   54REM   INPUT "minimum = [retn 1]  (BUG V. slow.) ", min
   55REM   IF min< 1 THEN min =1
   56
   57s=1  :  st=2*n  :  min=1
   58   REM  INPUT ' "test condition  for pause..cont.   " '' "test = TRUE"'"test = -1"'"ABS test < 20"' "e.g.  test^2 = 1.    ", test$
   59   test$ = "test = 1 "  :  REM  PRINT test$
   60   test = 1  :  REM  0
   61REM   INPUT "only Coprime 2,3,5,7,11,13.  Yy"; c$
   62REM      IF INSTR("Yy", c$) THEN          PRINT  "Yes."
   63   NOW = TIME
   64
   65   FOR p = s TO 2^30 STEP st
   66      IF p < min THEN NEXT p :
   67      REM p = p + st*INT((min-p)/st  ) :P."cont."p GET$ :NEXT p
   68
   69      IF INSTR("Yy", "y" ) THEN
   70         IF p MOD 2 = 0  THEN ENDIF : NEXT p
   71         FOR q = 3 TO 13 STEP 2
   72            IF p MOD q = 0  THEN q = 17
   73            NEXT q
   74         IF q = 17 THEN ENDIF  : NEXT p
   75         
   76
   77      test = FNPower(a, n, p)
   78      IF EVAL( test$)  THEN  P% = n:   PRINT  ' ; "k = ";(p-1)/n  " . LOG = ";INT(10*LOG p)" ,  ";p " divides  ";a " ^ "; n "  -1  ..  csec = "TIME-NOW : ERROR :   ELSE PRINT ","; :  IF POS > 66 THEN PRINT ';p "  csec= ";TIME-NOW"  ";
   79      test = 0:  REM clear flag
   80      ENDIF
   81      NEXT p
   82OTHERWISE
   83REM "s", "S"
   84p = EVAL p$
   85PRINT  p$ "  =  "; p
   86
   87PRINT ; "result = "; FNPower(a, n, p)
   88
   89   OTHERWISE
   90
   91   ENDCASE
   92UNTIL FALSE
   93END
   94REM (c) Now use a good algorithm to compute a power modulus 
   95REM  for example :
   96    DEF FNPower( a, n1, p)
   97       r = 1
   98       x = a
   99       WHILE ( n1>0)
  100    IF n1 MOD 2 THEN
  101            r = FNmod(r*x,p)
  102            REM  :  PRINT r
  103            ENDIF
  104              x = FNmod(x*x,p)
  105              REM x = x - p*INT(x/p+.5)
  106              REM  MOD p
  107               n1 = n1 DIV 2
  108    ENDWHILE :  REM           }
  109    REM    return(r);
  110   = r
  111
  112DEF FNmod(z,y)
  113IF z = z+1  THEN ERROR
  114 = z - y*INT(z/ y+.5)
  115       END
  116 
  117   REM  *> Good Luck...
  118   REM   
  119   REM   -- 
  120   REM   Philippe Langevin,
  121   REM   http://www.univ-tln.fr/~langevin/  */