← Back to branch summary

~ubuntu-branches/ubuntu/vivid/gcl/vivid

« back to all changes in this revision

Viewing changes to o/pari_big.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2002-03-04 14:29:59 UTC
  • Revision ID: james.westby@ubuntu.com-20020304142959-dey14w08kr7lldu3
Tags: upstream-2.5.0.cvs20020219
Import upstream version 2.5.0.cvs20020219

Show diffs side-by-side

added added

removed removed

Lines of Context:
o/pari_big.c
 
1
  /* Copyright William F. Schelter 1991
 
2
   Bignum routines.
 
3
 
 
4
 
 
5
   
 
6
num_arith.c: add_int_big
 
7
num_arith.c: big_minus
 
8
num_arith.c: big_plus
 
9
num_arith.c: big_quotient_remainder
 
10
num_arith.c: big_sign
 
11
num_arith.c: big_times
 
12
num_arith.c: complement_big
 
13
num_arith.c: copy_big
 
14
num_arith.c: div_int_big
 
15
num_arith.c: mul_int_big
 
16
num_arith.c: normalize_big
 
17
num_arith.c: normalize_big_to_object
 
18
num_arith.c: stretch_big
 
19
num_arith.c: sub_int_big
 
20
num_comp.c: big_compare
 
21
num_comp.c: big_sign
 
22
num_log.c: big_sign
 
23
num_log.c: copy_to_big
 
24
num_log.c: normalize_big
 
25
num_log.c: normalize_big_to_object
 
26
num_log.c: stretch_big
 
27
num_pred.c: big_sign
 
28
number.c: big_to_double
 
29
predicate.c: big_compare
 
30
typespec.c: big_sign
 
31
print.d: big_minus
 
32
print.d: big_sign
 
33
print.d: big_zerop
 
34
print.d: copy_big
 
35
print.d: div_int_big
 
36
read.d: add_int_big
 
37
read.d: big_to_double
 
38
read.d: complement_big
 
39
read.d: mul_int_big
 
40
read.d: normalize_big
 
41
read.d: normalize_big_to_object
 
42
 
 
43
 */
 
44
 
 
45
#define BCOPY_BODY(x,y) \
 
46
do { int *ucop = (int *)(x); \
 
47
    int *vcop = (int *) (y); \
 
48
  {int j = lgef(ucop); \
 
49
    while(--j >= 0) \
 
50
      { *vcop++ = *ucop++;}}}while (0)
 
51
 
 
52
bcopy_body(x,y)
 
53
    GEN x,y;
 
54
{BCOPY_BODY(x,y);}
 
55
 
 
56
 
 
57
 
 
58
/* make a bignum with (most <<32 + least) */
 
59
object
 
60
bignum2(most, least)
 
61
int most, least;
 
62
{ static  plong u [4] 
 
63
   = {0x01010004 ,0x01010004, 0,0};
 
64
  GEN w;
 
65
  int l;
 
66
  if(most) {setlgef(u,4),l=4;}
 
67
    else {l=3; setlgef(u,3);}
 
68
  MP_START_LOW(w,u,l);
 
69
  MP_NEXT_UP(w) = least;
 
70
  if (most) MP_NEXT_UP(w) = most;
 
71
 return  make_integer(u);
 
72
}
 
73
 
 
74
 
 
75
 
 
76
/* coerce a pari GEN to a bignum or fixnum */
 
77
 
 
78
object
 
79
make_integer(u)
 
80
GEN u;
 
