PRINT
PRINT "＜フェルマーの２平方の定理＞"
PRINT "     ４で割って１余る素数はすべて、２つの平方数の和で表すことができる。"
PRINT
PRINT "４で割って１余る素数を小さい順に記述します："
PRINT "     ";
DIM W1(1000),W2(1000),W3(1000),H(1000)
LET e=0
FOR i=1 TO 2610
   IF PrimeQ(i)=1 THEN
      IF MOD(i,4)=1 THEN
         LET H(e+1)=i
         PRINT USING "######":i;
         LET e=e+1
         IF INT(e/10)=e/10 THEN
            PRINT
            PRINT "     ";
         END IF
      END IF
   END IF
NEXT i
PRINT
PRINT
PRINT "検索結果："
LET c=1
FOR a=0 TO 50
   FOR b=a TO 50
      LET x=a^2+b^2
      IF PrimeQ(x)=1 THEN
         IF MOD(x,4)=1 THEN
            LET W1(c)=x
            LET W2(c)=a
            LET W3(c)=b
            LET c=c+1
         END IF
      END IF
   NEXT b
NEXT a
CALL sort(W1,W2,W3,c-1)
FOR i=1 TO c-1
   IF W1(i)>2620 THEN GOTO 10
   IF W1(i)<>H(i) THEN
      PRINT "MISS!"
   END IF
   PRINT "    ";W1(i);" ＝ 4 *";INT(W1(i)/4);"＋";MOD(W1(i),4);" ＝ ";W2(i);"^ 2 ＋ ";W3(i);"^ 2" 
10 NEXT i
   PRINT
   PRINT "  以上 少なくとも 2609 までは、"
   PRINT "フェルマーの２平方の定理が成り立っていることがわかりました。"
   PRINT
END

EXTERNAL FUNCTION PrimeQ(n)       ! 素数判定   1 : 素数、  0 : 素数でない
   LET PrimeQ=0
   IF n<2 OR n<>INT(n) THEN EXIT FUNCTION 
   IF MOD(n,2)=0 THEN
      IF n=2 THEN LET PrimeQ=1
   ELSE
      LET k=3
      DO WHILE k*k<=n      ! ３〜√n の奇数のみで検証する
         IF MOD(n,k)=0 THEN EXIT FUNCTION 
         LET k=k+2
      LOOP
      LET PrimeQ=1
   END IF
END FUNCTION

EXTERNAL SUB sort(a(),b(),c(),m)
   FOR i=1 TO m-1
      FOR j=i+1 TO m
         IF a(i)>a(j) THEN
            LET t1=a(j) 
            LET a(j)=a(i)
            LET a(i)=t1
            LET t2=b(j) 
            LET b(j)=b(i)
            LET b(i)=t2
            LET t3=c(j) 
            LET c(j)=c(i)
            LET c(i)=t3
         END IF
      NEXT j
   NEXT i
END SUB