素数の振りして素数に非ず

*ST
INPUT A
L=SQR(A)
L=INT(L)
T=A
FOR N=2 TO L
	IF T MOD N=0 THEN T=T/N:PRINT N
	IF T=1 THEN GOTO *ST
NEXT N
PRINT T
GOTO *ST

*LP1
INPUT A
*LP2
GOSUB *ST
IF T=1 OR A=T THEN PRINT:GOTO *LP1
A=T
GOTO *LP2

*ST
L=SQR(A)
L=INT(L)
T=A
FOR N=2 TO L
	IF T MOD N=0 THEN T=T/N:PRINT N;:RETURN
NEXT N
IF T=A THEN PRINT T;
RETURN

OPTION ARITHMETIC DECIMAL_HIGH

10 INPUT A
20 GOSUB 100
   IF T=1 OR A=T THEN 
      PRINT
      GOTO 10
   END IF
   LET A=T
   GOTO 20
    
100 LET L=SQR(A)
    LET L=INT(L)
    LET T=A
    FOR N=2 TO L
       IF MOD(T , N) =0 THEN
          LET T=T/N
          PRINT N;
          RETURN
       END IF
    NEXT N
    IF T=A THEN PRINT T
    RETURN
 END

    OPTION ARITHMETIC RATIONAL

    Y=100
    FOR X=1 TO 30
        Y=Y*10
        A=Y+1
        PRINT A;
        GOSUB 20
    NEXT X
    STOP

20 GOSUB 100
   IF T=1 OR A=T THEN 
      RETURN
   END IF
   LET A=T
   GOTO 20
    
100
    LET L=INTSQR(A)
    LET T=A
    FOR N=2 TO L
       IF MOD(T , N) =0 THEN
          LET T=T/N
          PRINT N;
          RETURN
       END IF
    NEXT N
    IF T=A THEN PRINT T
    RETURN
 END

    OPTION ARITHMETIC RATIONAL

    Y=100
    FOR X=4 TO 1000
        Y=Y*10
        A=Y+1
        PRINT X;
        IF X<10 THEN PRINT A;
        GOSUB 100
    NEXT X
    STOP

100
    LET L=INTSQR(A)
    LET T=A
    FOR N=2 TO L
       IF MOD(T , N) =0 THEN
          PRINT N
          RETURN
       END IF
    NEXT N
    PRINT "S"
    RETURN
 END

    OPTION ARITHMETIC RATIONAL

    INPUT PROMPT "START ":S
    INPUT PROMPT "ADD   ":A
    LET E=S+A
    FOR X=S TO E
        GOSUB 100
    NEXT X
    PRINT
    STOP

100
    LET L=INTSQR(X)
    FOR N=2 TO L
       IF MOD(X , N) =0 THEN RETURN
    NEXT N
    PRINT X;
    RETURN
 END

OPTION ARITHMETIC RATIONAL

10 INPUT num
   LET n = num
20!'
   CALL factor(n, f)
   IF f = 1 THEN!' Call by reference
      PRINT num; " は素数です"
      GOTO 10
   ELSE
      PRINT num; " は合成数です（"; f; "で割り切れて商は"; n; "）"
      LET n = num
      CALL adumaQ(n)
      GOTO 10
   END IF
END

EXTERNAL SUB factor(n, f)
   OPTION ARITHMETIC RATIONAL
   LET m=INTSQR(n)
   FOR i=2 TO m
      IF MOD(n, i) = 0 THEN
         LET f = i
         LET n = n / i
         EXIT SUB
      END IF
   NEXT i
   LET f = 1
END SUB
!'
EXTERNAL SUB adumaQ(n)
   OPTION ARITHMETIC RATIONAL
   LET N_NEIGHBOR = 3
   IF n <= 2 THEN
      EXIT SUB
   END IF
   LET d = 0
   PRINT " 隣接する素数は近い順に（大小"; N_NEIGHBOR; "つまで）"
   PRINT " 大きな方は"
   FOR i = 1 TO N_NEIGHBOR
      LET f = n
      DO WHILE f > 1
         LET d = d + 1
         CALL factor(n+d, f)
      LOOP
      PRINT "  "; n+d
   NEXT I 
   LET d = 0
   PRINT " 小さな方は"
   FOR i = 1 TO N_NEIGHBOR
      LET f = n
      DO WHILE f > 1
         LET d = d + 1
         CALL factor(n-d, f)
      LOOP
      IF n-d >= 2 THEN
         PRINT "  "; n-d
      ELSE
         PRINT " もうありません"
      END IF
   NEXT I 
END SUB

? 1001
 1001  は合成数です（ 7 で割り切れて商は 143 ）
 隣接する素数は近い順に（大小 3 つまで）
 大きな方は
   1009
   1013
   1019
 小さな方は
   997
   991
   983
