← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to mp/mpi.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
mp/mpi.c
 
1
 
 
2
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 
3
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 
4
/*~                                                                     ~*/
 
5
/*~                    OPERATIONS DE BASE (NOYAU)                       ~*/
 
6
/*~                                                                     ~*/
 
7
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 
8
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 
9
/*   This file was modified by W. Schelter to be suitable for optimization
 
10
     and inlining of assembler for maximum speed
 
11
 */
 
12
 
 
13
#include "config.h"
 
14
#include "genpari.h"
 
15
#include "arith.h"
 
16
 
 
17
GEN   mulsi(x,y)
 
18
     plong x;
 
19
     GEN y;
 
20
{ TEMPVARS
 
21
  plong s=signe(y),ly=lgef(y),i;
 
22
  GEN z,zp,yp; ulong hiremainder;
 
23
  
 
24
  if((!x)||(!s)) return gzero;
 
25
  if(x<0) {s= -s;
 
26
           x= -x;
 
27
           if (x < 0) /* -2^31 */
 
28
             {return mulii(stoi(1<<31),y);}
 
29
         }
 
30
  z=cgeti(ly+1);
 
31
  hiremainder=0;
 
32
  MP_START_LOW(yp,y,ly);    MP_START_LOW(zp,z,ly+1);
 
33
  i = MP_COUNT_LG(ly);
 
34
  WHILE_COUNT(--i) { MP_NEXT_UP(zp) = addmul(x,MP_NEXT_UP(yp));}
 
35
  if(hiremainder) {MP_NEXT_UP(zp)=hiremainder;
 
36
                   setlgef(z,ly+1);}
 
37
  else {avma+=4;z[1]=z[0]-1;z++;setlgef(z,ly);}
 
38
  setsigne(z,s);return z;
 
39
}
 
40
 
 
41
int expi(x)
 
42
     GEN x;
 
43
{
 
44
  plong lx=x[1]&0xffff;
 
45
  
 
46
  return lx==2 ? -8388608 : ((lx-2)<<5)-bfffo(x[2])-1;
 
47
}
 
48
 
 
49
GEN addsi(x,y)
 
50
     plong x;
 
51
     GEN y;
 
52
{
 
53
  plong sx,sy,ly,p,i;
 
54
  ulong overflow;  
 
55
  GEN z; TEMPVARS
 
56
 
 
57
  
 
58
  if(!x) return icopy(y);
 
59
  sy=signe(y);if(!sy) return stoi(x);
 
60
  if(x<0) {sx= -1;
 
61
           x= -x;
 
62
           if (x < 0)  /* x=-2^31 */
 
63
             return addii(MOST_NEGS,y);
 
64
         } else sx=1;
 
65
  ly=lgef(y);
 
66
  if(sx==sy)
 
67
    {
 
68
      p=addll(x,y[ly-1]);
 
69
      if(overflow)
 
70
        {
 
71
          z=cgeti(ly+1);z[ly]=p;
 
72
          for(i=ly-1;(i>2)&&(y[i-1]==0xffffffff);i--) z[i]=0;
 
73
          if(i>2)
 
74
            {
 
75
              z[i]=y[i-1]+1;i--;while(i>=3) {z[i]=y[i-1];i--;}
 
76
              z[2]=z[1]=z[0]-1;z++;avma+=4;
 
77
            }
 
78
          else {z[2]=1;z[1]=z[0];}
 
79
        }
 
80
      else
 
81
        {
 
82
          z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
 
83
        }
 
84
      setsigne(z,sx);
 
85
    }
 
86
  else
 
87
    {
 
88
      if(ly==3)
 
89
        {
 
90
          if((ulong)y[2]>(ulong)x)
 
91
            {
 
92
              z=cgeti(3);z[1]=(sy<<24)+3;z[2]=y[2]-x;return z;
 
93
            }
 
94
          if(y[2]==x) return gzero;
 
95
          z=cgeti(3);z[1]=((-sy)<<24)+3;z[2]=x-y[2];return z;
 
96
        }
 
97
      p=subll(y[ly-1],x);
 
98
      if(overflow)
 
99
        {
 
100
          z=cgeti(ly);z[ly-1]=p;
 
101
          for(i=ly-2;!(y[i]);i--) z[i]=0xffffffff;
 
102
          z[i]=y[i]-1;
 
103
          if((i>2)||z[i]) {i--;for(;i>=1;i--) z[i]=y[i];}
 
104
          else
 
105
            {
 
106
              z[2]=z[1]=z[0]-1;z++;avma+=4;setsigne(z,sy);
 
107
            }
 
108
        }
 
109
      else
 
110
        {
 
111
          z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
 
112
        }    
 
113
    }
 
114
  return z;
 
115
}
 
