URL
https://opencores.org/ocsvn/eco32/eco32/trunk
Subversion Repositories eco32
[/] [eco32/] [trunk/] [fp/] [implementation/] [arith/] [mmix-arith.c] - Rev 116
Go to most recent revision | Compare with Previous | Blame | View Log
#define sign_bit ((unsigned) 0x80000000) \ #define ROUND_OFF 1 #define ROUND_UP 2 #define ROUND_DOWN 3 #define ROUND_NEAR 4 \ #define X_BIT (1<<8) #define Z_BIT (1<<9) #define U_BIT (1<<10) #define O_BIT (1<<11) #define I_BIT (1<<12) #define W_BIT (1<<13) #define V_BIT (1<<14) #define D_BIT (1<<15) #define E_BIT (1<<18) \ #define zero_exponent (-1000) \ #define bignum_prec 157 \ #define magic_offset 2112 #define origin 37 \ #define buf0 (buf+8) #define buf_max (buf+777) \ /*1:*/ #line 32 "./mmix-arith.w" #include <stdio.h> #include <string.h> #include <ctype.h> /*2:*/ #line 49 "./mmix-arith.w" #ifdef __STDC__ #define ARGS(list) list #else #define ARGS(list) () #endif /*:2*/ #line 36 "./mmix-arith.w" typedef enum{false,true}bool; /*3:*/ #line 60 "./mmix-arith.w" typedef unsigned int tetra; typedef struct{tetra h,l;}octa; /*:3*/ #line 38 "./mmix-arith.w" /*36:*/ #line 605 "./mmix-arith.w" typedef enum{zro,num,inf,nan}ftype; /*:36*//*59:*/ #line 1110 "./mmix-arith.w" typedef struct{ int a; int b; tetra dat[bignum_prec]; }bignum; /*:59*/ #line 39 "./mmix-arith.w" /*4:*/ #line 67 "./mmix-arith.w" octa zero_octa; octa neg_one= {-1,-1}; octa inf_octa= {0x7ff00000,0}; octa standard_NaN= {0x7ff80000,0}; octa aux; bool overflow; /*:4*//*9:*/ #line 174 "./mmix-arith.w" extern octa aux; extern bool overflow; /*:9*//*30:*/ #line 464 "./mmix-arith.w" int cur_round; /*:30*//*32:*/ #line 528 "./mmix-arith.w" int exceptions; /*:32*//*69:*/ #line 1359 "./mmix-arith.w" octa val; char*next_char; /*:69*//*75:*/ #line 1432 "./mmix-arith.w" static char buf[785]= "00000000"; /*:75*/ #line 40 "./mmix-arith.w" /*5:*/ #line 78 "./mmix-arith.w" octa oplus ARGS((octa,octa)); octa oplus(y,z) octa y,z; {octa x; x.h= y.h+z.h; x.l= y.l+z.l; if(x.l<y.l)x.h++; return x; } octa ominus ARGS((octa,octa)); octa ominus(y,z) octa y,z; {octa x; x.h= y.h-z.h; x.l= y.l-z.l; if(x.l> y.l)x.h--; return x; } /*:5*//*6:*/ #line 102 "./mmix-arith.w" octa incr ARGS((octa,int)); octa incr(y,delta) octa y; int delta; {octa x; x.h= y.h;x.l= y.l+delta; if(delta>=0){ if(x.l<y.l)x.h++; }else if(x.l> y.l)x.h--; return x; } /*:6*//*7:*/ #line 117 "./mmix-arith.w" octa shift_left ARGS((octa,int)); octa shift_left(y,s) octa y; int s; { while(s>=32)y.h= y.l,y.l= 0,s-= 32; if(s){register tetra yhl= y.h<<s,ylh= y.l>>(32-s); y.h= yhl+ylh;y.l<<= s; } return y; } octa shift_right ARGS((octa,int,int)); octa shift_right(y,s,u) octa y; int s,u; { while(s>=32)y.l= y.h,y.h= (u?0:-(y.h>>31)),s-= 32; if(s){register tetra yhl= y.h<<(32-s),ylh= y.l>>s; y.h= (u?0:(-(y.h>>31))<<(32-s))+(y.h>>s);y.l= yhl+ylh; } return y; } /*:7*//*8:*/ #line 150 "./mmix-arith.w" octa omult ARGS((octa,octa)); octa omult(y,z) octa y,z; { register int i,j,k; tetra u[4],v[4],w[8]; register tetra t; octa acc; /*10:*/ #line 178 "./mmix-arith.w" u[3]= y.h>>16,u[2]= y.h&0xffff,u[1]= y.l>>16,u[0]= y.l&0xffff; v[3]= z.h>>16,v[2]= z.h&0xffff,v[1]= z.l>>16,v[0]= z.l&0xffff; /*:10*/ #line 159 "./mmix-arith.w" ; for(j= 0;j<4;j++)w[j]= 0; for(j= 0;j<4;j++) if(!v[j])w[j+4]= 0; else{ for(i= k= 0;i<4;i++){ t= u[i]*v[j]+w[i+j]+k; w[i+j]= t&0xffff,k= t>>16; } w[j+4]= k; } /*11:*/ #line 182 "./mmix-arith.w" aux.h= (w[7]<<16)+w[6],aux.l= (w[5]<<16)+w[4]; acc.h= (w[3]<<16)+w[2],acc.l= (w[1]<<16)+w[0]; /*:11*/ #line 170 "./mmix-arith.w" ; return acc; } /*:8*//*12:*/ #line 191 "./mmix-arith.w" octa signed_omult ARGS((octa,octa)); octa signed_omult(y,z) octa y,z; { octa acc; acc= omult(y,z); if(y.h&sign_bit)aux= ominus(aux,z); if(z.h&sign_bit)aux= ominus(aux,y); overflow= (aux.h!=aux.l||(aux.h^(aux.h>>1)^(acc.h&sign_bit))); return acc; } /*:12*//*13:*/ #line 215 "./mmix-arith.w" octa odiv ARGS((octa,octa,octa)); octa odiv(x,y,z) octa x,y,z; { register int i,j,k,n,d; tetra u[8],v[4],q[4],mask,qhat,rhat,vh,vmh; register tetra t; octa acc; /*14:*/ #line 234 "./mmix-arith.w" if(x.h> z.h||(x.h==z.h&&x.l>=z.l)){ aux= y;return x; } /*:14*/ #line 224 "./mmix-arith.w" ; /*15:*/ #line 239 "./mmix-arith.w" u[7]= x.h>>16,u[6]= x.h&0xffff,u[5]= x.l>>16,u[4]= x.l&0xffff; u[3]= y.h>>16,u[2]= y.h&0xffff,u[1]= y.l>>16,u[0]= y.l&0xffff; v[3]= z.h>>16,v[2]= z.h&0xffff,v[1]= z.l>>16,v[0]= z.l&0xffff; /*:15*/ #line 225 "./mmix-arith.w" ; /*16:*/ #line 244 "./mmix-arith.w" for(n= 4;v[n-1]==0;n--); /*:16*/ #line 226 "./mmix-arith.w" ; /*17:*/ #line 250 "./mmix-arith.w" vh= v[n-1]; for(d= 0;vh<0x8000;d++,vh<<= 1); for(j= k= 0;j<n+4;j++){ t= (u[j]<<d)+k; u[j]= t&0xffff,k= t>>16; } for(j= k= 0;j<n;j++){ t= (v[j]<<d)+k; v[j]= t&0xffff,k= t>>16; } vh= v[n-1]; vmh= (n> 1?v[n-2]:0); /*:17*/ #line 227 "./mmix-arith.w" ; for(j= 3;j>=0;j--)/*20:*/ #line 276 "./mmix-arith.w" { /*21:*/ #line 284 "./mmix-arith.w" t= (u[j+n]<<16)+u[j+n-1]; qhat= t/vh,rhat= t-vh*qhat; if(n> 1)while(qhat==0x10000||qhat*vmh> (rhat<<16)+u[j+n-2]){ qhat--,rhat+= vh; if(rhat>=0x10000)break; } /*:21*/ #line 278 "./mmix-arith.w" ; /*22:*/ #line 296 "./mmix-arith.w" for(i= k= 0;i<n;i++){ t= u[i+j]+0xffff0000-k-qhat*v[i]; u[i+j]= t&0xffff,k= 0xffff-(t>>16); } /*:22*/ #line 279 "./mmix-arith.w" ; /*23:*/ #line 305 "./mmix-arith.w" if(u[j+n]!=k){ qhat--; for(i= k= 0;i<n;i++){ t= u[i+j]+v[i]+k; u[i+j]= t&0xffff,k= t>>16; } } /*:23*/ #line 280 "./mmix-arith.w" ; q[j]= qhat; } /*:20*/ #line 228 "./mmix-arith.w" ; /*18:*/ #line 264 "./mmix-arith.w" mask= (1<<d)-1; for(j= 3;j>=n;j--)u[j]= 0; for(k= 0;j>=0;j--){ t= (k<<16)+u[j]; u[j]= t>>d,k= t&mask; } /*:18*/ #line 229 "./mmix-arith.w" ; /*19:*/ #line 272 "./mmix-arith.w" acc.h= (q[3]<<16)+q[2],acc.l= (q[1]<<16)+q[0]; aux.h= (u[3]<<16)+u[2],aux.l= (u[1]<<16)+u[0]; /*:19*/ #line 230 "./mmix-arith.w" ; return acc; } /*:13*//*24:*/ #line 317 "./mmix-arith.w" octa signed_odiv ARGS((octa,octa)); octa signed_odiv(y,z) octa y,z; { octa yy,zz,q; register int sy,sz; if(y.h&sign_bit)sy= 2,yy= ominus(zero_octa,y); else sy= 0,yy= y; if(z.h&sign_bit)sz= 1,zz= ominus(zero_octa,z); else sz= 0,zz= z; q= odiv(zero_octa,yy,zz); overflow= false; switch(sy+sz){ case 2+1:aux= ominus(zero_octa,aux); if(q.h==sign_bit)overflow= true; case 0+0:return q; case 2+0:if(aux.h||aux.l)aux= ominus(zz,aux); goto negate_q; case 0+1:if(aux.h||aux.l)aux= ominus(aux,zz); negate_q:if(aux.h||aux.l)return ominus(neg_one,q); else return ominus(zero_octa,q); } } /*:24*//*25:*/ #line 346 "./mmix-arith.w" octa oand ARGS((octa,octa)); octa oand(y,z) octa y,z; {octa x; x.h= y.h&z.h;x.l= y.l&z.l; return x; } octa oandn ARGS((octa,octa)); octa oandn(y,z) octa y,z; {octa x; x.h= y.h&~z.h;x.l= y.l&~z.l; return x; } octa oxor ARGS((octa,octa)); octa oxor(y,z) octa y,z; {octa x; x.h= y.h^z.h;x.l= y.l^z.l; return x; } /*:25*//*26:*/ #line 387 "./mmix-arith.w" int count_bits ARGS((tetra)); int count_bits(x) tetra x; { register int xx= x; xx= xx-((xx>>1)&0x55555555); xx= (xx&0x33333333)+((xx>>2)&0x33333333); xx= (xx+(xx>>4))&0x0f0f0f0f; xx= xx+(xx>>8); return(xx+(xx>>16))&0xff; } /*:26*//*27:*/ #line 403 "./mmix-arith.w" tetra byte_diff ARGS((tetra,tetra)); tetra byte_diff(y,z) tetra y,z; { register tetra d= (y&0x00ff00ff)+0x01000100-(z&0x00ff00ff); register tetra m= d&0x01000100; register tetra x= d&(m-(m>>8)); d= ((y>>8)&0x00ff00ff)+0x01000100-((z>>8)&0x00ff00ff); m= d&0x01000100; return x+((d&(m-(m>>8)))<<8); } /*:27*//*28:*/ #line 421 "./mmix-arith.w" tetra wyde_diff ARGS((tetra,tetra)); tetra wyde_diff(y,z) tetra y,z; { register tetra a= ((y>>16)-(z>>16))&0x10000; register tetra b= ((y&0xffff)-(z&0xffff))&0x10000; return y-(z^((y^z)&(b-a-(b>>16)))); } /*:28*//*29:*/ #line 434 "./mmix-arith.w" octa bool_mult ARGS((octa,octa,bool)); octa bool_mult(y,z,xor) octa y,z; bool xor; { octa o,x; register tetra a,b,c; register int k; for(k= 0,o= y,x= zero_octa;o.h||o.l;k++,o= shift_right(o,8,1)) if(o.l&0xff){ a= ((z.h>>k)&0x01010101)*0xff; b= ((z.l>>k)&0x01010101)*0xff; c= (o.l&0xff)*0x01010101; if(xor)x.h^= a&c,x.l^= b&c; else x.h|= a&c,x.l|= b&c; } return x; } /*:29*//*31:*/ #line 503 "./mmix-arith.w" octa fpack ARGS((octa,int,char,int)); octa fpack(f,e,s,r) octa f; int e; char s; int r; { octa o; if(e> 0x7fd)e= 0x7ff,o= zero_octa; else{ if(e<0){ if(e<-54)o.h= 0,o.l= 1; else{octa oo; o= shift_right(f,-e,1); oo= shift_left(o,-e); if(oo.l!=f.l||oo.h!=f.h)o.l|= 1; } e= 0; }else o= f; } /*33:*/ #line 533 "./mmix-arith.w" if(o.l&3)exceptions|= X_BIT; switch(r){ case ROUND_DOWN:if(s=='-')o= incr(o,3);break; case ROUND_UP:if(s!='-')o= incr(o,3); case ROUND_OFF:break; case ROUND_NEAR:o= incr(o,o.l&4?2:1);break; } o= shift_right(o,2,1); o.h+= e<<20; if(o.h>=0x7ff00000)exceptions|= O_BIT+X_BIT; else if(o.h<0x100000)exceptions|= U_BIT; if(s=='-')o.h|= sign_bit; return o; /*:33*/ #line 525 "./mmix-arith.w" ; } /*:31*//*34:*/ #line 551 "./mmix-arith.w" tetra sfpack ARGS((octa,int,char,int)); tetra sfpack(f,e,s,r) octa f; int e; char s; int r; { register tetra o; if(e> 0x47d)e= 0x47f,o= 0; else{ o= shift_left(f,3).h; if(f.l&0x1fffffff)o|= 1; if(e<0x380){ if(e<0x380-25)o= 1; else{register tetra o0,oo; o0= o; o= o>>(0x380-e); oo= o<<(0x380-e); if(oo!=o0)o|= 1; } e= 0x380; } } /*35:*/ #line 579 "./mmix-arith.w" if(o&3)exceptions|= X_BIT; switch(r){ case ROUND_DOWN:if(s=='-')o+= 3;break; case ROUND_UP:if(s!='-')o+= 3; case ROUND_OFF:break; case ROUND_NEAR:o+= (o&4?2:1);break; } o= o>>2; o+= (e-0x380)<<23; if(o>=0x7f800000)exceptions|= O_BIT+X_BIT; else if(o<0x100000)exceptions|= U_BIT; if(s=='-')o|= sign_bit; return o; /*:35*/ #line 576 "./mmix-arith.w" ; } /*:34*//*37:*/ #line 608 "./mmix-arith.w" ftype funpack ARGS((octa,octa*,int*,char*)); ftype funpack(x,f,e,s) octa x; octa*f; int*e; char*s; { register int ee; exceptions= 0; *s= (x.h&sign_bit?'-':'+'); *f= shift_left(x,2); f->h&= 0x3fffff; ee= (x.h>>20)&0x7ff; if(ee){ *e= ee-1; f->h|= 0x400000; return(ee<0x7ff?num:f->h==0x400000&&!f->l?inf:nan); } if(!x.l&&!f->h){ *e= zero_exponent;return zro; } do{ee--;*f= shift_left(*f,1);}while(!(f->h&0x400000)); *e= ee;return num; } /*:37*//*38:*/ #line 634 "./mmix-arith.w" ftype sfunpack ARGS((tetra,octa*,int*,char*)); ftype sfunpack(x,f,e,s) tetra x; octa*f; int*e; char*s; { register int ee; exceptions= 0; *s= (x&sign_bit?'-':'+'); f->h= (x>>1)&0x3fffff,f->l= x<<31; ee= (x>>23)&0xff; if(ee){ *e= ee+0x380-1; f->h|= 0x400000; return(ee<0xff?num:(x&0x7fffffff)==0x7f800000?inf:nan); } if(!(x&0x7fffffff)){ *e= zero_exponent;return zro; } do{ee--;*f= shift_left(*f,1);}while(!(f->h&0x400000)); *e= ee+0x380;return num; } /*:38*//*39:*/ #line 663 "./mmix-arith.w" octa load_sf ARGS((tetra)); octa load_sf(z) tetra z; { octa f,x;int e;char s;ftype t; t= sfunpack(z,&f,&e,&s); switch(t){ case zro:x= zero_octa;break; case num:return fpack(f,e,s,ROUND_OFF); case inf:x= inf_octa;break; case nan:x= shift_right(f,2,1);x.h|= 0x7ff00000;break; } if(s=='-')x.h|= sign_bit; return x; } /*:39*//*40:*/ #line 680 "./mmix-arith.w" tetra store_sf ARGS((octa)); tetra store_sf(x) octa x; { octa f;tetra z;int e;char s;ftype t; t= funpack(x,&f,&e,&s); switch(t){ case zro:z= 0;break; case num:return sfpack(f,e,s,cur_round); case inf:z= 0x7f800000;break; case nan:if(!(f.h&0x200000)){ f.h|= 0x200000;exceptions|= I_BIT; } z= 0x7f800000|(f.h<<1)|(f.l>>31);break; } if(s=='-')z|= sign_bit; return z; } /*:40*//*41:*/ #line 705 "./mmix-arith.w" octa fmult ARGS((octa,octa)); octa fmult(y,z) octa y,z; { ftype yt,zt; int ye,ze; char ys,zs; octa x,xf,yf,zf; register int xe; register char xs; yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); xs= ys+zs-'+'; switch(4*yt+zt){ /*42:*/ #line 731 "./mmix-arith.w" case 4*nan+nan:if(!(y.h&0x80000))exceptions|= I_BIT; case 4*zro+nan:case 4*num+nan:case 4*inf+nan: if(!(z.h&0x80000))exceptions|= I_BIT,z.h|= 0x80000; return z; case 4*nan+zro:case 4*nan+num:case 4*nan+inf: if(!(y.h&0x80000))exceptions|= I_BIT,y.h|= 0x80000; return y; /*:42*/ #line 720 "./mmix-arith.w" ; case 4*zro+zro:case 4*zro+num:case 4*num+zro:x= zero_octa;break; case 4*num+inf:case 4*inf+num:case 4*inf+inf:x= inf_octa;break; case 4*zro+inf:case 4*inf+zro:x= standard_NaN; exceptions|= I_BIT;break; case 4*num+num:/*43:*/ #line 740 "./mmix-arith.w" xe= ye+ze-0x3fd; x= omult(yf,shift_left(zf,9)); if(aux.h>=0x400000)xf= aux; else xf= shift_left(aux,1),xe--; if(x.h||x.l)xf.l|= 1; return fpack(xf,xe,xs,cur_round); /*:43*/ #line 725 "./mmix-arith.w" ; } if(xs=='-')x.h|= sign_bit; return x; } /*:41*//*44:*/ #line 748 "./mmix-arith.w" octa fdivide ARGS((octa,octa)); octa fdivide(y,z) octa y,z; { ftype yt,zt; int ye,ze; char ys,zs; octa x,xf,yf,zf; register int xe; register char xs; yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); xs= ys+zs-'+'; switch(4*yt+zt){ /*42:*/ #line 731 "./mmix-arith.w" case 4*nan+nan:if(!(y.h&0x80000))exceptions|= I_BIT; case 4*zro+nan:case 4*num+nan:case 4*inf+nan: if(!(z.h&0x80000))exceptions|= I_BIT,z.h|= 0x80000; return z; case 4*nan+zro:case 4*nan+num:case 4*nan+inf: if(!(y.h&0x80000))exceptions|= I_BIT,y.h|= 0x80000; return y; /*:42*/ #line 763 "./mmix-arith.w" ; case 4*zro+inf:case 4*zro+num:case 4*num+inf:x= zero_octa;break; case 4*num+zro:exceptions|= Z_BIT; case 4*inf+num:case 4*inf+zro:x= inf_octa;break; case 4*zro+zro:case 4*inf+inf:x= standard_NaN; exceptions|= I_BIT;break; case 4*num+num:/*45:*/ #line 775 "./mmix-arith.w" xe= ye-ze+0x3fd; xf= odiv(yf,zero_octa,shift_left(zf,9)); if(xf.h>=0x800000){ aux.l|= xf.l&1; xf= shift_right(xf,1,1); xe++; } if(aux.h||aux.l)xf.l|= 1; return fpack(xf,xe,xs,cur_round); /*:45*/ #line 769 "./mmix-arith.w" ; } if(xs=='-')x.h|= sign_bit; return x; } /*:44*//*46:*/ #line 790 "./mmix-arith.w" octa fplus ARGS((octa,octa)); octa fplus(y,z) octa y,z; { ftype yt,zt; int ye,ze; char ys,zs; octa x,xf,yf,zf; register int xe,d; register char xs; yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); switch(4*yt+zt){ /*42:*/ #line 731 "./mmix-arith.w" case 4*nan+nan:if(!(y.h&0x80000))exceptions|= I_BIT; case 4*zro+nan:case 4*num+nan:case 4*inf+nan: if(!(z.h&0x80000))exceptions|= I_BIT,z.h|= 0x80000; return z; case 4*nan+zro:case 4*nan+num:case 4*nan+inf: if(!(y.h&0x80000))exceptions|= I_BIT,y.h|= 0x80000; return y; /*:42*/ #line 804 "./mmix-arith.w" ; case 4*zro+num:return fpack(zf,ze,zs,ROUND_OFF);break; case 4*num+zro:return fpack(yf,ye,ys,ROUND_OFF);break; case 4*inf+inf:if(ys!=zs){ exceptions|= I_BIT;x= standard_NaN;xs= zs;break; } case 4*num+inf:case 4*zro+inf:x= inf_octa;xs= zs;break; case 4*inf+num:case 4*inf+zro:x= inf_octa;xs= ys;break; case 4*num+num:if(y.h!=(z.h^0x80000000)||y.l!=z.l) /*47:*/ #line 821 "./mmix-arith.w" {octa o,oo; if(ye<ze||(ye==ze&&(yf.h<zf.h||(yf.h==zf.h&&yf.l<zf.l)))) /*48:*/ #line 839 "./mmix-arith.w" { o= yf,yf= zf,zf= o; d= ye,ye= ze,ze= d; d= ys,ys= zs,zs= d; } /*:48*/ #line 824 "./mmix-arith.w" ; d= ye-ze; xs= ys,xe= ye; if(d)/*49:*/ #line 859 "./mmix-arith.w" { if(d<=2)zf= shift_right(zf,d,1); else if(d> 53)zf.h= 0,zf.l= 1; else{ if(ys!=zs)d--,xe--,yf= shift_left(yf,1); o= zf; zf= shift_right(o,d,1); oo= shift_left(zf,d); if(oo.l!=o.l||oo.h!=o.h)zf.l|= 1; } } /*:49*/ #line 827 "./mmix-arith.w" ; if(ys==zs){ xf= oplus(yf,zf); if(xf.h>=0x800000)xe++,d= xf.l&1,xf= shift_right(xf,1,1),xf.l|= d; }else{ xf= ominus(yf,zf); if(xf.h>=0x800000)xe++,d= xf.l&1,xf= shift_right(xf,1,1),xf.l|= d; else while(xf.h<0x400000)xe--,xf= shift_left(xf,1); } return fpack(xf,xe,xs,cur_round); } /*:47*/ #line 813 "./mmix-arith.w" ; case 4*zro+zro:x= zero_octa; xs= (ys==zs?ys:cur_round==ROUND_DOWN?'-':'+');break; } if(xs=='-')x.h|= sign_bit; return x; } /*:46*//*50:*/ #line 883 "./mmix-arith.w" int fepscomp ARGS((octa,octa,octa,int)); int fepscomp(y,z,e,s) octa y,z,e; int s; { octa yf,zf,ef,o,oo; int ye,ze,ee; char ys,zs,es; register int yt,zt,et,d; et= funpack(e,&ef,&ee,&es); if(es=='-')return 2; switch(et){ case nan:return 2; case inf:ee= 10000; case num:case zro:break; } yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); switch(4*yt+zt){ case 4*nan+nan:case 4*nan+inf:case 4*nan+num:case 4*nan+zro: case 4*inf+nan:case 4*num+nan:case 4*zro+nan:return 2; case 4*inf+inf:return(ys==zs||ee>=1023); case 4*inf+num:case 4*inf+zro:case 4*num+inf:case 4*zro+inf: return(s&&ee>=1022); case 4*zro+zro:return 1; case 4*zro+num:case 4*num+zro:if(!s)return 0; case 4*num+num:break; } /*51:*/ #line 919 "./mmix-arith.w" /*52:*/ #line 934 "./mmix-arith.w" if(ye<0&&yt!=zro)yf= shift_left(y,2),ye= 0; if(ze<0&&zt!=zro)zf= shift_left(z,2),ze= 0; /*:52*/ #line 920 "./mmix-arith.w" ; if(ye<ze||(ye==ze&&(yf.h<zf.h||(yf.h==zf.h&&yf.l<zf.l)))) /*48:*/ #line 839 "./mmix-arith.w" { o= yf,yf= zf,zf= o; d= ye,ye= ze,ze= d; d= ys,ys= zs,zs= d; } /*:48*/ #line 922 "./mmix-arith.w" ; if(ze==zero_exponent)ze= ye; d= ye-ze; if(!s)ee-= d; if(ee>=1023)return 1; /*53:*/ #line 956 "./mmix-arith.w" if(d> 54)o= zero_octa,oo= zf; else o= shift_right(zf,d,1),oo= shift_left(o,d); if(oo.h!=zf.h||oo.l!=zf.l){ if(ee<1020)return 0; o= incr(o,ys==zs?0:1); } o= (ys==zs?ominus(yf,o):oplus(yf,o)); /*:53*/ #line 927 "./mmix-arith.w" ; if(!o.h&&!o.l)return 1; if(ee<968)return 0; if(ee>=1021)ef= shift_left(ef,ee-1021); else ef= shift_right(ef,1021-ee,1); return o.h<ef.h||(o.h==ef.h&&o.l<=ef.l); /*:51*/ #line 912 "./mmix-arith.w" ; } /*:50*//*54:*/ #line 972 "./mmix-arith.w" static void bignum_times_ten ARGS((bignum*)); static void bignum_dec ARGS((bignum*,bignum*,tetra)); static int bignum_compare ARGS((bignum*,bignum*)); void print_float ARGS((octa)); void print_float(x) octa x; { /*56:*/ #line 1035 "./mmix-arith.w" octa f,g; register int e; register int j,k; /*:56*//*66:*/ #line 1281 "./mmix-arith.w" bignum ff,gg; bignum tt; char s[18]; register char*p; /*:66*/ #line 980 "./mmix-arith.w" ; if(x.h&sign_bit)printf("-"); /*55:*/ #line 1019 "./mmix-arith.w" f= shift_left(x,1); e= f.h>>21; f.h&= 0x1fffff; if(!f.h&&!f.l)/*57:*/ #line 1045 "./mmix-arith.w" { if(!e){ printf("0.");return; } if(e==0x7ff){ printf("Inf");return; } e--; f.h= 0x3fffff,f.l= 0xffffffff; g.h= 0x400000,g.l= 2; } /*:57*/ #line 1023 "./mmix-arith.w" else{ g= incr(f,1); f= incr(f,-1); if(!e)e= 1; else if(e==0x7ff){ printf("NaN"); if(g.h==0x100000&&g.l==1)return; e= 0x3ff; }else f.h|= 0x200000,g.h|= 0x200000; } /*:55*/ #line 983 "./mmix-arith.w" ; /*63:*/ #line 1195 "./mmix-arith.w" k= (magic_offset-e)/28; ff.dat[k-1]= shift_right(f,magic_offset+28-e-28*k,1).l&0xfffffff; gg.dat[k-1]= shift_right(g,magic_offset+28-e-28*k,1).l&0xfffffff; ff.dat[k]= shift_right(f,magic_offset-e-28*k,1).l&0xfffffff; gg.dat[k]= shift_right(g,magic_offset-e-28*k,1).l&0xfffffff; ff.dat[k+1]= shift_left(f,e+28*k-(magic_offset-28)).l&0xfffffff; gg.dat[k+1]= shift_left(g,e+28*k-(magic_offset-28)).l&0xfffffff; ff.a= (ff.dat[k-1]?k-1:k); ff.b= (ff.dat[k+1]?k+1:k); gg.a= (gg.dat[k-1]?k-1:k); gg.b= (gg.dat[k+1]?k+1:k); /*:63*/ #line 984 "./mmix-arith.w" ; /*64:*/ #line 1223 "./mmix-arith.w" if(e> 0x401)/*65:*/ #line 1254 "./mmix-arith.w" {register int open= x.l&1; tt.dat[origin]= 10; tt.a= tt.b= origin; for(e= 1;bignum_compare(&gg,&tt)>=open;e++) bignum_times_ten(&tt); p= s; while(1){ bignum_times_ten(&ff); bignum_times_ten(&gg); for(j= '0';bignum_compare(&ff,&tt)>=0;j++) bignum_dec(&ff,&tt,0x10000000),bignum_dec(&gg,&tt,0x10000000); if(bignum_compare(&gg,&tt)>=open)break; *p++= j; if(ff.a==bignum_prec-1&&!open) goto done; } for(k= j;bignum_compare(&gg,&tt)>=open;k++)bignum_dec(&gg,&tt,0x10000000); *p++= (j+1+k)>>1; done:; } /*:65*/ #line 1224 "./mmix-arith.w" else{ if(ff.a> origin)ff.dat[origin]= 0; for(e= 1,p= s;gg.a> origin||ff.dat[origin]==gg.dat[origin];){ if(gg.a> origin)e--; else*p++= ff.dat[origin]+'0',ff.dat[origin]= 0,gg.dat[origin]= 0; bignum_times_ten(&ff); bignum_times_ten(&gg); } *p++= ((ff.dat[origin]+1+gg.dat[origin])>>1)+'0'; } *p= '\0'; /*:64*/ #line 985 "./mmix-arith.w" ; /*67:*/ #line 1296 "./mmix-arith.w" if(e> 17||e<(int)strlen(s)-17) printf("%c%s%se%d",s[0],(s[1]?".":""),s+1,e-1); else if(e<0)printf(".%0*d%s",-e,0,s); else if(strlen(s)>=e)printf("%.*s.%s",e,s,s+e); else printf("%s%0*d.",s,e-(int)strlen(s),0); /*:67*/ #line 986 "./mmix-arith.w" ; } /*:54*//*60:*/ #line 1120 "./mmix-arith.w" static void bignum_times_ten(f) bignum*f; { register tetra*p,*q;register tetra x,carry; for(p= &f->dat[f->b],q= &f->dat[f->a],carry= 0;p>=q;p--){ x= *p*10+carry; *p= x&0xfffffff; carry= x>>28; } *p= carry; if(carry)f->a--; if(f->dat[f->b]==0&&f->b> f->a)f->b--; } /*:60*//*61:*/ #line 1138 "./mmix-arith.w" static int bignum_compare(f,g) bignum*f,*g; { register tetra*p,*pp,*q,*qq; if(f->a!=g->a)return f->a> g->a?-1:1; pp= &f->dat[f->b],qq= &g->dat[g->b]; for(p= &f->dat[f->a],q= &g->dat[g->a];p<=pp;p++,q++){ if(*p!=*q)return*p<*q?-1:1; if(q==qq)return p<pp; } return-1; } /*:61*//*62:*/ #line 1155 "./mmix-arith.w" static void bignum_dec(f,g,r) bignum*f,*g; tetra r; { register tetra*p,*q,*qq; register int x,borrow; while(g->b> f->b)f->dat[++f->b]= 0; qq= &g->dat[g->a]; for(p= &f->dat[g->b],q= &g->dat[g->b],borrow= 0;q>=qq;p--,q--){ x= *p-*q-borrow; if(x>=0)borrow= 0,*p= x; else borrow= 1,*p= x+r; } for(;borrow;p--) if(*p)borrow= 0,*p= *p-1; else*p= r-1; while(f->dat[f->a]==0){ if(f->a==f->b){ f->a= f->b= bignum_prec-1,f->dat[bignum_prec-1]= 0; return; } f->a++; } while(f->dat[f->b]==0)f->b--; } /*:62*//*68:*/ #line 1340 "./mmix-arith.w" static void bignum_double ARGS((bignum*)); int scan_const ARGS((char*)); int scan_const(s) char*s; { /*70:*/ #line 1363 "./mmix-arith.w" register char*p,*q; register bool NaN; int sign; /*:70*//*76:*/ #line 1435 "./mmix-arith.w" register char*dec_pt; register int exp; register int zeros; /*:76*//*81:*/ #line 1503 "./mmix-arith.w" register int k,x; register char*pp; bignum ff,tt; /*:81*/ #line 1346 "./mmix-arith.w" ; val.h= val.l= 0; p= s; if(*p=='+'||*p=='-')sign= *p++;else sign= '+'; if(strncmp(p,"NaN",3)==0)NaN= true,p+= 3; else NaN= false; if((isdigit(*p)&&!NaN)||(*p=='.'&&isdigit(*(p+1)))) /*73:*/ #line 1396 "./mmix-arith.w" { for(q= buf0,dec_pt= (char*)0;isdigit(*p);p++){ val= oplus(val,shift_left(val,2)); val= incr(shift_left(val,1),*p-'0'); if(q> buf0||*p!='0') if(q<buf_max)*q++= *p; else if(*(q-1)=='0')*(q-1)= *p; } if(NaN)*q++= '1'; if(*p=='.')/*74:*/ #line 1415 "./mmix-arith.w" { dec_pt= q; p++; for(zeros= 0;isdigit(*p);p++) if(*p=='0'&&q==buf0)zeros++; else if(q<buf_max)*q++= *p; else if(*(q-1)=='0')*(q-1)= *p; } /*:74*/ #line 1406 "./mmix-arith.w" ; next_char= p; if(*p=='e'&&!NaN)/*77:*/ #line 1447 "./mmix-arith.w" {register char exp_sign; p++; if(*p=='+'||*p=='-')exp_sign= *p++;else exp_sign= '+'; if(isdigit(*p)){ for(exp= *p++-'0';isdigit(*p);p++) if(exp<1000)exp= 10*exp+*p-'0'; if(!dec_pt)dec_pt= q,zeros= 0; if(exp_sign=='-')exp= -exp; next_char= p; } } /*:77*/ #line 1408 "./mmix-arith.w" else exp= 0; if(dec_pt)/*78:*/ #line 1460 "./mmix-arith.w" { /*79:*/ #line 1477 "./mmix-arith.w" x= buf+341+zeros-dec_pt-exp; if(q==buf0||x>=1413){ make_it_zero:exp= -99999;goto packit; } if(x<0){ make_it_infinite:exp= 99999;goto packit; } ff.a= x/9; for(p= q;p<q+8;p++)*p= '0'; q= q-1-(q+341+zeros-dec_pt-exp)%9; for(p= buf0-x%9,k= ff.a;p<=q&&k<=156;p+= 9,k++) /*80:*/ #line 1497 "./mmix-arith.w" { for(x= *p-'0',pp= p+1;pp<p+9;pp++)x= 10*x+*pp-'0'; ff.dat[k]= x; } /*:80*/ #line 1490 "./mmix-arith.w" ; ff.b= k-1; for(x= 0;p<=q;p+= 9)if(strncmp(p,"000000000",9)!=0)x= 1; ff.dat[156]+= x; while(ff.dat[ff.b]==0)ff.b--; /*:79*/ #line 1462 "./mmix-arith.w" ; /*83:*/ #line 1526 "./mmix-arith.w" val= zero_octa; if(ff.a> 36){ for(exp= 0x3fe;ff.a> 36;exp--)bignum_double(&ff); for(k= 54;k;k--){ if(ff.dat[36]){ if(k>=32)val.h|= 1<<(k-32);else val.l|= 1<<k; ff.dat[36]= 0; if(ff.b==36)break; } bignum_double(&ff); } }else{ tt.a= tt.b= 36,tt.dat[36]= 2; for(exp= 0x3fe;bignum_compare(&ff,&tt)>=0;exp++)bignum_double(&tt); for(k= 54;k;k--){ bignum_double(&ff); if(bignum_compare(&ff,&tt)>=0){ if(k>=32)val.h|= 1<<(k-32);else val.l|= 1<<k; bignum_dec(&ff,&tt,1000000000); if(ff.a==bignum_prec-1)break; } } } if(k==0)val.l|= 1; /*:83*/ #line 1463 "./mmix-arith.w" ; packit:/*84:*/ #line 1559 "./mmix-arith.w" val= fpack(val,exp,sign,ROUND_NEAR); if(NaN){ if((val.h&0x7fffffff)==0x40000000)val.h|= 0x7fffffff,val.l= 0xffffffff; else if((val.h&0x7fffffff)==0x3ff00000&&!val.l)val.h|= 0x40000000,val.l= 1; else val.h|= 0x40000000; } /*:84*/ #line 1464 "./mmix-arith.w" ; return 1; } /*:78*/ #line 1410 "./mmix-arith.w" ; if(sign=='-')val= ominus(zero_octa,val); return 0; } /*:73*/ #line 1353 "./mmix-arith.w" ; if(NaN)/*71:*/ #line 1368 "./mmix-arith.w" { next_char= p; val.h= 0x600000,exp= 0x3fe; goto packit; } /*:71*/ #line 1354 "./mmix-arith.w" ; if(strncmp(p,"Inf",3)==0)/*72:*/ #line 1375 "./mmix-arith.w" { next_char= p+3; goto make_it_infinite; } /*:72*/ #line 1355 "./mmix-arith.w" ; no_const_found:next_char= s;return-1; } /*:68*//*82:*/ #line 1511 "./mmix-arith.w" static void bignum_double(f) bignum*f; { register tetra*p,*q;register int x,carry; for(p= &f->dat[f->b],q= &f->dat[f->a],carry= 0;p>=q;p--){ x= *p+*p+carry; if(x>=1000000000)carry= 1,*p= x-1000000000; else carry= 0,*p= x; } *p= carry; if(carry)f->a--; if(f->dat[f->b]==0&&f->b> f->a)f->b--; } /*:82*//*85:*/ #line 1575 "./mmix-arith.w" int fcomp ARGS((octa,octa)); int fcomp(y,z) octa y,z; { ftype yt,zt; int ye,ze; char ys,zs; octa yf,zf; register int x; yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); switch(4*yt+zt){ case 4*nan+nan:case 4*zro+nan:case 4*num+nan:case 4*inf+nan: case 4*nan+zro:case 4*nan+num:case 4*nan+inf:return 2; case 4*zro+zro:return 0; case 4*zro+num:case 4*num+zro:case 4*zro+inf:case 4*inf+zro: case 4*num+num:case 4*num+inf:case 4*inf+num:case 4*inf+inf: if(ys!=zs)x= 1; else if(y.h> z.h)x= 1; else if(y.h<z.h)x= -1; else if(y.l> z.l)x= 1; else if(y.l<z.l)x= -1; else return 0; break; } return(ys=='-'?-x:x); } /*:85*//*86:*/ #line 1608 "./mmix-arith.w" octa fintegerize ARGS((octa,int)); octa fintegerize(z,r) octa z; int r; { ftype zt; int ze; char zs; octa xf,zf; zt= funpack(z,&zf,&ze,&zs); if(!r)r= cur_round; switch(zt){ case nan:if(!(z.h&0x80000)){exceptions|= I_BIT;z.h|= 0x80000;} case inf:case zro:return z; case num:/*87:*/ #line 1627 "./mmix-arith.w" if(ze>=1074)return fpack(zf,ze,zs,ROUND_OFF); if(ze<=1020)xf.h= 0,xf.l= 1; else{octa oo; xf= shift_right(zf,1074-ze,1); oo= shift_left(xf,1074-ze); if(oo.l!=zf.l||oo.h!=zf.h)xf.l|= 1; } switch(r){ case ROUND_DOWN:if(zs=='-')xf= incr(xf,3);break; case ROUND_UP:if(zs!='-')xf= incr(xf,3); case ROUND_OFF:break; case ROUND_NEAR:xf= incr(xf,xf.l&4?2:1);break; } xf.l&= 0xfffffffc; if(ze>=1022)return fpack(shift_left(xf,1074-ze),ze,zs,ROUND_OFF); if(xf.l)xf.h= 0x3ff00000,xf.l= 0; if(zs=='-')xf.h|= sign_bit; return xf; /*:87*/ #line 1623 "./mmix-arith.w" ; } } /*:86*//*88:*/ #line 1650 "./mmix-arith.w" octa fixit ARGS((octa,int)); octa fixit(z,r) octa z; int r; { ftype zt; int ze; char zs; octa zf,o; zt= funpack(z,&zf,&ze,&zs); if(!r)r= cur_round; switch(zt){ case nan:case inf:exceptions|= I_BIT;return z; case zro:return zero_octa; case num:if(funpack(fintegerize(z,r),&zf,&ze,&zs)==zro)return zero_octa; if(ze<=1076)o= shift_right(zf,1076-ze,1); else{ if(ze> 1085||(ze==1085&&(zf.h> 0x400000|| (zf.h==0x400000&&(zf.l||zs!='-')))))exceptions|= W_BIT; if(ze>=1140)return zero_octa; o= shift_left(zf,ze-1076); } return(zs=='-'?ominus(zero_octa,o):o); } } /*:88*//*89:*/ #line 1681 "./mmix-arith.w" octa floatit ARGS((octa,int,int,int)); octa floatit(z,r,u,p) octa z; int r; int u; int p; { int e;char s; register int t; exceptions= 0; if(!z.h&&!z.l)return zero_octa; if(!r)r= cur_round; if(!u&&(z.h&sign_bit))s= '-',z= ominus(zero_octa,z);else s= '+'; e= 1076; while(z.h<0x400000)e--,z= shift_left(z,1); while(z.h>=0x800000){ e++; t= z.l&1; z= shift_right(z,1,1); z.l|= t; } if(p)/*90:*/ #line 1707 "./mmix-arith.w" { register int ex;register tetra t; t= sfpack(z,e,s,r); ex= exceptions; sfunpack(t,&z,&e,&s); exceptions= ex; } /*:90*/ #line 1703 "./mmix-arith.w" ; return fpack(z,e,s,r); } /*:89*//*91:*/ #line 1718 "./mmix-arith.w" octa froot ARGS((octa,int)); octa froot(z,r) octa z; int r; { ftype zt; int ze; char zs; octa x,xf,rf,zf; register int xe,k; if(!r)r= cur_round; zt= funpack(z,&zf,&ze,&zs); if(zs=='-'&&zt!=zro)exceptions|= I_BIT,x= standard_NaN; else switch(zt){ case nan:if(!(z.h&0x80000))exceptions|= I_BIT,z.h|= 0x80000; return z; case inf:case zro:x= z;break; case num:/*92:*/ #line 1750 "./mmix-arith.w" xf.h= 0,xf.l= 2; xe= (ze+0x3fe)>>1; if(ze&1)zf= shift_left(zf,1); rf.h= 0,rf.l= (zf.h>>22)-1; for(k= 53;k;k--){ rf= shift_left(rf,2);xf= shift_left(xf,1); if(k>=43)rf= incr(rf,(zf.h>>(2*(k-43)))&3); else if(k>=27)rf= incr(rf,(zf.l>>(2*(k-27)))&3); if((rf.l> xf.l&&rf.h>=xf.h)||rf.h> xf.h){ xf.l++;rf= ominus(rf,xf);xf.l++; } } if(rf.h||rf.l)xf.l++; return fpack(xf,xe,'+',r); /*:92*/ #line 1736 "./mmix-arith.w" ; } if(zs=='-')x.h|= sign_bit; return x; } /*:91*//*93:*/ #line 1774 "./mmix-arith.w" octa fremstep ARGS((octa,octa,int)); octa fremstep(y,z,delta) octa y,z; int delta; { ftype yt,zt; int ye,ze; char xs,ys,zs; octa x,xf,yf,zf; register int xe,thresh,odd; yt= funpack(y,&yf,&ye,&ys); zt= funpack(z,&zf,&ze,&zs); switch(4*yt+zt){ /*42:*/ #line 731 "./mmix-arith.w" case 4*nan+nan:if(!(y.h&0x80000))exceptions|= I_BIT; case 4*zro+nan:case 4*num+nan:case 4*inf+nan: if(!(z.h&0x80000))exceptions|= I_BIT,z.h|= 0x80000; return z; case 4*nan+zro:case 4*nan+num:case 4*nan+inf: if(!(y.h&0x80000))exceptions|= I_BIT,y.h|= 0x80000; return y; /*:42*/ #line 1788 "./mmix-arith.w" ; case 4*zro+zro:case 4*num+zro:case 4*inf+zro: case 4*inf+num:case 4*inf+inf:x= standard_NaN; exceptions|= I_BIT;break; case 4*zro+num:case 4*zro+inf:case 4*num+inf:return y; case 4*num+num:/*94:*/ #line 1809 "./mmix-arith.w" odd= 0; thresh= ye-delta; if(thresh<ze)thresh= ze; while(ye>=thresh)/*95:*/ #line 1830 "./mmix-arith.w" { if(yf.h==zf.h&&yf.l==zf.l)goto zero_out; if(yf.h<zf.h||(yf.h==zf.h&&yf.l<zf.l)){ if(ye==ze)goto try_complement; ye--,yf= shift_left(yf,1); } yf= ominus(yf,zf); if(ye==ze)odd= 1; while(yf.h<0x400000)ye--,yf= shift_left(yf,1); } /*:95*/ #line 1815 "./mmix-arith.w" ; if(ye>=ze){ exceptions|= E_BIT;return fpack(yf,ye,ys,ROUND_OFF); } if(ye<ze-1)return fpack(yf,ye,ys,ROUND_OFF); yf= shift_right(yf,1,1); try_complement:xf= ominus(zf,yf),xe= ze,xs= '+'+'-'-ys; if(xf.h> yf.h||(xf.h==yf.h&&(xf.l> yf.l||(xf.l==yf.l&&!odd)))) xf= yf,xs= ys; while(xf.h<0x400000)xe--,xf= shift_left(xf,1); return fpack(xf,xe,xs,ROUND_OFF); /*:94*/ #line 1793 "./mmix-arith.w" ; zero_out:x= zero_octa; } if(ys=='-')x.h|= sign_bit; return x; } /*:93*/ #line 41 "./mmix-arith.w" /*:1*/
Go to most recent revision | Compare with Previous | Blame | View Log