81
{ int l = lgef(u);
 
82
  if (l > (MP_CODE_WORDS+1) ||
 
83
      ( l == (MP_CODE_WORDS+1)  &&
 
84
       (MP_ONLY_WORD(u) & (1<<31)) != 0
 
85
       && (MP_ONLY_WORD(u) == ( 1<<31) ? signe(u) > 0 : 1)))
 
86
    { object ans ;
 
87
      GEN w;
 
88
 
 
89
      { BEGIN_NO_INTERRUPT;
 
90
      big_register_1->big.big_length = lg(u);
 
91
      big_register_1->big.big_self =  u;
 
92
      ans = alloc_object(t_bignum);
 
93
      ans->big.big_self = 0;
 
94
      w = (plong *)alloc_relblock(lg(u)*sizeof(plong));
 
95
      /* may have been relocated */
 
96
      u = (GEN)   big_register_1->big.big_self ;
 
97
      ans->big.big_self = w;
 
98
      ans->big.big_length = l;
 
99
      BCOPY_BODY(u , w);
 
100
      setlg(w,l);
 
101
      END_NO_INTERRUPT;}        
 
102
      return ans;
 
103
    }
 
104
  else
 
105
    if (signe(u) > 0) return make_fixnum(MP_ONLY_WORD(u));
 
106
  else
 
107
    if (signe(u) < 0) return make_fixnum(-MP_ONLY_WORD(u));
 
108
  else
 
109
    return(small_fixnum(0));
 
110
 }
 
111
 
 
112
 
 
113
object
 
114
make_bignum(u)
 
115
GEN u;
 
116
{  BEGIN_NO_INTERRUPT;
 
117
    { object ans = alloc_object(t_bignum);
 
118
      GEN w;
 
119
      ans->big.big_length = lg(u);
 
120
      /* save u */
 
121
      ans->big.big_self = u;
 
122
      w = (plong *)alloc_relblock(lg(u)*sizeof(plong));
 
123
      /* restore  u */
 
124
      u = ans->big.big_self ;
 
125
      ans->big.big_self = w;
 
126
      BCOPY_BODY(u ,  ans->big.big_self);
 
127
      END_NO_INTERRUPT;
 
128
      return ans;
 
129
     }}
 
130
 
 
131
big_zerop(x)
 
132
 object x;
 
133
{ return (signe(MP(x))== 0);}
 
134
 
 
135
big_compare(x, y)
 
136
     object x,y;
 
137
{return cmpii(MP(x),MP(y));}
 
138
 
 
139
object
 
140
big_minus(x)
 
141
     object x;
 
142
{ object y; 
 
143
  setsigne(MP(x),-(signe(MP(x))));
 
144
  y = make_integer(MP(x));
 
145
  setsigne(MP(x),-(signe(MP(x))));
 
146
  return  y;
 
147
}
 
148
 
 
149
gcopy_to_big(res,x)
 
150
     GEN res;
 
151
     object x;
 
152
{int l = (x)->big.big_length;
 
153
 int lgres = lg(res);
 
154
 if (l< lgres)
 
155
    { BEGIN_NO_INTERRUPT;
 
156
      big_register_1->big.big_length = lgres;
 
157
      big_register_1->big.big_self = res;
 
158
      (x)->big.big_self = (GEN) alloc_relblock(lgres*sizeof(int));
 
159
      (x)->big.big_length = lgres; 
 
160
      res =    big_register_1->big.big_self ;
 
161
      END_NO_INTERRUPT;
 
162
    }
 
163
 BCOPY_BODY(res,(x)->big.big_self);
 
164
  if (l>lgres)
 
165
    { setlg((x)->big.big_self, l);}
 
166
} 
 
167
            
 
168
 
 
169
add_int_big(i, x)
 
170
int i;
 
171
object x;
 
172
{
 
173
       MPOP_DEST(x,addsi,i,MP(x));
 
174
}
 
175
 
 
176
sub_int_big(i, x)
 
177
int i;
 
178
object x;
 
179
{ MPOP_DEST(x,subsi,i,MP(x));
 
180
}
 
181
 
 
182
mul_int_big(i, x)
 
183
int i;
 
184
object x;
 
185
{ MPOP_DEST(x,mulsi,i,MP(x));
 
186
}    
 
187
 
 
188
/*
 
189
        Div_int_big(i, x) destructively divides non-negative bignum x
 
190
        by positive int i.
 
191
        X will hold the quotient from  the division.
 
192
        Div_int_big(i, x) returns the remainder of the division.
 
193
        I should be positive.
 
194
        X should be non-negative.
 
195
*/
 
196
 
 
197
div_int_big(i, x)
 
198
int i;
 
199
object x;
 
200
{ save_avma;
 
201
  GEN res = divis(MP(x),i);
 
202
  gcopy_to_big(res,x);
 
203
  restore_avma;
 
204
  return hiremainder;
 
205
}
 
