リュカ・レーマー・テスト最終版netwalker(ARM)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define N 10000
#define A 800
#define B 1000
#define TXT "lucas7a.txt"

/*
技術評論社　奥村晴彦　著
C言語による最新アルゴリズム事典　２３３０円＋税
P160〜161より引用・改変

A=始める数字
B=計算をここまでの数字までやる
TXT=実行結果を書き込むファイル名

AからBまでの素数のみリュカ・レーマー・テストをする
メルセンヌ素数は、2,3,5,7,13,17,19,31,61,89,107,127,
521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,
19937,21701,23209,44497,86243,110503,132049,216091


リュカ・レーマー・テストのアルゴリズムは

M(P)=2^P-1;　として、

x(1)=4;
x(i+1)=(x(i)^2-2) mod M(P)
x(P-1)=0　なら素数


*/

int prime(int p);
double gettime();
double st,pt;

main(){
   int p,a=A,b=B,c,m,s,i;
   int d,e,f,g,h;
   double da,db,dc,dd;

   char s1[15]="";
   char s2[ 5]="a";
   char s3[10]="fail\n";
   char s4[25]="from A=%d to B=%d\n";
   char s5[15]="file: %s\n";
   char s6[15]="CALC START:\n";
   char s7[40]="p=%d はM素数 t=%4.2f pt=%4.2f\n";
   char s8[50]="p=%d はM素数でない t=%4.2f pt=%4.2f\n";
   char s9[35]="CALC END: t=%4.2f pt=%4.2f\n";
   char s10[5]="%s";
   char s11[19]="\33[0;0HSave File: ";
   char s12[13]="From Anum: ";
   char s13[13]="To   Bnum: ";
   char s14[5]="%d";
   char s15[10]="\33[2J";
   char s16[10]="\n\n\n";
   char s17[50]="\33[8;0H\33[Kp=%d はM素数 t=%4.2f pt=%4.2f\n";
   char s18[60]="\33[8;0H\33[Kp=%d はM素数でない t=%4.2f pt=%4.2f\n";
   FILE *fp;

   pt=0.0;
   st=0.0;

  /* 上の２行はアドレスがころころ変わるのを防ぐ     */
  /* 今の時点で pt=.L3+68 st=.L3+72 ファイル名はs1に */

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#133");
   asm("   bl printf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#92");
   asm("   bl printf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#73");
   asm("   sub r1,fp,#123");
   asm("   bl scanf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#105");
   asm("   bl printf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#123");
   asm("   sub r1,fp,#104");
   asm("   bl scanf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#118");
   asm("   bl printf");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#123");
   asm("   sub r1,fp,#100");
   asm("   bl scanf");

   asm("   bl    gettime");
   asm("   ldr   r2, .L3+68");
   asm("   stmia r2, {r0-r1}");

   asm("   sub r0,fp,#123");
   asm("   sub r1,fp,#128");
   asm("   bl  fopen");

   asm("   mov r10,r0");
   asm("   cmp r10,#0");
   asm("   bne T2");
   asm("   sub r0,fp,#138");
   asm("   bl  puts");
   asm("   mov r0,#1");
   asm("   bl  exit");
   asm("T2:");

   asm("   mov r0,r10");
   asm("   sub r1,fp,#163");
   asm("   ldr r2,[fp,#-104]");
   asm("   ldr r3,[fp,#-100]");
   asm("   bl  fprintf");
   asm("   mov r0,r10");
   asm("   sub r1,fp,#178");
   asm("   sub r2,fp,#123");
   asm("   bl  fprintf");
   asm("   mov r0,r10");
   asm("   sub r1,fp,#193");
   asm("   bl  fprintf");

   asm("   mov r0,r10");
   asm("   bl  fclose");

   asm("   ldr r4,[fp,#-104]");
   asm("   ldr r9,[fp,#-100]");
   asm("   b U3");
   asm("U12:");

   asm("   mov r5,r4");
   asm("   mov r6,#2");

   asm("   b U4");
   asm("U8:   ");
   asm("   mov r3,r6");
   asm("   mul r2,r3,r6");
   asm("   cmp r2,r5");
   asm("   bgt U5");

   asm("   mov r0,r5");
   asm("   mov r1,r6");
   asm("   bl __aeabi_idivmod");
   asm("   cmp r1,#0");
   asm("   bne U6");

   asm("   mov r0,r5");
   asm("   mov r1,r6");
   asm("   bl __aeabi_idiv");
   asm("   mov r5,r0");
   asm("   mov r6,#2");
   asm("   b U4");
   asm("U6:   ");

   asm("   cmp r6,#2");
   asm("   bne U7");
   asm("   add r6,r6,#1");
   asm("   b U4");
   asm("U7:   ");
   asm("   add r6,r6,#2");

   asm("U4:");
   asm("   cmp r6,r5");
   asm("   ble U8");
   asm("U5:   ");

   asm("   cmp r5,r4");
   asm("   bne U9");

   asm("   mov r0,r5");
   asm("   bl prime");
   asm("   cmp r0,#0");
   asm("   beq U10");

   asm("   bl  gettime");
   asm("   mov r7,r0");
   asm("   mov r8,r1");

   asm("   ldr r2,.L3+68");
   asm("   ldmia r2,{r2-r3}");
   asm("   bl  __aeabi_dsub");

   asm("   str r0,[sp,#0]");
   asm("   str r1,[sp,#4]");
   asm("   sub r2,fp,#60");
   asm("   stmia r2,{r0-r1}");

   asm("   mov r0,r7");
   asm("   mov r1,r8");
   asm("   ldr r2,.L3+72");
   asm("   ldmia r2,{r2-r3}");

   asm("   bl __aeabi_dsub");
   asm("   sub r2,fp,#52");
   asm("   stmia r2,{r0-r1}");

   asm("   mov r2,r0");
   asm("   mov r3,r1");

   asm("   sub r0,fp,#233");
   asm("   sub r0,r0,#210");
   asm("   mov r1,r4");
   asm("   bl  printf");

   asm("   sub r0,fp,#123");
   asm("   sub r1,fp,#128");
   asm("   bl  fopen");

   asm("   mov r10,r0");
   asm("   cmp r10,#0");
   asm("   bne T3");
   asm("   sub r0,fp,#138");
   asm("   bl  puts");
   asm("   mov r0,#1");
   asm("   bl  exit");
   asm("T3:");

   asm("   mov r0,r10");
   asm("   sub r1,fp,#233");
   asm("   mov r2,r4");
   asm("   sub r7,fp,#52");
   asm("   ldmia r7,{r7-r8}");
   asm("   stmia sp,{r7-r8}");

   asm("   sub r7,fp,#60");
   asm("   ldmia r7,{r7-r8}");
   asm("   str r7,[sp,#8]");
   asm("   str r8,[sp,#12]");
   asm("   bl  fprintf");

   asm("   mov r0,r10");
   asm("   bl  fclose");

   asm("   b U9");
   asm("U10:   ");

   asm("   bl  gettime");
   asm("   mov r7,r0");
   asm("   mov r8,r1");

   asm("   ldr r2,.L3+68");
   asm("   ldmia r2,{r2-r3}");
   asm("   bl  __aeabi_dsub");

   asm("   stmia sp,{r0-r1}");
   asm("   mov r0,r7");
   asm("   mov r1,r8");
   asm("   ldr r2,.L3+72");
   asm("   ldmia r2,{r2-r3}");

   asm("   bl __aeabi_dsub");
   asm("   mov r2,r0");
   asm("   mov r3,r1");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#253");
   asm("   mov r1,r4");
   asm("   bl  printf");

   asm("U9:");
   asm("   add r4,r4,#1");
   asm("U3:");
   asm("   cmp r4,r9");
   asm("   ble U12");

   asm("   sub r0,fp,#123");
   asm("   sub r1,fp,#128");
   asm("   bl  fopen");

   asm("   mov r10,r0");
   asm("   cmp r10,#0");
   asm("   bne T4");
   asm("   sub r0,fp,#138");
   asm("   bl  puts");
   asm("   mov r0,#1");
   asm("   bl  exit");
   asm("T4:");

   asm("   bl  gettime");
   asm("   mov r5,r0");
   asm("   mov r6,r1");

   asm("   ldr r2,.L3+68");
   asm("   ldmia r2,{r2-r3}");
   asm("   bl  __aeabi_dsub");

   asm("   stmia sp,{r0-r1}");
   asm("   mov r0,r5");
   asm("   mov r1,r6");
   asm("   ldr r2,.L3+72");
   asm("   ldmia r2,{r2-r3}");

   asm("   bl __aeabi_dsub");
   asm("   mov r2,r0");
   asm("   mov r3,r1");

   asm("   mov r0,r10");
   asm("   sub r1,fp,#250");
   asm("   sub r1,r1,#68");
   asm("   bl  fprintf");

   asm("   mov r0,r10");
   asm("   bl  fclose");

   asm("   sub r0,fp,#250");
   asm("   sub r0,r0,#143");
   asm("   bl printf");

}

