2013年12月15日日曜日

32bit整数レジスタを自己流実数レジスタとして扱い平方根を求める

akira@dynabookAZ:~/arm1312$ ./c014
Number: 2   2.0を入力
l=200000    レジスタの中でこれを2.0とする
org. sqrt(dd)=1.4142135624  C言語のmath.hで計算させた参考値
sqrt=1.4142135624     
@@@  sqrt1=  @@@   l=  16a09e   +1.4142131805  自己流で計算させた値
Number: 3
l=300000
org. sqrt(dd)=1.7320508076
sqrt=1.7320508076
@@@  sqrt1=  @@@   l=  1bb67a   +1.7320499420
Number: 0.001
l=418
org. sqrt(dd)=0.0316227766
sqrt=0.0316227766
@@@  sqrt1=  @@@   l=    817d   +0.0316133499
Number: 1678
l=68e00000
org. sqrt(dd)=40.9633982965
sqrt=40.9633982965
@@@  sqrt1=  @@@   l= 28f72e9   +40.9655542373
Number: 2.3333
l=255532
org. sqrt(dd)=1.5275143207
sqrt=1.5275143207
@@@  sqrt1=  @@@   l=  1870b2   +1.5275135040
Number: ^C
akira@dynabookAZ:~/arm1312$ cat c014.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define  x8  0x80000000
#define  x7f 0x7fffffff
#define  xpi 0x3243f6

long main1b();
long add0(long l1,long l2);
long sub0(long l1,long l2);
long add1(long l1,long l2);
long sub1(long l1,long l2);
long mul1(long l1,long l2);
long div1(long l1,long l2);
long mod1(long l1,long l2);
void printl(char *txt,long l);
long sin0(long l1);
long cos0(long l1);
long sin1(long l1);
long cos1(long l1);
long atan1(long l1);
long sqrt1(long l1);

main(){
   long l,l1;

   while(1){
      l=main1b();
      l1=sqrt1(l);
      printl("sqrt1=",l1);
   }
}

long sqrt1(long l1){
   long s1,ld,l5,l2;

   s1=l1&x8;
   l1&=x7f;
   if(s1==x8){
      printf("Minus Over !\n");
      exit(0);
   }
   if(l1==0){ return l1;}
   else if(l1==0x100000){ return l1; }
   else if(0<l1&&l1<0x100000){
      l1=mul1(l1,0x2cb9000);
      ld=l2=l1;
      l5=0;
      while(l5<5){
         l5++;
         ld=mul1(l2,l2);
         ld=add1(ld,l1);
         ld=div1(ld,l2)>>1;
         l2=ld;
      }
      l2=div1(l2,0x6b0000);
      return l2;
   }
   else if(0x100000<l1&&l1<0x2cb9000){
      ld=l2=l1;
      l5=0;
      while(l5<5){
         l5++;
         ld=mul1(l2,l2);
         ld=add1(ld,l1);
         ld=div1(ld,l2)>>1;
         l2=ld;
      }
      return l2;
   }
   else{
      l1=div1(l1,0x2cb9000);
      ld=l2=l1;
      l5=0;
      while(l5<5){
         l5++;
         ld=mul1(l2,l2);
         ld=add1(ld,l1);
         ld=div1(ld,l2)>>1;
         l2=ld;
      }
      l2=mul1(l2,0x6b0000);
      return l2;
   }
}


long sin1(long l1){

   l1=add1(l1,0x3f9e0358);
   l1=mod1(l1,0x6487ec);

   if(0<=l1&&l1<=0x1921fb) return sin0(l1);
   else if(0x1921fb<l1&&l1<=0x3243f6){
      l1=sub1(0x3243f6,l1);
      return sin0(l1);
   }
   else if(0x3243f6<l1&&l1<=0x4b65f1){
      l1=sub1(l1,0x3243f6);
      return sin0(l1)^x8;
   }
   else if(0x4b65f1<l1&&l1<=0x6487ec){
      l1=sub1(0x6487ec,l1);
      return sin0(l1)^x8;
   }

}