206
 
 
207
 
 
208
object
 
209
big_plus(x, y)
 
210
object x,y;
 
211
{ MPOP(return,addii,MP(x),MP(y));
 
212
}
 
213
 
 
214
object
 
215
big_times(x, y)
 
216
object x,y;
 
217
{
 
218
  MPOP(return,mulii,MP(x),MP(y));
 
219
}
 
220
 
 
221
 
 
222
 
 
223
big_quotient_remainder(x0, y0, qp, rp)
 
224
     object x0,y0,*qp,*rp;
 
225
{
 
226
  GEN res,quot;
 
227
  save_avma;
 
228
  res = dvmdii(MP(x0),MP(y0),&quot);
 
229
  *qp = make_integer(res);
 
230
  *rp = make_integer(quot);
 
231
  restore_avma;
 
232
  return;
 
233
 
 
234
}
 
235
 
 
236
        
 
237
double
 
238
big_to_double(x)
 
239
     object x;
 
240
{
 
241
        double d, e;
 
242
        GEN u = MP(x);
 
243
        unsigned int *w;
 
244
        int l;
 
245
        e =  4.294967296e9;
 
246
 
 
247
        l = lgef(u);
 
248
        MP_START_HIGH(w,(unsigned int *) u,l);
 
249
        l = l - MP_CODE_WORDS;
 
250
 
 
251
        if (l == 0) return 0.0;
 
252
 
 
253
        d = (double) MP_NEXT_DOWN(w);
 
254
        while (--l > 0)
 
255
          {d = e*d + (double)(MP_NEXT_DOWN(w));}
 
256
        if (signe(u)>0) return d;
 
257
          else return -d;
 
258
      }
 
259
        
 
260
 
 
261
object
 
262
normalize_big_to_object(x)
 
263
 object x;
 
264
{ return make_integer(MP(x));}
 
265
  
 
266
 
 
267
object copy_big(x)
 
268
     object x;
 
269
{
 
270
  if (type_of(x)==t_bignum)
 
271
    return make_bignum(MP(x));
 
272
  else FEerror("bignum expected",0);
 
273
 
 
274
}
 
275
 
 
276
 
 
277
object
 
278
copy_to_big(x)
 
279
     object x;
 
280
{object y;
 
281
 
 
282
        if (type_of(x) == t_fixnum) {
 
283
          save_avma;
 
284
          y = make_bignum(stoi(fix(x)));
 
285
          restore_avma;
 
286
        } else if (type_of(x) == t_bignum)
 
287
                y = copy_big(x);
 
288
        else
 
289
                FEerror("integer expected",0);
 
290
        return(y);
 
291
      }
 
292
  
 
293
 
 
294
/* return the power of x */
 
295
GEN
 
296
powerii(x,y)
 
297
     GEN x, y;
 
298
{  GEN ans = gun;
 
299
   if (signe(y) < 0) FEerror("bad",0);
 
300
   while (lgef(y) > 2){
 
301
     if (MP_LOW(y,lgef(y)) & 1)
 
302
       { ans = mulii(ans,x);}
 
303
     x = mulii(x,x);
 
304
     y = shifti(y,-1);}
 
305
   return ans;
 
306
 }
 
307
 
 
308
 
 
309
 
 
310
replace_copy1(x,y)
 
311
     GEN y,x;
 
312
{ int j = lgef(x);
 
313
 if (y && j <= lg(y))
 
314
    { x++; y++;
 
315
      while (--j >0)
 
316
      {*y++ = *x++;}
 
317
     return 0;}
 
318
 END:
 
319
 return j*2*sizeof(GEN);
 
320
}
 
321
 
 
322
/* doubles the length ! */
 
323
GEN
 
324
replace_copy2(x,y)
 
325
     GEN y,x;
 
326
{GEN yp = y;  
 
327
 int k,j = lgef(x);
 
328
 k = j;
 
329
 while (--j >=0)
 
330
   {*yp++ = *x++;}
 
331
 y[0] = INT_FLAG + k*2;
 
332
 return y;}
 
333
 
 
334
#define STOI(x,y) do{ \
 
335
  if (x ==0) { y[1]=2;} \
 
