home *** CD-ROM | disk | FTP | other *** search
- #include "genpari.h"
-
- void printversion()
- {
- pariputs(" GP/PARI CALCULATOR Version 1.37\n");
- pariputs(" (portable 32-bit version)\n");
- }
-
- ulong overflow,hiremainder;
-
- int addll(x,y)
- int x,y;
- {
- int z;
-
- z=x+y;
- if((x>=0)&&(y>=0)) {overflow=0;return z;}
- if((x<0)&&(y<0)) {overflow=1;return z;}
- overflow=(z>=0)?1:0; return z;
- }
-
- int addllx(x,y)
- int x,y;
- {
- int z;
-
- z=x+y+overflow;
- if((x>=0)&&(y>=0)) {overflow=0;return z;}
- if((x<0)&&(y<0)) {overflow=1;return z;}
- overflow=(z>=0)?1:0; return z;
- }
-
- int subll(x,y)
- int x,y;
- {
- int z;
-
- z=x-y;
- if((x>=0)&&(y<0)) {overflow=1;return z;}
- if((x<0)&&(y>=0)) {overflow=0;return z;}
- overflow=(z>=0)?0:1; return z;
- }
-
- int subllx(x,y)
- int x,y;
- {
- int z;
-
- z=x-y-overflow;
- if((x>=0)&&(y<0)) {overflow=1;return z;}
- if((x<0)&&(y>=0)) {overflow=0;return z;}
- overflow=(z>=0)?0:1; return z;
- }
-
- int shiftl(x,y)
- ulong x,y;
- {
- hiremainder=x>>(32-y);return (x<<y);
- }
-
- int shiftlr(x,y)
- ulong x,y;
- {
- hiremainder=x<<(32-y);return (x>>y);
- }
-
- int bfffo(x)
- ulong x;
- {
- int sc;
- static int tabshi[16]={4,3,2,2,1,1,1,1,0,0,0,0,0,0,0,0};
-
- if(x&(0xffff0000)) sc=0;else {sc=16;x<<=16;}
- if(!(x&(0xff000000))) {sc+=8;x<<=8;}
- if(x&(0xf0000000)) x>>=28;else {sc+=4;x>>=24;}
- sc+=tabshi[x];return sc;
- }
-
- int mulll(x,y)
- ulong x,y;
- {
- ulong xlo,xhi,ylo,yhi;
- ulong z;
-
- xlo=x&65535;xhi=x>>16;ylo=y&65535;yhi=y>>16;
- z=addll(xlo*yhi,xhi*ylo);
- hiremainder=(overflow)?xhi*yhi+65536+(z>>16):xhi*yhi+(z>>16);
- z=addll(xlo*ylo,(z<<16));hiremainder+=overflow;
- return z;
- }
-
- int addmul(x,y)
- ulong x,y;
- {
- ulong xlo,xhi,ylo,yhi;
- ulong z,z2;
-
- xlo=x&65535;xhi=x>>16;ylo=y&65535;yhi=y>>16;
- z=addll(xlo*yhi,xhi*ylo);
- z2=(overflow)?xhi*yhi+65536+(z>>16):xhi*yhi+(z>>16);
- z=addll(xlo*ylo,(z<<16));z2+=overflow;
- z=addll(z,hiremainder);hiremainder=z2+overflow;
- return z;
- }
-
- /* Ancienne version de divll, en generale plus lente sauf si le processeur flottant
- est tres rapide
-
- int divll(x,y)
- ulong x,y;
- {
- ulong p,p1,p2;
- long p3,p4;
- double dbl;
-
- if((ulong)hiremainder>=y) err(divller1);
- p1=hiremainder;dbl=(C4*p1+x)/y;
- if(dbl < C4/2) p=(ulong)dbl;
- else
- {
- if(dbl>C4-1) p=0xffffffff;
- else{dbl-=(C4/2);p=(ulong)dbl;p|=0x80000000;}
- }
- p2=mulll(p,y);p3=x-p2;
- if(x<p2) hiremainder++;
- if((p1==(ulong)hiremainder)&&((ulong)p3<y))
- {hiremainder=p3;return p;}
- if(p1<(ulong)hiremainder) {hiremainder=p3+y;return p-1;}
- hiremainder=p3-y;return p+1;
- }
-
- */
-
- int divll(x,y)
- ulong x,y;
- {
- #define HIBIT 0x80000000
- #define HIMASK 0xffff0000
- #define LOMASK 0xffff
- #define HIWORD(a) (a >> 16)
- /* si le compilateur est bugge, il faut mettre (a >> 16) & LOMASK) */
- #define LOWORD(a) (a & LOMASK)
- #define GLUE(hi, lo) ((hi << 16) + lo)
- #define SPLIT(a, b, c) b = HIWORD(a); c = LOWORD(a)
-
- ulong v1, v2, u3, u4, q1, q2, aux, aux1, aux2;
- int k;
-
- for(k = 0; !(y & HIBIT); k++)
- {
- hiremainder <<= 1;
- if (x & HIBIT) hiremainder++;
- x <<= 1;
- y <<= 1;
- }
-
- SPLIT(y, v1, v2);
- SPLIT(x, u3, u4);
-
- q1 = hiremainder / v1; if (q1 & HIMASK) q1 = LOMASK;
- hiremainder -= q1 * v1;
- aux = v2 * q1;
- again:
- SPLIT(aux, aux1, aux2);
- if (aux2 > u3) aux1++;
- if (aux1 > hiremainder) {q1--; hiremainder += v1; aux -= v2; goto again;}
- u3 -= aux2;
- hiremainder -= aux1;
- hiremainder <<= 16; hiremainder += u3 & LOMASK;
-
- q2 = hiremainder / v1; if (q2 & HIMASK) q2 = LOMASK;
- hiremainder -= q2 * v1;
- aux = v2 * q2;
- again2:
- SPLIT(aux, aux1, aux2);
- if (aux2 > u4) aux1++;
- if (aux1 > hiremainder) {q2--; hiremainder += v1; aux -= v2; goto again2;}
- u4 -= aux2;
- hiremainder -= aux1;
- hiremainder <<= 16; hiremainder += u4 & LOMASK;
- hiremainder >>= k;
- return GLUE(q1, q2);
- }
-
-
-
- long mulmodll(a,b,c)
- ulong a,b,c;
- {
- divll(mulll(a,b),c);
- return hiremainder;
- }
-