long cos1(long l1){

   l1=add1(l1,0x3f9e0358);
   l1=mod1(l1,0x6487ec);

  /*printf("171 l1=%lx   ",l1);*/

   if(0<=l1&&l1<=0x1921fb) return cos0(l1);
   else if(0x1921fb<l1&&l1<=0x3243f6){
      l1=sub1(0x3243f6,l1);
      return cos0(l1)^x8;
   }
   else if(0x3243f6<l1&&l1<=0x4b65f1){
      l1=sub1(l1,0x3243f6);
      return cos0(l1)^x8;
   }
   else if(0x4b65f1<l1&&l1<=0x6487ec){
      l1=sub1(0x6487ec,l1);
      return cos0(l1);
   }

}


long sin0(long l1){
   long l,a1,a2,a3,a4,a5;
 
   a1=div1(l1,6<<20);
   a1=mul1(a1,l1);
   a1=mul1(a1,l1);

   a2=div1(a1,20<<20);
   a2=mul1(a2,l1);
   a2=mul1(a2,l1);

   a3=div1(a2,42<<20);
   a3=mul1(a3,l1);
   a3=mul1(a3,l1);

   a4=div1(a3,72<<20);
   a4=mul1(a4,l1);
   a4=mul1(a4,l1);

   a5=div1(a4,110<<20);
   a5=mul1(a5,l1);
   a5=mul1(a5,l1);

   l=sub1(l1,a1);
   l=add1(l,a2);
   l=sub1(l,a3);
   l=add1(l,a4);
   l=sub1(l,a5);
   return l;
}

long cos0(long l1){
   long l,a1,a2,a3,a4,a5,a6;
 
   a1=div1(l1,2<<20);
   a1=mul1(a1,l1);

   a2=div1(a1,12<<20);
   a2=mul1(a2,l1);
   a2=mul1(a2,l1);

   a3=div1(a2,30<<20);
   a3=mul1(a3,l1);
   a3=mul1(a3,l1);

   a4=div1(a3,56<<20);
   a4=mul1(a4,l1);
   a4=mul1(a4,l1);

   a5=div1(a4,90<<20);
   a5=mul1(a5,l1);
   a5=mul1(a5,l1);
/*
   a6=div1(a5,132<<20);
   a6=mul1(a6,l1);
   a6=mul1(a6,l1);
*/
   l=1<<20;
   l=sub1(l,a1);
   l=add1(l,a2);
   l=sub1(l,a3);
   l=add1(l,a4);
   l=sub1(l,a5);
   /*l=add1(l,a6);*/

   return l;
}

long atan1(long l1){
   long l,b,a1,a2,a3,a4,a5;
 
   b=mul1(l1,l1);
   b=mul1(b,l1);
   a1=div1(b,3<<20);

   b=mul1(b,l1);
   b=mul1(b,l1);
   a2=div1(b,5<<20);

   b=mul1(b,l1);
   b=mul1(b,l1);
   a3=div1(b,7<<20);

   b=mul1(b,l1);
   b=mul1(b,l1);
   a4=div1(b,9<<20);

   b=mul1(b,l1);
   b=mul1(b,l1);
   a5=div1(b,11<<20);

   l=sub1(l1,a1);
   l=add1(l,a2);
   l=sub1(l,a3);
   l=add1(l,a4);
   l=sub1(l,a5);
   return l;
}



long add0(long l1,long l2){
   return (l1+l2);
}

long sub0(long l1,long l2){
   if(l1>=l2) return (l1-l2);
   else       return ((l2-l1)^x8);
}

long add1(long l1,long l2){
   long s1,s2;

   s1=l1&x8;
   s2=l2&x8;
   l1&=x7f;
   l2&=x7f;
   if(s1!=x8&&s2!=x8) return add0(l1,l2);
   if(s1!=x8&&s2==x8) return sub0(l1,l2);
   if(s1==x8&&s2!=x8) return sub0(l1,l2)^x8;
   if(s1==x8&&s2==x8) return add0(l1,l2)^x8;
}