? 10001
 10001  は合成数です（ 73 で割り切れて商は 137 ）
 隣接する素数は近い順に（大小 3 つまで）
 大きな方は
   10007
   10009
   10037
 小さな方は
   9973
   9967
   9949
? 100001
 100001  は合成数です（ 11 で割り切れて商は 9091 ）
 隣接する素数は近い順に（大小 3 つまで）
 大きな方は
   100003
   100019
   100043
 小さな方は
   99991
   99989
   99971

OPTION ARITHMETIC RATIONAL

LET n_term = 18

DIM n(20)
FOR i = 1 TO 20
   LET n(i) = 10^i +1
NEXT I 

LET ts = TIME
FOR i = 1 TO n_term
   PRINT i; ")"; n(i); " = ";
   CALL factorize(n(i))
   PRINT
NEXT I
LET te = TIME
PRINT "Elapsed time: "; te - ts; "sec"
END

EXTERNAL SUB factorize(n)
OPTION ARITHMETIC RATIONAL
LET f = 0
DO
   CALL factor(f, n)
   PRINT f;
   IF n = 1 THEN
      EXIT DO
   END IF
LOOP
END SUB

EXTERNAL SUB factor(f, n)
OPTION ARITHMETIC RATIONAL
LET m=SQR(n)
FOR f=2 TO m
   IF MOD(n, f) =0 THEN
      LET n = n/f
      EXIT SUB
   END IF
NEXT f

LET f = n
LET n = 1
END SUB

var TBI1 = new TinyBigInt(1);
var TBI2 = new TinyBigInt(2);
var N_TERM = 18;

var a = new Array();
for (var i = 0; i < N_TERM; i++){
  for (var j = 0, z = ""; j < i; j++){
    z = z + "0";
  }
  a[i] = new TinyBigInt("1" + z + "1");
}

var ts = Date.now();
var i = 0;
for (var ae of a){
  i++;
  var msg = i +") "+ ae +" = ";
  factorize(ae);
  print(msg);
}
var te = Date.now();
print("Ealpsed time: "+(te - ts)/1000 +" sec.");

function factorize(ae){
  do {
    var r = factor(ae);
    msg = msg +" "+ r.f;
    ae = r.n;
  } while( r.n.compare(TBI1));
}

function factor(n){// n: TynyBigInt
  var m = n.sqrt();
  var i = TBI2;
  while(i.compare(m) <= 0){
 　 if (!n.mod(i).compare(0)){
      return {f: i, n: n.div(i)};
    }
    i = i.add(1);
  }
  return {f: n, n: TBI1};
}

#include <stdio.h>
#include <sys/time.h>
#include <gmp.h>

#define N_TERM 18

extern void factorize(mpz_t n);
extern void factor(mpz_t f, mpz_t n);
mpz_t C0, C1, C2;

int main()
{
  unsigned long a = 1001;
  mpz_t mp[N_TERM];
  char mp_s[1000];
  struct timeval ts, te;

  mpz_init_set_str(C0, "0", 10);
  mpz_init_set_str(C1, "1", 10);
  mpz_init_set_str(C2, "2", 10);

  mp_s[0] = '1';
  for (int i = 0; i < N_TERM; ++i){
    for (int j = 1; j < i+1; ++j){ mp_s[j] = '0';}
    mp_s[i+1] = '1';
    mp_s[i+2] = '\0';
    mpz_init_set_str(mp[i], mp_s, 10);
  }

  gettimeofday(&ts, NULL);
  for (int i = 0; i < N_TERM; ++i){
    printf("%03d) ", i+1);
    mpz_out_str(NULL, 10, mp[i]);
    printf(" = ");
    factorize(mp[i]);
    printf("\n");
  }
  gettimeofday(&te, NULL);
  printf("Elapsed time: %f sec.\n", (te.tv_sec - ts.tv_sec) + (te.tv_usec - ts.tv_usec)/1000000.0);
}

void factorize(mpz_t n)
{
  mpz_t f;

  mpz_init(f);
  do {
    factor(f, n);
    mpz_out_str(NULL, 10, f);
    printf(" ");
  } while (mpz_cmp(n, C1));
}

void factor(mpz_t f, mpz_t n)
{
  mpz_t ub, q, r; 
  mpz_t i, ii;
  mpz_init(q);
  mpz_init(r);
  mpz_init(ub);
  mpz_init(ii);
  mpz_sqrt(ub, n);

  mpz_set(i, C2);
  while (mpz_cmp(i, ub) <= 0){
    mpz_fdiv_qr(q, r, n, i);
    if (!mpz_cmp(r, C0)){
      mpz_set(f, i);
      mpz_set(n, q);
      return;
    }
    mpz_add(ii, i , C1);
    mpz_set(i, ii);
  }
  mpz_set(f, n);
  mpz_set(n, C1);
}