« back to all changes in this revision
Viewing changes to mpz_extras.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Tim Abbott
- Date: 2008-05-30 00:34:12 UTC
- Revision ID: james.westby@ubuntu.com-20080530003412-84wwi86s4iy2gji1
Tags: upstream-1.06
Import upstream version 1.06
- files added:
- QS
- doc
- graphing
- magma-profiles
- pari-profiles
added
removed
mpz_extras.c
1
/*============================================================================
2
3
Copyright (C) 2007, William Hart
4
5
This file is part of FLINT.
6
7
FLINT is free software; you can redistribute it and/or modify
8
it under the terms of the GNU General Public License as published by
9
the Free Software Foundation; either version 2 of the License, or
10
(at your option) any later version.
11
12
FLINT is distributed in the hope that it will be useful,
13
but WITHOUT ANY WARRANTY; without even the implied warranty of
14
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15
GNU General Public License for more details.
16
17
You should have received a copy of the GNU General Public License
18
along with FLINT; if not, write to the Free Software
19
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21
===============================================================================*/
22
#include <stdio.h>
23
#include <stdlib.h>
24
#include <string.h>
25
#include <gmp.h>
26
#include "mpz_extras.h"
27
#include "flint.h"
28
#include "mpn_extras.h"
29
#include "F_mpn_mul-tuning.h"
30
#include "memory-manager.h"
31
#include "longlong_wrapper.h"
32
#include "longlong.h"
33
34
#define DEBUG2 1
35
#define DEBUG 0
36
37
/*
38
Memory manager to allocate a single mpz_t. It returns a pointer to the mpz_t.
39
mpz_t's should be released in the order they were allocated.
40
*/
41
42
#define RESALLOC 100 //allocate this many mpz_t's at once to save on overheads
43
44
mpz_t** reservoir; // Array of pointers to mpz_t's in the reservoir
45
unsigned long rescount=0; //Next available mpz_t in reservoir
46
unsigned long currentalloc=0; //total number of mpz_t's in reservoir
47
48
mpz_t* F_mpz_alloc(void)
49
{
50
static int initialised = 0;
51
static mpz_t** tempres;
52
mpz_t* alloc_d;
53
54
//allocate another block of mpz_t's if none are currently allocated, or the reservoir is depleted
55
if (rescount==currentalloc) // need more limb_memp_t's
56
{
57
if (!initialised)
58
{
59
reservoir = (mpz_t**)malloc(RESALLOC*sizeof(mpz_t*)); //allocate space for the array of pointers
60
reservoir[0] = (mpz_t*)malloc(RESALLOC*sizeof(mpz_t)); //allocate space for the mpz_t's
61
for (unsigned long i=0; i<RESALLOC-1; i++)
62
{
63
reservoir[i+1]=reservoir[i]+1; //initialise the array
64
mpz_init(*reservoir[i]); //initialise the mpz_t's
65
}
66
mpz_init(*reservoir[RESALLOC-1]);
67
rescount=0;
68
initialised = 1;
69
currentalloc = RESALLOC;
70
} else
71
{
72
//copy old reservoir into larger one
73
tempres = reservoir;
74
reservoir = (mpz_t**)malloc((currentalloc+RESALLOC)*sizeof(mpz_t*));
75
reservoir[currentalloc] = (mpz_t*)malloc(RESALLOC*sizeof(mpz_t));
76
memcpy(reservoir,tempres,currentalloc*sizeof(mpz_t*));
77
for (unsigned long i=currentalloc; i<RESALLOC+currentalloc-1; i++)
78
{
79
reservoir[i+1]=reservoir[i]+1; //initialise the array
80
mpz_init(*reservoir[i]); //initialise the mpz_t's
81
}
82
mpz_init(*reservoir[currentalloc+RESALLOC-1]);
83
84
currentalloc+=RESALLOC;
85
//free old reservoir
86
free(tempres);
87
}
88
}
89
90
alloc_d = reservoir[rescount];
91
rescount++;
92
return alloc_d;
93
}
94
95
/*
96
Release a mpz_t back into the reservoir
97
*/
98
99
void F_mpz_release(void)
100
{
101
rescount--;
102
}
103
104
/*
105
sets res to a*b modulo p
106
assumes res is not p
107
Does not assume a and b are reduced mod p
108
*/
109
110
void F_mpz_mulmod(mpz_t res, mpz_t a, mpz_t b, mpz_t p)
111
{
112
mpz_t* temp = F_mpz_alloc();
113
114
mpz_fdiv_r(*temp,a,p);
115
mpz_fdiv_r(res,b,p);
116
mpz_mul(res,*temp,res);
117
mpz_fdiv_r(res,res,p);
118
119
F_mpz_release();
120
121
return;
122
}
123
124
/*
125
Sets res to a*b modulo p
126
Does not assume a and b are reduced mod p
127
*/
128
129
unsigned long F_mpz_mulmod_ui(mpz_t res, mpz_t a, mpz_t b, unsigned long p)
130
{
131
unsigned long p1, p2, al, bl;
132
133
al = mpz_fdiv_r_ui(res, a, p);
134
bl = mpz_fdiv_r_ui(res, b, p);
135
136
umul_ppmm(p2, p1, al, bl);
137
138
unsigned long norm, q, r;
139
140
if (p2 >= p) p2 %= p;
141
142
#if UDIV_NEEDS_NORMALIZATION
143
count_lead_zeros(norm, p);
144
udiv_qrnnd(q, r, (p2<<norm) + (p1>>(FLINT_BITS-norm)), p1<<norm, p<<norm);
145
#else
146
udiv_qrnnd(q, r, p2, p1, p);
147
#endif
148
149
mpz_set_ui(res, r);
150
151
return r;
152
}
153
154
/*
155
Sets res to the square root of a modulo p for a prime p
156
Returns 0 if a is not a square modulo p
157
*/
158
159
int F_mpz_sqrtmod(mpz_t res, mpz_t a, mpz_t p)
160
{
161
int r,k,m;
162
mpz_t* p1 = F_mpz_alloc();
163
mpz_t* two = F_mpz_alloc();
164
mpz_t* b = F_mpz_alloc();
165
mpz_t* g = F_mpz_alloc();
166
mpz_t* bpow = F_mpz_alloc();
167
mpz_t* gpow = F_mpz_alloc();
168
mpz_t* mk = F_mpz_alloc();
169
170
if (mpz_kronecker(a,p)!=1)
171
{
172
mpz_set_ui(res,0);
173
return 0; //return 0 if a is not a square mod p
174
}
175
176
if ((mpz_cmp_ui(a,0)==0)||(mpz_cmp_ui(a,1)==0))
177
{
178
mpz_set(res,a);
179
return 1;
180
}
181
182
if ((mpz_tstbit(p,0)==1)||(mpz_tstbit(p,1)==1))
183
{
184
mpz_add_ui(*p1,p,1);
185
mpz_fdiv_q_2exp(*p1,*p1,2);
186
mpz_powm(res,a,*p1,p);
187
return 1;
188
}
189
190
mpz_set_ui(*two,2);
191
192
mpz_sub_ui(*p1,p,1);
193
r = mpz_remove(*p1,*p1,*two);
194
mpz_powm(*b,a,*p1,p);
195
for (k=2; ;k++)
196
{
197
if (mpz_ui_kronecker(k,p) == -1) break;
198
}
199
mpz_set_ui(*mk,k);
200
mpz_powm(*g,*mk,*p1,p);
201
mpz_add_ui(*p1,*p1,1);
202
mpz_divexact_ui(*p1,*p1,2);
203
mpz_powm(res,a,*p1,p);
204
if (!mpz_cmp_ui(*b,1)) return 1;
205
206
while (mpz_cmp_ui(*b,1))
207
{
208
mpz_set(*bpow,*b);
209
for (m=1; (m<=r-1) && (mpz_cmp_ui(*bpow,1));m++)
210
{
211
mpz_powm_ui(*bpow,*bpow,2,p);
212
}
213
mpz_set(*gpow,*g);
214
for (int i = 1;i<= r-m-1;i++)
215
{
216
mpz_powm_ui(*gpow,*gpow,2,p);
217
};
218
F_mpz_mulmod(res,res,*gpow,p);
219
mpz_powm_ui(*gpow,*gpow,2,p);
220
F_mpz_mulmod(*b,*b,*gpow,p);
221
mpz_set(*gpow,*g);
222
r = m;
223
}
224
225
F_mpz_release();F_mpz_release();F_mpz_release();F_mpz_release();
226
F_mpz_release();F_mpz_release();F_mpz_release();
227
228
return 1;
229
}
230
231
/*
232
Computes the square root of a modulo p^k when given z, the square root mod p^(k-1)
233
*/
234
235
void __sqrtmodpow(mpz_t res, mpz_t z, mpz_t a, mpz_t pk, mpz_t tempsqpow, mpz_t inv)
236
{
237
mpz_mul_ui(inv,z,2);
238
mpz_invert(inv,inv,pk);
239
mpz_set(tempsqpow,a);
240
mpz_submul(tempsqpow,z,z);
241
mpz_fdiv_r(tempsqpow,tempsqpow,pk);
242
F_mpz_mulmod(tempsqpow,tempsqpow,inv,pk);
243
mpz_add(tempsqpow,tempsqpow,z);
244
mpz_set(res,tempsqpow);
245
246
return;
247
}
248
249
/*
250
Computes the square root of a modulo p^k when given z, the square root mod p^{k-1}
251
*/
252
253
void F_mpz_sqrtmodpklift(mpz_t res, mpz_t z, mpz_t a, mpz_t pk)
254
{
255
mpz_t* tempsqpow = F_mpz_alloc();
256
mpz_t* inv = F_mpz_alloc();
257
258
__sqrtmodpow(res, z, a, pk, *tempsqpow, *inv);
259
260
F_mpz_release();F_mpz_release();
261
}
262
263
/*
264
computes the square root of a modulo p^k when given the square root modulo p
265
*/
266
267
void F_mpz_sqrtmodptopk(mpz_t res, mpz_t sqrt, mpz_t a, mpz_t p, int k)
268
{
269
mpz_t* tempsqpow = F_mpz_alloc();
270
mpz_t* inv = F_mpz_alloc();
271
mpz_t* pk = F_mpz_alloc();
272
273
mpz_set(res,sqrt);
274
mpz_set(*pk,p);
275
for (int i = 2; i<=k; i++)
276
{
277
mpz_mul(*pk,*pk,p);
278
__sqrtmodpow(res,res,a,*pk, *tempsqpow, *inv);
279
}
280
281
F_mpz_release();F_mpz_release();F_mpz_release();
282
}
283
284
/*
285
Computes the square root of a modulo p^k
286
Returns 0 if the square root of a does not exist mod p
287
*/
288
289
int F_mpz_sqrtmodpk(mpz_t res, mpz_t a, mpz_t p, int k)
290
{
291
int exists;
292
mpz_t* sqrtmodp = F_mpz_alloc();
293
294
exists = F_mpz_sqrtmod(*sqrtmodp,a,p);
295
F_mpz_sqrtmodptopk(res,*sqrtmodp,a,p,k);
296
297
F_mpz_release();
298
299
return exists;
300
}
301
302
303
/*
304
Find res mod n=n1*n2 such that res = x1 mod n1 and res = x2 mod n2
305
*/
306
307
void F_mpz_CRT(mpz_t res, mpz_t n, mpz_t x1, mpz_t x2, mpz_t n1, mpz_t n2)
308
{
309
mpz_t* ntemp = F_mpz_alloc();
310
mpz_t* restemp = F_mpz_alloc();
311
mpz_t* chtemp = F_mpz_alloc();
312
313
mpz_mul(*ntemp,n1,n2);
314
mpz_invert(*restemp,n2,n1);
315
F_mpz_mulmod(*restemp,res,n2,*ntemp);
316
F_mpz_mulmod(*restemp,*restemp,x1,*ntemp);
317
318
mpz_invert(*chtemp,n1,n2);
319
F_mpz_mulmod(*chtemp,*chtemp,n1,*ntemp);
320
F_mpz_mulmod(*chtemp,*chtemp,x2,*ntemp);
321
322
mpz_add(res,*restemp,*chtemp);
323
mpz_fdiv_r(res,res,*ntemp);
324
mpz_set(n,*ntemp);
325
326
F_mpz_release();F_mpz_release();F_mpz_release();
327
328
return;
329
}
330
331
/*
332
Compute the Montgomery reduced form of a mod m
333
Returns n such that m < 2^n (n will be divisible by FLINT_BITS)
334
Assumes a is already reduced mod m
335
*/
336
337
unsigned long F_mpz_mont_red(mpz_t res, mpz_t a, mpz_t m)
338
{
339
unsigned long n = mpz_size(m)*FLINT_BITS;
340
341
mpz_mul_2exp(res, a, n);
342
mpz_mod(res, res, m);
343
344
return n;
345
}
346
347
/*
348
Compute the Montgomery multiplication r = a*b mod m assuming a and b are in
349
Montgomery form with respect to 2^n where m < 2^n and R is -m mod 2^n
350
*/
351
352
void F_mpz_mont_mul(mpz_t res, mpz_t a, mpz_t b, mpz_t m, mpz_t R, unsigned long n)
353
{
354
mpz_t x, s;
355
mpz_init(x);
356
mpz_init(s);
357
358
mpz_mul(x, a, b);
359
mpz_fdiv_r_2exp(s, x, n);
360
mpz_mul(s, s, R);
361
mpz_fdiv_r_2exp(s, s, n);
362
mpz_mul(res, s, m);
363
mpz_add(res, res, x);
364
mpz_fdiv_q_2exp(res, res, n);
365
366
if (mpz_cmp(res, m) >= 0) mpz_sub(res, res, m);
367
368
mpz_clear(x);
369
mpz_clear(s);
370
}
371
372
/*
373
Compute a^exp mod m using Montgomery reduction
374
Requires that m is odd and positive and that exp is positive
375
*/
376
377
void F_mpz_expmod_mont(mpz_t res, mpz_t a, mpz_t exp, mpz_t m)
378
{
379
unsigned long n;
380
unsigned long bits = mpz_sizeinbase(exp, 2);
381
mpz_t aRED;
382
mpz_t powRED;
383
mpz_t R;
384
mpz_t temp;
385
int flag = 0;
386
387
mpz_init(aRED);
388
mpz_init(powRED);
389
mpz_init(R);
390
mpz_init(temp);
391
392
n = F_mpz_mont_red(aRED, a, m);
393
394
mpz_set_ui(temp, 1);
395
mpz_mul_2exp(temp, temp, n);
396
mpz_invert(R, m, temp);
397
mpz_sub(R, temp, R);
398
if (mpz_cmp(R, temp) == 0) mpz_sub(R, R, temp);
399
400
mpz_set(powRED, aRED);
401
#ifdef DEBUG
402
gmp_printf("powRED = %Zd\n", powRED);
403
#endif
404
405
for (unsigned long i = 0; i < bits - 1; i++)
406
{
407
if (mpz_tstbit(exp, i))
408
{
409
if (flag) F_mpz_mont_mul(res, res, powRED, m, R, n);
410
else
411
{
412
mpz_set(res, powRED);
413
flag = 1;
414
}
415
}
416
F_mpz_mont_mul(powRED, powRED, powRED, m, R, n);
417
#ifdef DEBUG
418
gmp_printf("powRED = %Zd\n", powRED);
419
#endif
420
}
421
422
if (flag) F_mpz_mont_mul(res, res, powRED, m, R, n);
423
else mpz_set(res, powRED);
424
425
mpz_set_ui(temp, 1);
426
F_mpz_mont_mul(res, res, temp, m, R, n);
427
428
mpz_clear(temp);
429
mpz_clear(R);
430
mpz_clear(powRED);
431
mpz_clear(aRED);
432
}
433
434
void F_mpz_divrem_BZ(mpz_t Q, mpz_t R, mpz_t A, mpz_t B)
435
{
436
unsigned long n = mpz_size(B);
437
unsigned long m = mpz_size(A) - n;
438
439
if ((long) m < 0)
440
{
441
mpz_set_ui(Q, 0);
442
mpz_set(R, A);
443
return;
444
}
445
446
if (m < 64)
447
{
448
mpz_fdiv_qr(Q, R, A, B);
449
return;
450
}
451
452
unsigned long k = m/2;
453
454
mpz_t * B0 = F_mpz_alloc();
455
mpz_t * B1 = F_mpz_alloc();
456
mpz_t * A0 = F_mpz_alloc();
457
mpz_t * A1 = F_mpz_alloc();
458
mpz_t * Q0 = F_mpz_alloc();
459
mpz_t * Q1 = F_mpz_alloc();
460
mpz_t * R0 = F_mpz_alloc();
461
mpz_t * R1 = F_mpz_alloc();
462
mpz_t * temp = F_mpz_alloc();
463
mpz_t * temp2 = F_mpz_alloc();
464
mpz_t * temp3 = F_mpz_alloc();
465
466
467
mpz_fdiv_q_2exp(*B1, B, FLINT_BITS*k);
468
mpz_fdiv_q_2exp(*A1, A, FLINT_BITS*2*k);
469
470
F_mpz_divrem_BZ(*Q1, *R1, *A1, *B1);
471
mpz_fdiv_r_2exp(*B0, B, FLINT_BITS*k);
472
mpz_fdiv_r_2exp(*A0, A, FLINT_BITS*2*k);
473
mpz_mul_2exp(*temp, *R1, FLINT_BITS*2*k);
474
mpz_add(*temp, *temp, *A0);
475
mpz_mul_2exp(*temp2, *Q1, FLINT_BITS*k);
476
mpz_mul(*temp2, *temp2, *B0);
477
mpz_sub(*temp, *temp, *temp2);
478
mpz_mul_2exp(*temp2, B, FLINT_BITS*k);
479
480
while (mpz_cmp_ui(*temp, 0) < 0)
481
{
482
mpz_sub_ui(*Q1, *Q1, 1);
483
mpz_add(*temp, *temp, *temp2);
484
}
485
mpz_fdiv_q_2exp(*temp2, *temp, FLINT_BITS*k);
486
F_mpz_divrem_BZ(*Q0, *R0, *temp2, *B1);
487
488
mpz_fdiv_r_2exp(*temp2, *temp, FLINT_BITS*k);
489
mpz_mul_2exp(R, *R0, FLINT_BITS*k);
490
mpz_add(R, R, *temp2);
491
mpz_submul(R, *Q0, *B0);
492
while (mpz_cmp_ui(R, 0) < 0)
493
{
494
mpz_sub_ui(*Q0, *Q0, 1);
495
mpz_add(R, R, B);
496
}
497
mpz_mul_2exp(Q, *Q1, FLINT_BITS*k);
498
mpz_add(Q, Q, *Q0);
499
500
F_mpz_release(); F_mpz_release(); F_mpz_release(); F_mpz_release();
501
F_mpz_release(); F_mpz_release(); F_mpz_release(); F_mpz_release();
502
F_mpz_release(); F_mpz_release(); F_mpz_release();
503
}
504
505
void F_mpz_rem_BZ(mpz_t R, mpz_t A, mpz_t B)
506
{
507
unsigned long n = mpz_size(B);
508
unsigned long m = mpz_size(A) - n;
509
510
if ((long) m < 0)
511
{
512
mpz_set(R, A);
513
return;
514
}
515
516
if (m < 64)
517
{
518
mpz_fdiv_r(R, A, B);
519
return;
520
}
521
522
unsigned long k = m/2;
523
524
mpz_t * B0 = F_mpz_alloc();
525
mpz_t * B1 = F_mpz_alloc();
526
mpz_t * A0 = F_mpz_alloc();
527
mpz_t * A1 = F_mpz_alloc();
528
mpz_t * Q0 = F_mpz_alloc();
529
mpz_t * Q1 = F_mpz_alloc();
530
mpz_t * R0 = F_mpz_alloc();
531
mpz_t * R1 = F_mpz_alloc();
532
mpz_t * temp = F_mpz_alloc();
533
mpz_t * temp2 = F_mpz_alloc();
534
mpz_t * temp3 = F_mpz_alloc();
535
536
537
mpz_fdiv_q_2exp(*B1, B, FLINT_BITS*k);
538
mpz_fdiv_q_2exp(*A1, A, FLINT_BITS*2*k);
539
540
F_mpz_divrem_BZ(*Q1, *R1, *A1, *B1);
541
mpz_fdiv_r_2exp(*B0, B, FLINT_BITS*k);
542
mpz_fdiv_r_2exp(*A0, A, FLINT_BITS*2*k);
543
mpz_mul_2exp(*temp, *R1, FLINT_BITS*2*k);
544
mpz_add(*temp, *temp, *A0);
545
mpz_mul_2exp(*temp2, *Q1, FLINT_BITS*k);
546
mpz_mul(*temp2, *temp2, *B0);
547
mpz_sub(*temp, *temp, *temp2);
548
mpz_mul_2exp(*temp2, B, FLINT_BITS*k);
549
550
while (mpz_cmp_ui(*temp, 0) < 0)
551
{
552
mpz_sub_ui(*Q1, *Q1, 1);
553
mpz_add(*temp, *temp, *temp2);
554
}
555
mpz_fdiv_q_2exp(*temp2, *temp, FLINT_BITS*k);
556
F_mpz_divrem_BZ(*Q0, *R0, *temp2, *B1);
557
558
mpz_fdiv_r_2exp(*temp2, *temp, FLINT_BITS*k);
559
mpz_mul_2exp(R, *R0, FLINT_BITS*k);
560
mpz_add(R, R, *temp2);
561
mpz_submul(R, *Q0, *B0);
562
while (mpz_cmp_ui(R, 0) < 0)
563
{
564
mpz_add(R, R, B);
565
}
566
567
F_mpz_release(); F_mpz_release(); F_mpz_release(); F_mpz_release();
568
F_mpz_release(); F_mpz_release(); F_mpz_release(); F_mpz_release();
569
F_mpz_release(); F_mpz_release(); F_mpz_release();
570
}
571
572
void F_mpz_mulmod_BZ(mpz_t res, mpz_t a, mpz_t b, mpz_t m)
573
{
574
mpz_t * temp = F_mpz_alloc();
575
576
mpz_mul(*temp, a, b);
577
F_mpz_rem_BZ(res, *temp, m);
578
579
F_mpz_release();
580
}
581
582
583
void F_mpz_expmod_BZ(mpz_t res, mpz_t a, mpz_t exp, mpz_t m)
584
{
585
unsigned long n;
586
unsigned long bits = mpz_sizeinbase(exp, 2);
587
mpz_t aRED;
588
mpz_t powRED;
589
mpz_t temp;
590
int flag = 0;
591
592
mpz_init(aRED);
593
mpz_init(powRED);
594
mpz_init(temp);
595
596
mpz_set(powRED, a);
597
#if DEBUG
598
gmp_printf("powRED = %Zd\n", powRED);
599
#endif
600
601
for (unsigned long i = 0; i < bits - 1; i++)
602
{
603
if (mpz_tstbit(exp, i))
604
{
605
if (flag) F_mpz_mulmod_BZ(res, res, powRED, m);
606
else
607
{
608
mpz_set(res, powRED);
609
flag = 1;
610
}
611
}
612
F_mpz_mulmod_BZ(powRED, powRED, powRED, m);
613
#if DEBUG
614
gmp_printf("powRED = %Zd\n", powRED);
615
#endif
616
}
617
618
if (flag) F_mpz_mulmod_BZ(res, res, powRED, m);
619
else mpz_set(res, powRED);
620
621
mpz_clear(temp);
622
mpz_clear(powRED);
623
mpz_clear(aRED);
624
}
625
626
/*
627
Large integer multiplication code
628
*/
629
630
void __F_mpz_mul(mpz_t res, mpz_t a, mpz_t b, unsigned long twk)
631
{
632
unsigned long sa = mpz_size(a);
633
unsigned long sb = mpz_size(b);
634
635
if (sa+sb > FLINT_FFT_LIMBS_CROSSOVER)
636
{
637
unsigned long s1 = (FLINT_BIT_COUNT(a->_mp_d[sa-1]) + FLINT_BIT_COUNT(b->_mp_d[sb-1]) <= FLINT_BITS);
638
639
mp_limb_t* output =
640
(mp_limb_t*) flint_stack_alloc(sa + sb - s1);
641
__F_mpn_mul(output, a->_mp_d, sa, b->_mp_d, sb, twk);
642
mpz_import(res, sa+sb-s1, -1, sizeof(mp_limb_t), 0, 0, output);
643
if (mpz_sgn(res) != mpz_sgn(a)*mpz_sgn(b)) mpz_neg(res,res);
644
flint_stack_release();
645
} else mpz_mul(res, a, b);
646
}
647
648
void F_mpz_mul(mpz_t res, mpz_t a, mpz_t b)
649
{
650
unsigned long sa = mpz_size(a);
651
unsigned long sb = mpz_size(b);
652
653
if (sa+sb > FLINT_FFT_LIMBS_CROSSOVER)
654
{
655
unsigned long s1 = (FLINT_BIT_COUNT(a->_mp_d[sa-1]) + FLINT_BIT_COUNT(b->_mp_d[sb-1]) <= FLINT_BITS);
656
mp_limb_t* output =
657
(mp_limb_t*) flint_stack_alloc(sa + sb - s1);
658
F_mpn_mul(output, a->_mp_d, sa, b->_mp_d, sb);
659
mpz_import(res, sa+sb-s1, -1, sizeof(mp_limb_t), 0, 0, output);
660
if (mpz_sgn(res) != mpz_sgn(a)*mpz_sgn(b)) mpz_neg(res,res);
661
flint_stack_release();
662
} else mpz_mul(res, a, b);
663
}
664
Loggerhead is a web-based interface for Breezy
Version: 2.0.1