long sub1(long l1,long l2){
   long s1,s2;

   s1=l1&x8;
   s2=l2&x8;
   l1&=x7f;
   l2&=x7f;
   if(s1!=x8&&s2!=x8) return sub0(l1,l2);
   if(s1!=x8&&s2==x8) return add0(l1,l2);
   if(s1==x8&&s2!=x8) return add0(l1,l2)^x8;
   if(s1==x8&&s2==x8) return sub0(l1,l2)^x8;
}


long mul1(long l1,long l2){
   long long int ll;
   long s1,s2;

   s1=l1&x8;
   s2=l2&x8;
   l1&=x7f;
   l2&=x7f;
   ll=((long long int)l1*l2)>>20;
   return ((long)ll)^(s1^s2);
}

long div1(long l1,long l2){
   long long int ll;
   long s1,s2;

   s1=l1&x8;
   s2=l2&x8;
   l1&=x7f;
   l2&=x7f;
   ll=(((long long int)l1)<<20)/l2;
   return ((long)ll)^(s1^s2);
}


long mod1(long l1,long l2){
   l1&=x7f;
   l2&=x7f;

   l1%=l2;
   return l1;
}



long main1b(){
   char ss[100]="                            ";
   long l,l2;
   int  i,kd=0,ke=0,js=0,dotf=0,i5;
   double kf=0.0,kg=0.1,l1,dd;
   double ld;

   printf("Number: ");
   scanf("%s",ss);

   i=0;
   while(ss[i]!='\0'){
      if(ss[i]=='-') js=1;
      if(ss[i]=='.') dotf=1;
      if('0'<=ss[i]&&ss[i]<='9'){
         if(dotf==0){
            kd=ke*10+(ss[i]-48);
            ke=kd;
         }
         else{
            kf+=kg*(ss[i]-48);
            kg*=0.1;
         }
      }
      i++;
   }
   if(kd>=2048){
      printf("2048 OVER !\n");
      exit(0);
   }
   dd=kd+kf;
   l=kd<<20;
   l1=kf*16.0;
   l2=(int)l1;
   l+=l2<<16;
   l1=(l1-l2)*16.0;
   l2=(int)l1;
   l+=l2<<12;
   l1=(l1-l2)*16.0;
   l2=(int)l1;
   l+=l2<<8;
   l1=(l1-l2)*16.0;
   l2=(int)l1;
   l+=l2<<4;
   l1=(l1-l2)*16.0;
   l2=(int)l1;
   l+=l2;
   if(js==1) l^=x8;
   if(js==1) dd*=-1.0;
   printf("l=%lx\n",l);

   if(dd<0.0){
      printf("Minus OVER !\n");
      exit(0);
   }
   printf("org. sqrt(dd)=%3.10f\n",sqrt(dd));

   if(dd==0.0) printf("sqrt=0.0\n");
   else if(dd==1.0) printf("sqrt=1.0\n");
   else{
      ld=dd;
      i5=0;
      while(i5<10){
         ld=(ld*ld+dd)/(2*ld);
         i5++;
      }
      printf("sqrt=%3.10f\n",ld);
   }

   return l;
}

void printl(char *txt,long l){
   long i,j;
   int  a,dotf=0;

   printf("@@@  %s  @@@   ",txt);
   printf("l=%8lx   ",l);
   if(l<0){ dotf=1; l^=x8; }
   /*printf("%+5ld.",l);*/
   i=l>>20;
   if(dotf==1) printf("-");
   else        printf("+");
   printf("%ld.",i);
   a=0;
   while(a<10){
      i=l&0x000fffff;
      j=((i<<2)+i)<<1;
      l=j;
      j>>=20;
      printf("%lx",j);
      a++;
   }
   printf("\n");
}


akira@dynabookAZ:~/arm1312$

0 件のコメント:

コメントを投稿