336
  else if((x)>0) {y[1]=0x1000003;y[2]=x;} \
 
337
                  else{y[1]=0xff000003;y[2]= -x;}}while (0)
 
338
 
 
339
/* actually y == 0 is not supposed to happen !*/
 
340
                    
 
341
obj_replace_copy1(x,y)
 
342
     object x;
 
343
     GEN y;
 
344
{ int j ;
 
345
  GEN xp;
 
346
  { if (type_of(x) == t_bignum)
 
347
      {   j = lgef(MP(x));
 
348
          if (y && j <= lg(y))
 
349
            { xp=MP(x);
 
350
              xp++; y++;
 
351
              while (--j >0)
 
352
                {*y++ = *xp++;}
 
353
              return 0;}}
 
354
  else
 
355
    { if (y==0) return 3*2*sizeof(GEN) ;
 
356
      STOI(fix(x),y); return 0;}}
 
357
 END:
 
358
 return j*2*sizeof(GEN);
 
359
}
 
360
 
 
361
/* doubles the length ! */
 
362
GEN
 
363
obj_replace_copy2(x,y)
 
364
     object x;
 
365
     GEN y;
 
366
{GEN yp = y;
 
367
 GEN xp;
 
368
 int k,j;
 
369
 if (type_of(x) == t_bignum)
 
370
   { j = lgef(MP(x));
 
371
     k = j;
 
372
     xp = MP(x);
 
373
     while (--j >=0)
 
374
       {*yp++ = *xp++;}
 
375
     y[0] = INT_FLAG + k*2;}
 
376
 else  {STOI(fix(x),yp); y[0] = INT_FLAG+3*2;}
 
377
 return y;}
 
378
 
 
379
 
 
380
GEN
 
381
otoi(x)
 
382
     object x;
 
383
{if (type_of(x)==t_fixnum) return stoi(fix(x));
 
384
 if (type_of(x)==t_bignum)
 
385
   return (MP(x));
 
386
 FEwrong_type_argument(sLinteger,x);
 
387
 return 0;
 
388
}
 
389
 
 
390
object
 
391
alloc_bignum_static(len)
 
392
int len;
 
393
    { object ans = alloc_object(t_bignum);
 
394
      GEN w;
 
395
      ans->big.big_length = len;
 
396
      ans->big.big_self = 0;
 
397
      w = (GEN)AR_ALLOC(alloc_contblock,len,unsigned plong);
 
398
      ans->big.big_self = w;
 
399
      w[0] = INT_FLAG + len;
 
400
      return ans;
 
401
     }
 
402
 
 
403
 
 
404
GEN
 
405
setq_io(x,all,val)
 
406
     GEN x;
 
407
     object val;
 
408
     object *all;
 
409
{int n= obj_replace_copy1(val,x);
 
410
 if (n)
 
411
   { *all = alloc_bignum_static(n/sizeof(int));
 
412
     return obj_replace_copy2(val,MP(*all));
 
413
   }
 
414
 else return x;}
 
415
 
 
416
 
 
417
GEN
 
418
setq_ii(x,all,val)
 
419
     GEN x;
 
420
     GEN val;
 
421
     object *all;
 
422
{int n= replace_copy1(val,x);
 
423
 if (n)
 
424
   { *all = alloc_bignum_static(n/sizeof(int));
 
425
     return replace_copy2(val,MP(*all));
 
426
   }
 
427
 else return x;}
 
428
 
 
429
 
 
430
 
 
431
 
 
432
void
 
433
isetq_fix(var,s)
 
434
     GEN var;
 
435
     int s;
 
436
{/* if (var==0) FEerror("unitialized integer var",0); */
 
437
 STOI(s,var);
 
438
}
 
439
 
 
440
GEN
 
441
icopy_bignum(a,y)
 
442
     object a;
 
443
     GEN y;
 
444
{ int *ucop = (int *)MP(a); 
 
445
  int *vcop = (int *) (y);
 
446
  int j = lgef(ucop);
 
447
  {while(--j >= 0) 
 
448
     { *vcop++ = *ucop++;}
 
449
   setlg(y,a->big.big_length);
 
450
   return y;}}
 