116
 
 
117
GEN addii(x,y)
 
118
     GEN x,y;
 
119
{
 
120
  plong sx,sy,sz,lx,ly,i,j,p;
 
121
  GEN z,xp,yp,zp,xpp,xhigh; TEMPVARS
 
122
  ulong overflow;  
 
123
  lx = lgef(x);
 
124
  ly = lgef(y);
 
125
  if (lx< ly) {z=x;x=y;y=z; sx=lx; lx=ly; ly=sx;} 
 
126
       /* lx>=ly so sx==0 ==> sy==0 */
 
127
  
 
128
  if (0 == (sy=signe(y))) return icopy(x);
 
129
  sx = signe(x);
 
130
 
 
131
  if(sx==sy)
 
132
    {
 
133
      z=cgeti(lx+1);overflow=0;
 
134
      MP_START_LOW(zp,z,lx+1);
 
135
      MP_START_LOW(xp,x,lx);
 
136
      MP_START_LOW(yp,y,ly);
 
137
 
 
138
#ifdef QUICK_LOOP
 
139
       i = ly - 2;
 
140
      QUICK_LOOP(i,ADDXCC);
 
141
#else      
 
142
      i = MP_COUNT_LG(ly);
 
143
      WHILE_COUNT(--i)
 
144
        {ADDLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
 
145
#endif
 
146
 
 
147
     
 
148
      if(overflow)
 
149
        {  GEN xhigh = &MP_HIGH(x,lx);
 
150
         again:    
 
151
          { GEN xpp = &MP_NEXT_UP(xp);
 
152
           if (xpp >= xhigh)
 
153
             { if (*xpp == 0xffffffff)
 
154
                    { MP_NEXT_UP(zp)=0;
 
155
                      goto again;}
 
156
                  
 
157
             else
 
158
               { MP_NEXT_UP(zp) = *xpp + 1;
 
159
                 while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
 
160
                   { MP_NEXT_UP(zp) = *xpp ;}
 
161
                 z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;}}
 
162
           else {z[2]=1;z[1]=x[1]+1;}
 
163
        }}
 
164
      else
 
165
        { j = COUNT(lx - ly);
 
166
            WHILE_COUNT( --j)
 
167
              { MP_NEXT_UP(zp) = MP_NEXT_UP(xp);}
 
168
          z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;
 
169
        }
 
170
    }
 
171
  else
 
172
    {
 
173
      if(lx==ly)
 
174
        /* we have to compare x and y */
 
175
        { j = MP_COUNT_LG(lx);
 
176
          MP_START_HIGH(xp,x,lx);
 
177
          MP_START_HIGH(yp,y,lx);
 
178
          WHILE_COUNT(--j)
 
179
            { ulong tx = MP_NEXT_DOWN(xp);
 
180
              ulong ty = MP_NEXT_DOWN(yp);
 
181
              if ( ty > tx)
 
182
                {z=x;x=y;y=z;sz=sx;sx=sy;sy=sz;
 
183
                 goto DIFFER;}
 
184
              else
 
185
                if ( tx > ty)
 
186
                  {goto DIFFER;}}
 
187
          SAME:  return gzero;  
 
188
        DIFFER:;
 
189
        }
 
190
      
 
191
      z=cgeti(lx);overflow=0;
 
192
      MP_START_LOW(xp,x,lx);MP_START_LOW(yp,y,ly);MP_START_LOW(zp,z,lx);
 
193
      i = MP_COUNT_LG(ly);
 
194
 
 
195
#ifdef QUICK_LOOP
 
196
       i = ly - 2;
 
197
      QUICK_LOOP(i,SUBXCC);
 
198
#else      
 
199
      i = MP_COUNT_LG(ly);
 
200
      WHILE_COUNT(--i)
 
201
        {SUBLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
 
202
#endif
 
203
      
 
204
      if(overflow)
 
205
        { ulong tx ;  
 
206
          while((tx=MP_NEXT_UP(xp)) == 0)
 
207
            MP_NEXT_UP(zp) = 0xffffffff;
 
208
          if (xp >= (xhigh = &MP_HIGH(x,lx)))
 
209
            { MP_NEXT_UP(zp) = tx -1;
 
210
              while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
 
211
                { MP_NEXT_UP(zp) = *xpp;}}
 
212
        }
 
213
      else
 
214
        { i = COUNT(lx - ly);
 
215
          WHILE_COUNT(--i)
 
216
            MP_NEXT_UP(zp) = MP_NEXT_UP(xp);
 
217
        }
 
218
      if(z[2]) z[1]=x[1];
 
219
      else
 
220
        { zp = &z[3];
 
221
          while (*zp ==0){zp++;}  /* x was != y by above */
 
222
          zp -= 2;
 
223
          i = zp - z;
 
224
          zp[1] = (zp[0] = z[0]-i);
 
225
          z = zp;
 
226
          setsigne(z,sx);
 
227
          avma+=(i<<2);
 
228
        }
 
229
    }
 
230
  return z;
 
231
}      
 
232
 
 
233
GEN mulss(x,y)
 
234
     plong x,y;
 
235
{
 
236
  plong s,p1;
 
237
  GEN z; ulong hiremainder;
 
238
  
 
239
  if((!x)||(!y)) return gzero;
 
240
  s=1;
 
241
  if(x<0)
 
242
    {s= -1;
 
243
     x= -x;
 
244
     if (x<0)
 
245
       return mulsi(y,stoi(x));
 
246
   }
 
247
  if(y<0) {s= -s;
 
248
           y= -y;
 
249
           if(y<0)
 
250
             return mulsi((s > 0 ? x : -x),ABS_MOST_NEGS);
 
251
         }
 
252
  p1=mulll(x,y);
 
253
  if(hiremainder) {z=cgeti(4);z[2]=hiremainder;z[3]=p1;}
 
254
  else {z=cgeti(3);z[2]=p1;}
 
255
  z[1]=z[0];setsigne(z,s);return z;
 
256
}
 
257
 
 
258
 
 
259
GEN mulii(x,y)
 
260
     GEN x,y;
 
261
{
 
262
  plong i,j,lx=lgef(x),ly=lgef(y),sx,sy,lz,p1,p2;
 
263
  GEN z; TEMPVARS
 
264
    
 
265
    GEN zz,yy,zp,xx;
 
266
  GEN ylow;
 
267
 ulong hiremainder; 
 
268
  ulong overflow;
 
269
  
 
270
  sx=signe(x);if(!sx) return gzero;
 
271
  sy=signe(y);if(!sy) return gzero;
 
272
  if(sy<0) sx= -sx;
 
273
  if(lx>ly) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;}
 
274
  lz=lx+ly-2;if(lz>=0x10000) err(muler1);
 
275
  z=cgeti(lz);z[1]=z[0];setsigne(z,sx);
 
276
 
 
277
  MP_START_LOW(xx,x,lx);
 
278
  p1 = MP_NEXT_UP(xx);
 
279
 
 
280
  hiremainder=0;
 
281
  i = COUNT(ly-2);
 
282
  MP_START_LOW(yy,y,ly);
 
283
  MP_START_LOW(zz,z,lz);
 
284
 
 
285
  WHILE_COUNT (--i)
 
286
    { MP_NEXT_UP(zz) = addmul(p1,MP_NEXT_UP(yy));}
 
287
 
 
288
  MP_NEXT_UP(zz) = hiremainder;
 
289
 
 
290
  /* restart zz one above bottom */ 
 
291
  MP_START_LOW(zz,z,lz); 
 
292
  
 
293
  MP_START_LOW(ylow,y,ly);
 
294
  ly  = COUNT(ly - MP_CODE_WORDS);
 
295
  lx -= MP_CODE_WORDS;
 
296
  while (--lx > 0)              /* one less iteration first term of x, already used */
 
297
    { plong tem;
 
298
      register plong p11;
 
299
      p11 = MP_NEXT_UP(xx);
 
300
      i = ly;
 
301
      yy = ylow;
 
302
      zp = &MP_NEXT_UP(zz);    /* *zp = second from low word of z first time through */
 
303
      tem = 0;
 
304
      /* ZerO is just a 68k kludge to getit to keep 0 in a reg during this loop*/
 
305
#undef ZERO
 
306
#define ZERO ZerO
 
307
      {  int ZerO = 0; 
 
308
      WHILE_COUNT(--i)
 
309
        { p2 = MP_NEXT_UP(yy);
 
310
          p2 = mulul(p2,p11,hiremainder);
 
311
          MP_NEXT_UP(zp);
 
312
          p2 = add_carry(p2,*zp,hiremainder);
 
313
          p2 = add_carry(p2,tem,hiremainder);
 
314
          *zp = p2;
 
315
          tem = hiremainder;
 
316
        }
 
317
       }
 
318
      MP_NEXT_UP(zp) = hiremainder;
 
319
 
 
320
#undef ZERO
 
321
#define ZERO 0
 
322
 
 
323
    }
 
324
  if(!MP_HIGH(z,lz))
 
325
    {                           /* shift header one along  decreasing lg and lgef */
 
326
      z[2]=z[1]-1;z[1]=z[0]-1;z++;avma+=4;
 
327
    }
 
328
  return z;
 
329
}
 
330
 
 
331
GEN confrac(x)
 
332
     GEN x;
 
333
{
 
334
  plong lx=lg(x),ex= -expo(x)-1,ex1,av=avma,ly,ey;
 
335
  plong lr,nbdec,k,i,j; ulong hiremainder;
 
336
  GEN y,res; TEMPVARS
 
337
  
 
338
  ey=((lx-2)<<5)+ex;ly=(ey+63)>>5;y=cgeti(ly);ex1=ex>>5; /* 95 dans mp.s faux? */
 
339
  for(i=0;i<ex1;i++) y[i]=0;
 
340
  ex&=31;
 
341
  if(!ex) for(j=2;j<lx;j++) y[i++]=x[j];
 
342
  else
 
343
    {
 
344
      k=0;
 
345
      for(j=2;j<lx;j++) {y[i++]=shiftlr(x[j],ex)+k;k=hiremainder;}
 
346
      y[ly-2]=k;
 
347
    }
 
348
  y[ly-1]=0;
 
349
  nbdec=ey*0.30103+1;lr=(nbdec+17)/9;res=cgeti(lr);
 
350
  *res=nbdec;
 
351
  for(j=1;j<lr;j++)
 
352
    {
 
353
      hiremainder=0;
 
354
      for(i=ly-1;i>=0;i--) y[i]=addmul(y[i],1000000000);
 
355
      res[j]=hiremainder;
 
356
    }
 
357
  avma=av;return res;
 
358
}
 
359
 
 
360
/* x/y : uses hiremainder for return */
 
361
GEN divss(x,y)
 
362
     plong x,y;
 
363
{
 
364
  plong p1; 
 
365
  
 
366
  if(!y) err(diver1);
 
367
  if (x == (1<<31))
 
368
    return divis(stoi(x),y);
 
369
  hiremainder=0;p1=divll((ulong)abs(x),(ulong)abs(y));
 
370
  if(y<0) {hiremainder= -((plong)hiremainder);p1= -p1;}
 
371
  if(x<0) p1= -p1;
 
372
  return stoi(p1);
 
373
}
 
374
 
 
375
GEN modss(x,y)
 
376
     plong x,y;
 
377
{
 
378
  plong y1; ulong hiremainder;
 
379
  
 
380
  if(!y) err(moder1);
 
381
  if (x == (1<<31))
 
382
    return modis(stoi(x),y);
 
383
  hiremainder=0;divll(abs(x),y1=abs(y));
 
384
  if(!hiremainder) return gzero;
 
385
  return (((plong)hiremainder)<0) ? stoi(y1-hiremainder) : stoi(hiremainder);
 
386
}
 
387
 
 
388
GEN resss(x,y)
 
389
     plong x,y;
 
390
{ ulong hiremainder;
 
391
  if(!y) err(reser1);
 
392
  hiremainder=0;divll(abs(x),abs(y));
 
393
  return (y<0) ? stoi(-((plong)hiremainder)) : stoi(hiremainder);
 
394
}
 
395
 
 
396
/* uses hiremainder for return */
 
397
GEN divsi(x,y)
 
398
     plong x;
 
399
     GEN y;
 
400
{ 
 
401
  plong s=signe(y),ly=lgef(y),p1;
 
402
 
 
403
  if(!s) err(diver2);
 
404
  if((!x)||(ly>3)||(y[2]<0)) {hiremainder=x;return gzero;}
 
405
  if (x== 1<<31)
 
406
    return divii(stoi(x),y);
 
407
  hiremainder=0;p1=divll(abs(x),y[2]);
 
408
  if(signe(y)<0) {hiremainder= -((plong)hiremainder);p1= -p1;}
 
409
  if(x<0) p1= -p1;
 
410
  return stoi(p1);
 
411
}
 
412
 
 
413
/* this uses the GLOBAL hiremainder to return its remainder
 
414
   We cannot make it a local.
 
415
 */
 
416
GEN divis(y,x)
 
417
     plong x;
 
418
     GEN y;
 
419
{ ulong hi;
 
420
  plong s=signe(y),ly=lgef(y),i,d;
 
421
  GEN z;
 
422
  
 
423
  if(!x) err(diver4);
 
424
  if(!s) {hiremainder=0;return gzero;}
 
425
  if(x<0) {s= -s;x= -x;
 
426
           if (x < 0)
 
427
             return divii(y,stoi(x));
 
428
         } 
 
429
  if((ulong)x>(ulong)y[2])
 
430
    {
 
431
      if(ly==3) {hiremainder=itos(y);return gzero;}
 
432
      else {z=cgeti(ly-1);d=1;hi=y[2];}
 
433
    }
 
434
  else {z=cgeti(ly);d=0;hi=0;}
 
435
  for(i=d+2;i<ly;i++) z[i-d]=divul(y[i],x,hi);
 
436
  z[1]=z[0];setsigne(z,s);
 
437
  hiremainder= (s < 0 ? -((plong)hi) : hi) ;
 
438
  return z;
 
439
}
 
440
 
 
441
/*
 
442
#       entree : x pointe sur i2 de type I (dividende)              
 
443
#                y pointe sur i1 de type I (diviseur)               
 
444
#                z  contient un pointeur sur le reste si l'on       
 
445
#                veut a la fois q et r, 0 si l'on ne veut que le    
 
446
#                quotient, -1 si l'on ne veut que le reste          
 
447
#       sortie : resultat est  q si celui-ci est attendu, et sinon r 
 
448
#                r. z  pointe sur r si q et r sont attendus
 
449
#                (toutes les zones sont creees)                     
 
450
#       remarque : il s'agit de la 'fausse division' ; le reste est 
 
451
#                 du signe du dividende                             
 
452
*/
 
453
 
 
454
GEN dvmdii(x,y,z)
 
455
     GEN x,y,*z;
 
456
{
 
457
  plong av,lx,ly,lz,i,j,dec,sh,k,k1,sx=signe(x),sy=signe(y);
 
458
  plong saux,k3,k4,av1,flk4;
 
459
  ulong si,qp;
 
460
  GEN p1,p2,p3,p4,xp,p1p,p2p,pp; TEMPVARS
 
461
  ulong hiremainder,overflow;
 
462
  
 
463
  if(!sy) err(dvmer1);
 
464
  if(!sx)
 
465
    {
 
466
      if(((plong)z==0xffffffff)||((plong)z==0)) return gzero;
 
467
      *z=gzero;return gzero;
 
468
    }
 
469
  lx=lgef(x);ly=lgef(y);lz=lx-ly;
 
470
  if(lz<0)
 
471
    {
 
472
      if((plong)z==0xffffffff) return icopy(x);
 
473
      if(z==0) return gzero;
 
474
      *z=icopy(x);return gzero;
 
475
    }
 
476
  av=avma;if(sx<0) sy= -sy;
 
477
  if(ly==3)
 
478
    { int lgp1;
 
479
      si=y[2];
 
480
      MP_START_HIGH(xp,x,lx);
 
481
      if (si > (ulong)MP_HIGH(x,lx))
 
482
        { lgp1=lx-1; hiremainder= MP_NEXT_DOWN(xp);}
 
483
      else { lgp1=lx; hiremainder=0;}
 
484
      p1 = cgeti(lgp1); i = MP_COUNT_LG(lgp1);
 
485
      MP_START_HIGH(p1p,p1,lgp1);
 
486
      WHILE_COUNT(--i) { MP_NEXT_DOWN(p1p) = divll(MP_NEXT_DOWN(xp),si);}
 
487
 
 
488
      if((plong)z==0xffffffff)
 
489
        {
 
490
          avma=av;if(!hiremainder) return gzero;
 
491
          p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;return p2;
 
492
        }
 
493
      if(lgp1!= 2) {p1[1]=p1[0];setsigne(p1,sy);} else {avma=av;p1=gzero;}
 
494
      if(z==0) return p1;
 
495
      if(!hiremainder) *z=gzero;
 
496
      else {p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;*z=p2;}
 
497
      return p1;
 
498
    }
 
499
  else
 
500
    {
 
501
      p1=cgeti(lx);
 
502
      sh=bfffo(y[2]);
 
503
      if(sh)
 
504
        { GEN p2p,yp;
 
505
          MP_START_HIGH(yp,y,ly);
 
506
          p2=cgeti(ly);
 
507
          k=shiftl(MP_NEXT_DOWN(yp),sh);
 
508
          MP_START_HIGH(p2p,p2,ly);
 
509
          i = MP_COUNT_LG(ly-1);
 
510
          WHILE_COUNT(--i)
 
511
            {
 
512
              k1=shiftl(MP_NEXT_DOWN(yp),sh);
 
513
              MP_NEXT_DOWN(p2p) = k + hiremainder;
 
514
              k = k1;
 
515
            }
 
516
          MP_NEXT_DOWN(p2p) = k ; k=0;
 
517
 
 
518
          MP_START_HIGH(xp,x,lx);
 
519
          MP_START_HIGH(p1p,p1,lx);
 
520
          MP_NEXT_UP(p1p) ;    /* yes go out of range !! */
 
521
          i = MP_COUNT_LG(lx);
 
522
          WHILE_COUNT (--i)
 
523
            { k1 = shiftl(MP_NEXT_DOWN(xp),sh);
 
524
              MP_NEXT_DOWN(p1p) = k + hiremainder; k = k1;
 
525
            }
 
526
          MP_NEXT_DOWN(p1p) = k;
 
527
        }
 
528
      else {
 
529
        MP_START_HIGH(xp,x,lx);
 
530
        MP_START_HIGH(p1p,p1,lx);
 
531
        MP_NEXT_UP(p1p) ;    /* yes go out of range !! */
 
532
        MP_NEXT_DOWN(p1p) = 0;
 
533
        j = MP_COUNT_LG(lx);
 
534
        WHILE_COUNT (-- j)
 
535
          { MP_NEXT_DOWN(p1p) = MP_NEXT_DOWN(xp);}
 
536
        p2 = y;}
 
537
      si=p2[2];saux=p2[3];
 
538
      MP_START_HIGH(p1p,p1,lx); MP_NEXT_UP(p1p) ; /* out of bound */
 
539
      i = COUNT(lz+1);
 
540
      WHILE_COUNT(--i)       
 
541
 
 
542
        { GEN pp;
 
543
          if(MP_NEXT_DOWN(p1p)==si) 
 
544
            {  /* Using fact that next_down does post increment */
 
545
            qp=0xffffffff;k=addll(si,*p1p);    
 
546
                                                
 
547
            }
 
548
          else
 
549
            {
 
550
              hiremainder=p1p[-1];qp=divll(*p1p,si);
 
551
              overflow=0;k=hiremainder;
 
552
            }
 
553
          if(!overflow)
 
554
            {
 
555
/*            k1=mulll(qp,saux);k3=subll(k1,p1p[1]);k+=overflow;
 
556
              flk4=((ulong)hiremainder>(ulong)k);k4=subll(hiremainder,k);
 
557
              while(flk4) {qp--;k3=subll(k3,saux);
 
558
                           k4-=overflow;flk4=((ulong)k4>(ulong)si);
 
559
                           k4=subll(k4,si);}
 
560
*/
 
561
              k1=mulll(qp,saux);k3=subll(k1,p1p[1]);
 
562
              k4=subllx(hiremainder,k);
 
563
              while((!overflow)&&k4) {qp--;k3=subll(k3,saux);k4=subllx(k4,si);}
 
564
              
 
565
            }
 
566
          hiremainder=0;
 
567
 
 
568
          j = MP_COUNT_LG(ly);
 
569
          MP_START_LOW(pp,p1p,ly-2);
 
570
          MP_START_LOW(p2p,p2,ly);
 
571
          WHILE_COUNT(--j)        
 
572
            { GEN ppp;
 
573
              k1=addmul(qp,MP_NEXT_UP(p2p));
 
574
              ppp = &MP_NEXT_UP(pp);
 
575
              *ppp =subll(*ppp,k1);hiremainder+=overflow;
 
576
            }
 
577
          if((ulong)p1p[-1]<(ulong)hiremainder)
 
578
            {
 
579
              overflow=0;qp--;
 
580
              j = MP_COUNT_LG(ly);
 
581
              MP_START_LOW(pp,p1p,ly-2);
 
582
              MP_START_LOW(p2p,p2,ly);
 
583
              WHILE_COUNT(--j){ GEN ppp = &MP_NEXT_UP(pp);
 
584
                                ADDLLX(*ppp,MP_NEXT_UP(p2p),*ppp);}
 
585
 
 
586
            }
 
587
          p1p[-1] = qp;
 
588
        }
 
589
      av1=avma;
 
590
      if((plong)z!=0xffffffff)
 
591
        {ulong lgp3 = lz + 2;
 
592
         MP_START_LOW(p1p,p1,lgp3);
 
593
         if (p1[1]) {lgp3++;}
 
594
           else if (lz==0) sy=0;
 
595
         p3 = cgeti(lgp3);
 
596
         MP_START_LOW(pp,p3,lgp3);
 
597
         j = MP_COUNT_LG(lgp3);
 
598
         WHILE_COUNT(--j)
 
599
           {MP_NEXT_UP(pp) = MP_NEXT_UP(p1p) ;}
 
600
         if(lgp3<3) {p3[1]=2;} else {p3[1]=p3[0];setsigne(p3,sy);}
 
601
        }
 
602
      if(z!=0)
 
603
        {
 
604
          for(j=lz+2;(j<lx)&&(!p1[j]);j++);
 
605
          if(j==lx) p4=icopy(gzero);
 
606
          else
 
607
            {
 
608
              p4=cgeti(lx-j+2);p4[1]=p4[0];
 
609
              if(!sh) for(i=2;j<lx;j++,i++) p4[i]=p1[j];
 
610
              else
 
611
                {
 
612
                  hiremainder=0;k1=shiftlr(p1[j++],sh);k=hiremainder;
 
613
                  if(k1) {p4[2]=k1;dec=1;} else
 
614
                    {p4[1]=p4[0]-1;p4++;avma+=4;p4[1]=p4[0];dec=0;}
 
615
                  for(i=2+dec;j<lx;j++,i++)
 
616
                    {
 
617
                      p4[i]=shiftlr(p1[j],sh)+k;k=hiremainder;
 
618
                    }
 
619
                }
 
620
              setsigne(p4,sx);
 
621
            }
 
622
        }
 
623
      if((plong)z==0xffffffff) return gerepile(av,av1,p4);
 
624
      if((plong)z==0) return gerepile(av,av1,p3);
 
625
      dec=lpile(av,av1,0)>>2;*z=p4+dec;return p3+dec;     
 
626
    }
 
627
}
 
628
/* machines which provide an inline version of mulul need
 
629
 to provide a function for calls where that inlining can't take place
 
630
 */
 
631
#ifdef NEED_MULUL3
 
632
ulong
 
633
mulul3(a,b,c)
 
634
   ulong a,b,*c;
 
635
{ return mulul(a,b,*c);}
 
636
#endif
 
637
 
 
638
#ifdef NEED_DIVUL3
 
639
ulong
 
640
divul3(a,b,c)
 
641
   ulong a,b,*c;
 
642
{ return divul(a,b,*c);}
 
643
#endif
 
644
 
 
645
/*
 
646
;;- Local variables:
 
647
;;- version-control:t
 
648
;;- End:
 
649
*/