double gettime(){
   struct timeval tv;
   int a,b,c,d,e;
   double da=1e-6,db,dc,dd;

   asm("   push {r2,r3,r4}");
   asm("   push {r5,r6}");
   asm("   sub r0,fp,#80");
   asm("   mov r1,#0");
   asm("   bl gettimeofday");
   asm("   ldr r0,[fp,#-80]");
   asm("   bl __aeabi_i2d");
   asm("   mov r5,r0");
   asm("   mov r6,r1");
   asm("   ldr r0,[fp,#-76]");
   asm("   bl __aeabi_i2d");
   asm("   sub r2,fp,#52");
   asm("   ldmia r2,{r2-r3}");
   asm("   bl __aeabi_dmul");
   asm("   mov r2,r5");
   asm("   mov r3,r6");
   asm("   bl __aeabi_dadd");
   asm("   pop {r5,r6}");
   asm("   pop {r2,r3,r4}");


}

int prime(int p){
   char a[N+1],x[N];
   int h,i,j,k,s;
   int b=N,c=p,d=-10097,e=-20097,f=-20137,g;
   double da,db,dc,dd;

   char s1[40]="\33[6;0H\33[K%d/%d t=%4.2f pt=%4.2f\n";
   char s2[10]="%s";
   char s3[10]="s=%s\n";
   char s4[10]="%d";
   char s5[10]="f=%f\n";

   da=pt;
   db=st;

  /* 上の２行はアドレスがころころ変わるのを防ぐ */
  /* 今の時点で pt=.L11+56 st=.L11+60            */

   asm("   push {r4,r5,r6,r7,r8,r9,r10}");
   asm("   bl    gettime");
   asm("   ldr   r2, .L11+60");
   asm("   stmia r2, {r0-r1}");

   asm("   ldr r4,[fp,#-72]");
   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   mov r3,#0");
   asm("   mov r7,#0");
   asm("   b T22");

   asm("T23:");
   asm("   strb r3,[r1,r7]");
   asm("   add r7,r7,#1");
   asm("T22:");

   asm("   cmp r7,r4");
   asm("   blt T23");

   asm("   mov r3,#1");
   asm("   strb r3,[r1,#2]");

   asm("   mov r9,#2");
   asm("   b T24");
   asm("T34:");

   asm("   and r0,r9,#0x1f");
   asm("   cmp r0,#0");
   asm("   bne S4");

   asm("   bl  gettime");
   asm("   mov r5,r0");
   asm("   mov r6,r1");

   asm("   ldr r2,.L11+56");
   asm("   ldmia r2,{r2-r3}");
   asm("   bl  __aeabi_dsub");

   asm("   str r0,[sp,#4]");
   asm("   str r1,[sp,#8]");
   asm("   mov r0,r5");
   asm("   mov r1,r6");

   asm("   ldr r2,.L11+60");
   asm("   ldmia r2,{r2-r3}");

   asm("   bl __aeabi_dsub");
   asm("   mov r3,r0");
   asm("   str r1,[sp,#0]");

   asm("   ldr r0,[fp,#-60]");
   asm("   add r0,fp,r0");
   asm("   mov r1,r9");
   asm("   mov r2,r4");
   asm("   bl printf");

   asm("S4:");

   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   ldr r2,[fp,#-64]");
   asm("   add r0,fp,r2");
   asm("   mov r3,#1");
   asm("   mov r7,#0");

   asm("   b T25");
   asm("T26:");

   asm("   ldrb r2,[r1,r7]");
   asm("   strb r2,[r0,r7]");
   asm("   strb r3,[r1,r7]");
   asm("   add r7,r7,#1");

   asm("T25:");

   asm("   cmp r7,r4");
   asm("   blt T26");

   asm("   mov r3,#0");
   asm("   strb r3,[r1,#1]");

   asm("   mov r7,#0");
   asm("   b T27");
   asm("T33:");
   asm("   ldr r2,[fp,#-64]");
   asm("   add r0,fp,r2");
   asm("   ldrb r3,[r0,r7]");

   asm("   cmp r3,#0");
   asm("   beq T28");

   asm("   mov r5,#0");
   asm("   mov r6,r7");
   asm("   mov r8,#0");
   asm("   b T29");
   asm("T30:");
   asm("   ldr r2,[fp,#-64]");
   asm("   add r1,fp,r2");
   asm("   ldrb r0,[r1,r8]");
   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   ldrb r2,[r1,r6]");
   asm("   add r0,r0,r2");

   asm("   mov r2,r5,asr #1");
   asm("   add r5,r0,r2");

   asm("   and r3,r5,#1");

/*
   asm("   str r3,[fp,#-92]");
   asm("   push {r0,r1,r2,r3}");
   asm("   ldr r2,[fp,#-60]");
   asm("   add r0,fp,r2");
   asm("   sub r0,r0,#30");
   asm("   ldr r1,[fp,#-92]");
   asm("   bl printf");
   asm("   pop {r0,r1,r2,r3}");
*/


   asm("   strb r3,[r1,r6]");

   asm("   add r0,r6,#1");
   asm("   mov r1,r4");
   asm("   bl __aeabi_idivmod");
   asm("   mov r6,r1");
   asm("   add r8,r8,#1");
   asm("T29:   ");
   asm("   cmp r8,r4");
   asm("   blt T30");
   asm("   cmp r5,#1");
   asm("   ble T28");

   asm("   b   T31");

   asm("T32:");
   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   mov r0,#0");
   asm("   strb r0,[r1,r6]");

   asm("   add r0,r6,#1");
   asm("   mov r1,r4");
   asm("   bl __aeabi_idivmod");
   asm("   mov r6,r1");

   asm("T31:");
   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   ldrb r0,[r1,r6]");
   asm("   cmp r0,#0");
   asm("   bne T32");
   asm("   mov r0,#1");
   asm("   strb r0,[r1,r6]");
   asm("T28:");
   asm("   add r7,r7,#1");
   asm("T27:");
   asm("   cmp r7,r4");
   asm("   blt T33");
   asm("   add r9,r9,#1");
   asm("T24:");
   asm("   cmp r9,r4");
   asm("   blt T34");

   asm("   ldr r2,[fp,#-68]");
   asm("   add r1,fp,r2");
   asm("   ldrb r0,[r1,#0]");
   asm("   rsb r2,r0,#1");
   asm("   strb r2,[r1,r4]");

   asm("   mov r7,#1");
   asm("   b T35");
   asm("T36:");
   asm("   add r7,r7,#1");
   asm("T35:");
   asm("   ldrb r0,[r1,r7]");
   asm("   ldrb r2,[r1,#0]");
   asm("   cmp r0,r2");
   asm("   beq T36");

   asm("   cmp r7,r4");
   asm("   movne r0,#0");
   asm("   moveq r0,#1");

   asm("   pop {r4,r5,r6,r7,r8,r9,r10}");
}