451
     
 
452
 
 
453
GEN
 
454
icopy_fixnum(a,y)
 
455
     object a;
 
456
     GEN y;
 
457
       
 
458
{ int x= fix(a);
 
459
  if(!x) return gzero;
 
460
  y[0]=INT_FLAG+3;
 
461
  if(x>0) {y[1]=0x1000003;y[2]=x;}
 
462
  else{y[1]=0xff000003;y[2]= -x;}
 
463
  return y;
 
464
}
 
465
     
 
466
 
 
467
 
 
468
 
 
469
GEN     gnil,gzero,gun,gdeux,ghalf,gi;
 
470
plong    lontyp[30]={0,0x10000,0x10000,1,1,1,1,2,1,0,2,2,1,1,1,0,1,1,1,1};
 
471
unsigned plong hiremainder,overflow;
 
472
 
 
473
#ifdef STANDALONE
 
474
#define FEerror printf
 
475
#define make_si_sfun(a,b,c) 
 
476
#endif
 
477
 
 
478
#define INITIAL_PARI_STACK 400
 
479
char initial_pari_stack[400];
 
480
 
 
481
our_ulong bot= (our_ulong) initial_pari_stack;
 
482
our_ulong top = (our_ulong)(initial_pari_stack+INITIAL_PARI_STACK);
 
483
/* not initted */
 
484
our_ulong avma= 0;
 
485
 
 
486
 
 
487
void
 
488
err(s)
 
489
     int s;
 
490
{ switch (s) {
 
491
 case errpile:
 
492
  FEerror("Out of bignum stack space, (si::MULTIPLY-BIGNUM-STACK n) to grow",0);
 
493
 case dvmer1:
 
494
 case diver4:
 
495
 case diver2:
 
496
 case diver1:
 
497
  FEerror("Divide by zero",0);
 
498
 case muler1:
 
499
  FEerror("Multiply overflow",0);
 
500
 case moder1:
 
501
  FEerror("Mod by 0",0);
 
502
 default:
 
503
      FEerror("Integer Arithmetic error",0);
 
504
}}
 
505
 
 
506
 
 
507
 
 
508
 
 
509
multiply_bignum_stack(n)
 
510
     int n;
 
511
{ int parisize = n* (top - bot);
 
512
  in_saved_avma = 0;
 
513
  if (n> 1) 
 
514
    { if (bot != (our_ulong)initial_pari_stack) free(bot);
 
515
      set_pari_stack(parisize);
 
516
    }
 
517
  return parisize;
 
518
}
 
519
  
 
520
set_pari_stack(parisize)
 
521
     int parisize;
 
522
{
 
523
  bot=(plong)malloc(parisize);
 
524
  top = avma = bot + parisize;
 
525
}
 
526
 
 
527
/* things to be done every start */
 
528
init_big1()
 
529
{
 
530
 
 
531
}
 
532
 
 
533
init_big()
 
534
{
 
535
  if (avma==0)
 
536
    { 
 
537
        make_si_sfun("MULTIPLY-BIGNUM-STACK",multiply_bignum_stack,
 
538
                     ARGTYPE1(f_fixnum) | RESTYPE(f_fixnum));
 
539
        avma = top;
 
540
      }
 
541
 /* room for the permanent things */
 
542
 
 
543
  gnil = cgeti(2);gnil[1]=2; setpere(gnil,255);
 
544
  gzero = cgeti(2);gzero[1]=2; setpere(gzero, 255);
 
545
  gun = stoi(1); setpere(gun, 255);
 
546
  gdeux = stoi(2); setpere(gdeux, 255);
 
547
  ghalf = cgetg(3,4);ghalf[1]=un;ghalf[2]=deux; setpere(ghalf, 255);
 
548
  gi = cgetg(3,6); gi[1] = zero; gi[2] = un; setpere(gi, 255);
 
549
 
 
550
  /* set_pari_stack(BIGNUM_STACK_SIZE);*/
 
551
 }
 
552
 
 
553
 
 
554
 
 
555
     
 
556
 
 
557
 
 
558
  
 
559
 
 
560
 
 
561
 
 
562
 
 
563
 
 
564
 
 
565