1
/* Copyright (C) 2000  The PARI group.
2

3
This file is part of the PARI/GP package.
4

5
PARI/GP is free software; you can redistribute it and/or modify it under the
6
terms of the GNU General Public License as published by the Free Software
7
Foundation; either version 2 of the License, or (at your option) any later
8
version. It is distributed in the hope that it will be useful, but WITHOUT
9
ANY WARRANTY WHATSOEVER.
10

11
Check the License for details. You should have received a copy of it, along
12
with the package; see the file 'COPYING'. If not, write to the Free Software
13
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14

15
/*********************************************************************/
16
/**                                                                 **/
17
/**                     ARITHMETIC FUNCTIONS                        **/
18
/**                         (first part)                            **/
19
/**                                                                 **/
20
/*********************************************************************/
21
#include "pari.h"
22
#include "paripriv.h"
23

24
#define DEBUGLEVEL DEBUGLEVEL_arith
25

26
/******************************************************************/
27
/*                                                                */
28
/*                 GENERATOR of (Z/mZ)*                           */
29
/*                                                                */
30
/******************************************************************/
31
static GEN
32
remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
33
static ulong
34
u_remove2(ulong q) { return q >> vals(q); }
35
GEN
36
odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
37
static GEN
38
u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
39
/* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
40
 * (all prime divisors of q); return the q/l, l in L0 */
41
static GEN
42
is_gener_expo(GEN p, GEN L0)
43
{
44
  GEN L, q = shifti(p,-1);
45
  long i, l;
46
  if (L0) {
47
    l = lg(L0);
48
    L = cgetg(l, t_VEC);
49
  } else {
50
    L0 = L = odd_prime_divisors(q);
51
    l = lg(L);
52
  }
53
  for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
54
  return L;
55
}
56
static GEN
57
u_is_gener_expo(ulong p, GEN L0)
58
{
59
  const ulong q = p >> 1;
60
  long i;
61
  GEN L;
62
  if (!L0) L0 = u_odd_prime_divisors(q);
63
  L = cgetg_copy(L0,&i);
64
  while (--i) L[i] = q / uel(L0,i);
65
  return L;
66
}
67

68
int
69
is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
70
{
71
  long i;
72
  if (krouu(x, p) >= 0) return 0;
73
  for (i=lg(L)-1; i; i--)
74
  {
75
    ulong t = Fl_powu(x, uel(L,i), p);
76
    if (t == p_1 || t == 1) return 0;
77
  }
78
  return 1;
79
}
80
/* assume p prime */
81
ulong
82
pgener_Fl_local(ulong p, GEN L0)
83
{
84
  const pari_sp av = avma;
85
  const ulong p_1 = p-1;
86
  long x;
87
  GEN L;
88
  if (p <= 19) switch(p)
89
  { /* quick trivial cases */
90
    case 2:  return 1;
91
    case 7:
92
    case 17: return 3;
93
    default: return 2;
94
  }
95
  L = u_is_gener_expo(p,L0);
96
  for (x = 2;; x++)
97
    if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
98
}
99
ulong
100
pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
101

102
/* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
103
 * but wasteful) */
104
int
105
is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
106
{
107
  long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
108
  if (t >= 0) return 0;
109
  for (i = lg(L)-1; i; i--)
110
  {
111
    GEN t = Fp_pow(x, gel(L,i), p);
112
    if (equalii(t, p_1) || equali1(t)) return 0;
113
  }
114
  return 1;
115
}
116

117
/* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
118
GEN
119
pgener_Fp_local(GEN p, GEN L0)
120
{
121
  pari_sp av0 = avma;
122
  GEN x, p_1, L;
123
  if (lgefint(p) == 3)
124
  {
125
    ulong z;
126
    if (p[2] == 2) return gen_1;
127
    if (L0) L0 = ZV_to_nv(L0);
128
    z = pgener_Fl_local(uel(p,2), L0);
129
    set_avma(av0); return utoipos(z);
130
  }
131
  p_1 = subiu(p,1); L = is_gener_expo(p, L0);
132
  x = utoipos(2);
133
  for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
134
  set_avma(av0); return utoipos(uel(x,2));
135
}
136

137
GEN
138
pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
139

140
ulong
141
pgener_Zl(ulong p)
142
{
143
  if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
144
  /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
145
  if (p == 40487) return 10;
146
#ifndef LONG_IS_64BIT
147
  return pgener_Fl(p);
148
#else
149
  if (p < (1UL<<32)) return pgener_Fl(p);
150
  else
151
  {
152
    const pari_sp av = avma;
153
    const ulong p_1 = p-1;
154
    long x ;
155
    GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
156
    for (x=2;;x++)
157
      if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2)))
158
        return gc_ulong(av, x);
159
  }
160
#endif
161
}
162

163
/* p prime. Return a primitive root modulo p^e, e > 1 */
164
GEN
165
pgener_Zp(GEN p)
166
{
167
  if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
168
  else
169
  {
170
    const pari_sp av = avma;
171
    GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
172
    GEN x = utoipos(2);
173
    for (;; x[2]++)
174
      if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
175
    set_avma(av); return utoipos(uel(x,2));
176
  }
177
}
178

179
static GEN
180
gener_Zp(GEN q, GEN F)
181
{
182
  GEN p = NULL;
183
  long e = 0;
184
  if (F)
185
  {
186
    GEN P = gel(F,1), E = gel(F,2);
187
    long i, l = lg(P);
188
    for (i = 1; i < l; i++)
189
    {
190
      p = gel(P,i);
191
      if (absequaliu(p, 2)) continue;
192
      if (i < l-1) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
193
      e = itos(gel(E,i));
194
    }
195
    if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
196
  }
197
  else
198
    e = Z_isanypower(q, &p);
199
  return e > 1? pgener_Zp(p): pgener_Fp(q);
200
}
201

202
GEN
203
znprimroot(GEN N)
204
{
205
  pari_sp av = avma;
206
  GEN x, n, F;
207

208
  if ((F = check_arith_non0(N,"znprimroot")))
209
  {
210
    F = clean_Z_factor(F);
211
    N = typ(N) == t_VEC? gel(N,1): factorback(F);
212
  }
213
  N = absi_shallow(N);
214
  if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
215
  switch(mod4(N))
216
  {
217
    case 0: /* N = 0 mod 4 */
218
      pari_err_DOMAIN("znprimroot", "argument","=",N,N);
219
      x = NULL; break;
220
    case 2: /* N = 2 mod 4 */
221
      n = shifti(N,-1); /* becomes odd */
222
      x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
223
      break;
224
    default: /* N odd */
225
      x = gener_Zp(N,F);
226
      break;
227
  }
228
  return gerepilecopy(av, mkintmod(x, N));
229
}
230

231
/* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
232
GEN
233
rootsof1_Fp(GEN n, GEN p)
234
{
235
  pari_sp av = avma;
236
  GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
237
  GEN z = pgener_Fp_local(p, L);
238
  z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
239
  return gerepileuptoint(av, z);
240
}
241

242
GEN
243
rootsof1u_Fp(ulong n, GEN p)
244
{
245
  pari_sp av = avma;
246
  GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
247
  z = pgener_Fp_local(p, Flv_to_ZV(L));
248
  z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
249
  return gerepileuptoint(av, z);
250
}
251

252
ulong
253
rootsof1_Fl(ulong n, ulong p)
254
{
255
  pari_sp av = avma;
256
  GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
257
  ulong z = pgener_Fl_local(p, L);
258
  z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
259
  return gc_ulong(av,z);
260
}
261

262
/*********************************************************************/
263
/**                                                                 **/
264
/**                     INVERSE TOTIENT FUNCTION                    **/
265
/**                                                                 **/
266
/*********************************************************************/
267
/* N t_INT, L a ZV containing all prime divisors of N, and possibly other
268
 * primes. Return factor(N) */
269
GEN
270
Z_factor_listP(GEN N, GEN L)
271
{
272
  long i, k, l = lg(L);
273
  GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
274
  for (i = k = 1; i < l; i++)
275
  {
276
    GEN p = gel(L,i);
277
    long v = Z_pvalrem(N, p, &N);
278
    if (v)
279
    {
280
      gel(P,k) = p;
281
      gel(E,k) = utoipos(v);
282
      k++;
283
    }
284
  }
285
  setlg(P, k);
286
  setlg(E, k); return mkmat2(P,E);
287
}
288

289
/* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
290
 * L is a list of primes containing all prime divisors of n. */
291
static long
292
istotient_i(GEN n, GEN m, GEN L, GEN *px)
293
{
294
  pari_sp av = avma, av2;
295
  GEN k, D;
296
  long i, v;
297
  if (m && mod2(n))
298
  {
299
    if (!equali1(n)) return 0;
300
    if (px) *px = gen_1;
301
    return 1;
302
  }
303
  D = divisors(Z_factor_listP(shifti(n, -1), L));
304
  /* loop through primes p > m, d = p-1 | n */
305
  av2 = avma;
306
  if (!m)
307
  { /* special case p = 2, d = 1 */
308
    k = n;
309
    for (v = 1;; v++) {
310
      if (istotient_i(k, gen_2, L, px)) {
311
        if (px) *px = shifti(*px, v);
312
        return 1;
313
      }
314
      if (mod2(k)) break;
315
      k = shifti(k,-1);
316
    }
317
    set_avma(av2);
318
  }
319
  for (i = 1; i < lg(D); ++i)
320
  {
321
    GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
322
    if (m && cmpii(d, m) < 0) continue;
323
    p = addiu(d, 1);
324
    if (!isprime(p)) continue;
325
    k = diviiexact(n, d);
326
    for (v = 1;; v++) {
327
      GEN r;
328
      if (istotient_i(k, p, L, px)) {
329
        if (px) *px = mulii(*px, powiu(p, v));
330
        return 1;
331
      }
332
      k = dvmdii(k, p, &r);
333
      if (r != gen_0) break;
334
    }
335
    set_avma(av2);
336
  }
337
  return gc_long(av,0);
338
}
339

340
/* find x such that phi(x) = n */
341
long
342
istotient(GEN n, GEN *px)
343
{
344
  pari_sp av = avma;
345
  if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
346
  if (signe(n) < 1) return 0;
347
  if (mod2(n))
348
  {
349
    if (!equali1(n)) return 0;
350
    if (px) *px = gen_1;
351
    return 1;
352
  }
353
  if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
354
  {
355
    if (!px) set_avma(av);
356
    else
357
      *px = gerepileuptoint(av, *px);
358
    return 1;
359
  }
360
  return gc_long(av,0);
361
}
362

363
/*********************************************************************/
364
/**                                                                 **/
365
/**                     INTEGRAL LOGARITHM                          **/
366
/**                                                                 **/
367
/*********************************************************************/
368

369
/* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
370
 * e = floor(log_y B). Set *ptq = y^e if non-NULL */
371
long
372
ulogintall(ulong B, ulong y, ulong *ptq)
373
{
374
  ulong r, r2;
375
  long e;
376

377
  if (y == 2)
378
  {
379
    long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
380
    if (ptq) *ptq = 1UL << eB;
381
    return eB;
382
  }
383
  r = y, r2 = 1UL;
384
  for (e=1;; e++)
385
  { /* here, r = y^e, r2 = y^(e-1) */
386
    if (r >= B)
387
    {
388
      if (r != B) { e--; r = r2; }
389
      if (ptq) *ptq = r;
390
      return e;
391
    }
392
    r2 = r;
393
    r = umuluu_or_0(y, r);
394
    if (!r)
395
    {
396
      if (ptq) *ptq = r2;
397
      return e;
398
    }
399
  }
400
}
401

402
/* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
403
 * e = floor(log_y B). Set *ptq = y^e if non-NULL */
404
long
405
logintall(GEN B, GEN y, GEN *ptq)
406
{
407
  pari_sp av;
408
  long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
409
  GEN q, pow2;
410

411
  if (lgefint(B) == 3)
412
  {
413
    ulong q;
414
    if (lgefint(y) > 3)
415
    {
416
      if (ptq) *ptq = gen_1;
417
      return 0;
418
    }
419
    if (!ptq) return ulogintall(B[2], y[2], NULL);
420
    e = ulogintall(B[2], y[2], &q);
421
    *ptq = utoi(q); return e;
422
  }
423
  if (equaliu(y,2))
424
  {
425
    if (ptq) *ptq = int2n(eB);
426
    return eB;
427
  }
428
  av = avma;
429
  ey = expi(y);
430
  /* eB/(ey+1) - 1 < e <= eB/ey */
431
  emax = eB/ey;
432
  if (emax <= 13) /* e small, be naive */
433
  {
434
    GEN r = y, r2 = gen_1;
435
    for (e=1;; e++)
436
    { /* here, r = y^e, r2 = y^(e-1) */
437
      long fl = cmpii(r, B);
438
      if (fl >= 0)
439
      {
440
        if (fl) { e--; cgiv(r); r = r2; }
441
        if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
442
        return e;
443
      }
444
      r2 = r; r = mulii(r,y);
445
    }
446
  }
447
  /* e >= 13 ey / (ey+1) >= 6.5 */
448

449
  /* binary splitting: compute bits of e one by one */
450
  /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
451
  pow2 = new_chunk((long)log2(eB)+2);
452
  gel(pow2,0) = y;
453
  for (i=0, q=y;; )
454
  {
455
    GEN r = gel(pow2,i); /* r = y^2^i */
456
    long fl = cmpii(r,B);
457
    if (!fl)
458
    {
459
      e = 1L<<i;
460
      if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
461
      return e;
462
    }
463
    if (fl > 0) { i--; break; }
464
    q = r;
465
    if (1L<<(i+1) > emax) break;
466
    gel(pow2,++i) = sqri(q);
467
  }
468

469
  for (e = 1L<<i;;)
470
  { /* y^e = q < B < r = q * y^(2^i) */
471
    pari_sp av2 = avma;
472
    long fl;
473
    GEN r;
474
    if (--i < 0) break;
475
    r = mulii(q, gel(pow2,i));
476
    fl = cmpii(r, B);
477
    if (fl > 0) set_avma(av2);
478
    else
479
    {
480
      e += (1L<<i);
481
      q = r;
482
      if (!fl) break; /* B = r */
483
    }
484
  }
485
  if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
486
  return e;
487
}
488

489
long
490
logint0(GEN B, GEN y, GEN *ptq)
491
{
492
  if (typ(B) != t_INT) pari_err_TYPE("logint",B);
493
  if (signe(B) <= 0) pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
494
  if (typ(y) != t_INT) pari_err_TYPE("logint",y);
495
  if (cmpis(y, 2) < 0) pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
496
  return logintall(B,y,ptq);
497
}
498

499
/*********************************************************************/
500
/**                                                                 **/
501
/**                     INTEGRAL SQUARE ROOT                        **/
502
/**                                                                 **/
503
/*********************************************************************/
504
GEN
505
sqrtint(GEN a)
506
{
507
  if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
508
  switch (signe(a))
509
  {
510
    case 1: return sqrti(a);
511
    case 0: return gen_0;
512
    default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
513
  }
514
  return NULL; /* LCOV_EXCL_LINE */
515
}
516
GEN
517
sqrtint0(GEN a, GEN *r)
518
{
519
  if (!r) return sqrtint(a);
520
  if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
521
  switch (signe(a))
522
  {
523
    case 1: return sqrtremi(a, r);
524
    case 0: *r = gen_0; return gen_0;
525
    default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
526
  }
527
  return NULL; /* LCOV_EXCL_LINE */
528
}
529

530
/*********************************************************************/
531
/**                                                                 **/
532
/**                      PERFECT SQUARE                             **/
533
/**                                                                 **/
534
/*********************************************************************/
535
static int
536
carremod(ulong A)
537
{
538
  const int carresmod64[]={
539
    1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,1,0,0,0,0,
540
    0,0,0,1,0,0,1,0,0,0, 0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0, 0,0,0,0};
541
  const int carresmod63[]={
542
    1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0, 0,0,1,0,0,1,0,0,1,0,
543
    0,0,0,0,0,0,1,1,0,0, 0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0, 0,0,0};
544
  const int carresmod65[]={
545
    1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,1,0,0,1,
546
    1,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0, 0,1,0,0,1};
547
  const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
548
  return (carresmod64[A & 0x3fUL]
549
    && carresmod63[A % 63UL]
550
    && carresmod65[A % 65UL]
551
    && carresmod11[A % 11UL]);
552
}
553

554
/* emulate Z_issquareall on single-word integers */
555
long
556
uissquareall(ulong A, ulong *sqrtA)
557
{
558
  if (!A) { *sqrtA = 0; return 1; }
559
  if (carremod(A))
560
  {
561
    ulong a = usqrt(A);
562
    if (a * a == A) { *sqrtA = a; return 1; }
563
  }
564
  return 0;
565
}
566
long
567
uissquare(ulong A)
568
{
569
  if (!A) return 1;
570
  if (carremod(A))
571
  {
572
    ulong a = usqrt(A);
573
    if (a * a == A) return 1;
574
  }
575
  return 0;
576
}
577

578
long
579
Z_issquareall(GEN x, GEN *pt)
580
{
581
  pari_sp av;
582
  GEN y, r;
583

584
  switch(signe(x))
585
  {
586
    case -1: return 0;
587
    case 0: if (pt) *pt=gen_0; return 1;
588
  }
589
  if (lgefint(x) == 3)
590
  {
591
    ulong u = uel(x,2), a;
592
    if (!pt) return uissquare(u);
593
    if (!uissquareall(u, &a)) return 0;
594
    *pt = utoipos(a); return 1;
595
  }
596
  if (!carremod(umodiu(x, 64*63*65*11))) return 0;
597
  av = avma; y = sqrtremi(x, &r);
598
  if (r != gen_0) return gc_long(av,0);
599
  if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
600
  return 1;
601
}
602

603
/* a t_INT, p prime */
604
long
605
Zp_issquare(GEN a, GEN p)
606
{
607
  long v;
608
  GEN ap;
609

610
  if (!signe(a) || gequal1(a)) return 1;
611
  v = Z_pvalrem(a, p, &ap);
612
  if (v&1) return 0;
613
  return absequaliu(p, 2)? umodiu(ap, 8) == 1
614
                      : kronecker(ap,p) == 1;
615
}
616

617
static long
618
polissquareall(GEN x, GEN *pt)
619
{
620
  pari_sp av;
621
  long v;
622
  GEN y, a, b, p;
623

624
  if (!signe(x))
625
  {
626
    if (pt) *pt = gcopy(x);
627
    return 1;
628
  }
629
  if (odd(degpol(x))) return 0; /* odd degree */
630
  av = avma;
631
  v = RgX_valrem(x, &x);
632
  if (v & 1) return gc_long(av,0);
633
  a = gel(x,2); /* test constant coeff */
634
  if (!pt)
635
  { if (!issquare(a)) return gc_long(av,0); }
636
  else
637
  { if (!issquareall(a,&b)) return gc_long(av,0); }
638
  if (!degpol(x)) { /* constant polynomial */
639
    if (!pt) return gc_long(av,1);
640
    y = scalarpol(b, varn(x)); goto END;
641
  }
642
  p = characteristic(x);
643
  if (signe(p) && !mod2(p))
644
  {
645
    long i, lx;
646
    if (!absequaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
647
    x = gmul(x, mkintmod(gen_1, gen_2));
648
    lx = lg(x);
649
    if ((lx-3) & 1) return gc_long(av,0);
650
    for (i = 3; i < lx; i+=2)
651
      if (!gequal0(gel(x,i))) return gc_long(av,0);
652
    if (pt) {
653
      y = cgetg((lx+3) / 2, t_POL);
654
      for (i = 2; i < lx; i+=2)
655
        if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) return gc_long(av,0);
656
      y[1] = evalsigne(1) | evalvarn(varn(x));
657
      goto END;
658
    } else {
659
      for (i = 2; i < lx; i+=2)
660
        if (!issquare(gel(x,i))) return gc_long(av,0);
661
      return gc_long(av,1);
662
    }
663
  }
664
  else
665
  {
666
    long m = 1;
667
    x = RgX_Rg_div(x,a);
668
    /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
669
    if (!signe(p)) x = RgX_deflate_max(x,&m);
670
    y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
671
    if (!RgX_equal(RgX_sqr(y), x)) return gc_long(av,0);
672
    if (!pt) return gc_long(av,1);
673
    if (!gequal1(a)) y = gmul(b, y);
674
    if (m != 1) y = RgX_inflate(y,m);
675
  }
676
END:
677
  if (v) y = RgX_shift_shallow(y, v>>1);
678
  *pt = gerepilecopy(av, y); return 1;
679
}
680

681
/* b unit mod p */
682
static int
683
Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
684
{
685
  if (d == 1)
686
  { /* mod p: faster */
687
    if (!Fp_ispower(b, K, p)) return 0;
688
    if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
689
  }
690
  else
691
  { /* mod p^{2 +} */
692
    if (!ispower(cvtop(b, p, d), K, pt)) return 0;
693
    if (pt) *pt = gtrunc(*pt);
694
  }
695
  return 1;
696
}
697

698
/* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
699
 * Decide mod p^e, then reduce a mod q unless q = NULL. */
700
static int
701
handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
702
{
703
  GEN t, A;
704
  long v = Z_pvalrem(*pa, p, &A), d = e - v;
705
  if (d <= 0) t = gen_0;
706
  else
707
  {
708
    ulong r;
709
    v = uabsdivui_rem(v, K, &r);
710
    if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
711
    if (L && v) t = mulii(t, powiu(p, v));
712
  }
713
  if (q) *pa = modii(*pa, q);
714
  if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
715
  return 1;
716
}
717
long
718
Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
719
{
720
  GEN L, N;
721
  pari_sp av;
722
  long e, i, l;
723
  ulong pp;
724
  forprime_t S;
725

726
  if (!signe(a))
727
  {
728
    if (pt) {
729
      GEN t = cgetg(3, t_INTMOD);
730
      gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
731
    }
732
    return 1;
733
  }
734
  /* a != 0 */
735
  av = avma;
736

737
  if (typ(q) != t_INT) /* integer factorization */
738
  {
739
    GEN P = gel(q,1), E = gel(q,2);
740
    l = lg(P);
741
    L = pt? vectrunc_init(l): NULL;
742
    for (i = 1; i < l; i++)
743
    {
744
      GEN p = gel(P,i);
745
      long e = itos(gel(E,i));
746
      if (!handle_pe(&a, NULL, L, K, p, e)) return gc_long(av,0);
747
    }
748
    goto END;
749
  }
750
  if (!mod2(K)
751
      && kronecker(a, shifti(q,-vali(q))) == -1) return gc_long(av,0);
752
  L = pt? vectrunc_init(expi(q)+1): NULL;
753
  u_forprime_init(&S, 2, tridiv_bound(q));
754
  while ((pp = u_forprime_next(&S)))
755
  {
756
    int stop;
757
    e = Z_lvalrem_stop(&q, pp, &stop);
758
    if (!e) continue;
759
    if (!handle_pe(&a, q, L, K, utoipos(pp), e)) return gc_long(av,0);
760
    if (stop)
761
    {
762
      if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
763
      goto END;
764
    }
765
  }
766
  l = lg(primetab);
767
  for (i = 1; i < l; i++)
768
  {
769
    GEN p = gel(primetab,i);
770
    e = Z_pvalrem(q, p, &q);
771
    if (!e) continue;
772
    if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
773
    if (is_pm1(q)) goto END;
774
  }
775
  N = gcdii(a,q);
776
  if (!is_pm1(N))
777
  {
778
    if (ifac_isprime(N))
779
    {
780
      e = Z_pvalrem(q, N, &q);
781
      if (!handle_pe(&a, q, L, K, N, e)) return gc_long(av,0);
782
    }
783
    else
784
    {
785
      GEN part = ifac_start(N, 0);
786
      for(;;)
787
      {
788
        long e;
789
        GEN p;
790
        if (!ifac_next(&part, &p, &e)) break;
791
        e = Z_pvalrem(q, p, &q);
792
        if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
793
      }
794
    }
795
  }
796
  if (!is_pm1(q))
797
  {
798
    if (ifac_isprime(q))
799
    {
800
      if (!handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
801
    }
802
    else
803
    {
804
      GEN part = ifac_start(q, 0);
805
      for(;;)
806
      {
807
        long e;
808
        GEN p;
809
        if (!ifac_next(&part, &p, &e)) break;
810
        if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
811
      }
812
    }
813
  }
814
END:
815
  if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
816
  return 1;
817
}
818

819
static long
820
polmodispower(GEN x, GEN K, GEN *pt)
821
{
822
  pari_sp av = avma;
823
  GEN p = NULL, T = NULL;
824
  if (Rg_is_FpXQ(x, &T,&p) && p)
825
  {
826
    x = liftall_shallow(x);
827
    if (T) T = liftall_shallow(T);
828
    if (!Fq_ispower(x, K, T, p)) return gc_long(av,0);
829
    if (!pt) return gc_long(av,1);
830
    x = Fq_sqrtn(x, K, T,p, NULL);
831
    if (typ(x) == t_INT)
832
      x = Fp_to_mod(x,p);
833
    else
834
      x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
835
    *pt = gerepilecopy(av, x); return 1;
836
  }
837
  pari_err_IMPL("ispower for general t_POLMOD");
838
  return 0;
839
}
840
static long
841
rfracispower(GEN x, GEN K, GEN *pt)
842
{
843
  pari_sp av = avma;
844
  GEN n = gel(x,1), d = gel(x,2);
845
  long v = -RgX_valrem(d, &d), vx = varn(d);
846
  if (typ(n) == t_POL && varn(n) == vx) v += RgX_valrem(n, &n);
847
  if (!dvdsi(v, K)) return gc_long(av, 0);
848
  if (lg(d) >= 3)
849
  {
850
    GEN a = gel(d,2); /* constant term */
851
    if (!gequal1(a)) { d = RgX_Rg_div(d, a); n = gdiv(n, a); }
852
  }
853
  if (!ispower(d, K, pt? &d: NULL) || !ispower(n, K, pt? &n: NULL))
854
    return gc_long(av, 0);
855
  if (!pt) return gc_long(av, 1);
856
  x = gdiv(n, d);
857
  if (v) x = gmul(x, monomial(gen_1, v / itos(K), vx));
858
  *pt = gerepileupto(av, x); return 1;
859
}
860
long
861
issquareall(GEN x, GEN *pt)
862
{
863
  long tx = typ(x);
864
  GEN F;
865
  pari_sp av;
866

867
  if (!pt) return issquare(x);
868
  switch(tx)
869
  {
870
    case t_INT: return Z_issquareall(x, pt);
871
    case t_FRAC: av = avma;
872
      F = cgetg(3, t_FRAC);
873
      if (   !Z_issquareall(gel(x,1), &gel(F,1))
874
          || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
875
      *pt = F; return 1;
876

877
    case t_POLMOD:
878
      return polmodispower(x, gen_2, pt);
879
    case t_POL: return polissquareall(x,pt);
880
    case t_RFRAC: return rfracispower(x, gen_2, pt);
881

882
    case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
883
      if (!issquare(x)) return 0;
884
      *pt = gsqrt(x, DEFAULTPREC); return 1;
885

886
    case t_INTMOD:
887
      return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
888

889
    case t_FFELT: return FF_issquareall(x, pt);
890

891
  }
892
  pari_err_TYPE("issquareall",x);
893
  return 0; /* LCOV_EXCL_LINE */
894
}
895

896
long
897
issquare(GEN x)
898
{
899
  GEN a, p;
900
  long v;
901

902
  switch(typ(x))
903
  {
904
    case t_INT:
905
      return Z_issquare(x);
906

907
    case t_REAL:
908
      return (signe(x)>=0);
909

910
    case t_INTMOD:
911
      return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
912

913
    case t_FRAC:
914
      return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
915

916
    case t_FFELT: return FF_issquareall(x, NULL);
917

918
    case t_COMPLEX:
919
      return 1;
920

921
    case t_PADIC:
922
      a = gel(x,4); if (!signe(a)) return 1;
923
      if (valp(x)&1) return 0;
924
      p = gel(x,2);
925
      if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
926

927
      v = precp(x); /* here p=2, a is odd */
928
      if ((v>=3 && mod8(a) != 1 ) ||
929
          (v==2 && mod4(a) != 1)) return 0;
930
      return 1;
931

932
    case t_POLMOD:
933
      return polmodispower(x, gen_2, NULL);
934

935
    case t_POL:
936
      return polissquareall(x,NULL);
937

938
    case t_SER:
939
      if (!signe(x)) return 1;
940
      if (valp(x)&1) return 0;
941
      return issquare(gel(x,2));
942

943
    case t_RFRAC:
944
      return rfracispower(x, gen_2, NULL);
945
  }
946
  pari_err_TYPE("issquare",x);
947
  return 0; /* LCOV_EXCL_LINE */
948
}
949
GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
950
GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
951

952
long
953
ispolygonal(GEN x, GEN S, GEN *N)
954
{
955
  pari_sp av = avma;
956
  GEN D, d, n;
957
  if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
958
  if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
959
  if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
960
  if (signe(x) < 0) return 0;
961
  if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
962
  if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
963
  /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
964
  if (abscmpiu(S, 1<<16) < 0) /* common case ! */
965
  {
966
    ulong s = S[2], r;
967
    if (s == 4) return Z_issquareall(x, N);
968
    if (s == 3)
969
      D = addiu(shifti(x, 3), 1);
970
    else
971
      D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
972
    if (!Z_issquareall(D, &d)) return gc_long(av,0);
973
    if (s == 3)
974
      d = subiu(d, 1);
975
    else
976
      d = addiu(d, s - 4);
977
    n = absdiviu_rem(d, 2*s - 4, &r);
978
    if (r) return gc_long(av,0);
979
  }
980
  else
981
  {
982
    GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
983
    D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
984
    if (!Z_issquareall(D, &d)) return gc_long(av,0);
985
    d = addii(d, S_4);
986
    n = dvmdii(shifti(d,-1), S_2, &r);
987
    if (r != gen_0) return gc_long(av,0);
988
  }
989
  if (N) *N = gerepileuptoint(av, n); else set_avma(av);
990
  return 1;
991
}
992

993
/*********************************************************************/
994
/**                                                                 **/
995
/**                        PERFECT POWER                            **/
996
/**                                                                 **/
997
/*********************************************************************/
998
static long
999
polispower(GEN x, GEN K, GEN *pt)
1000
{
1001
  pari_sp av;
1002
  long v, d, k = itos(K);
1003
  GEN y, a, b;
1004
  GEN T = NULL, p = NULL;
1005

1006
  if (!signe(x))
1007
  {
1008
    if (pt) *pt = gcopy(x);
1009
    return 1;
1010
  }
1011
  d = degpol(x);
1012
  if (d % k) return 0; /* degree not multiple of k */
1013
  av = avma;
1014
  if (RgX_is_FpXQX(x, &T, &p) && p)
1015
  { /* over Fq */
1016
    if (T && typ(T) == t_FFELT)
1017
    {
1018
      if (!FFX_ispower(x, k, T, pt)) return gc_long(av,0);
1019
      return 1;
1020
    }
1021
    x = RgX_to_FqX(x,T,p);
1022
    if (!FqX_ispower(x, k, T,p, pt)) return gc_long(av,0);
1023
    if (pt) *pt = gerepileupto(av, FqX_to_mod(*pt, T, p));
1024
    return 1;
1025
  }
1026
  v = RgX_valrem(x, &x);
1027
  if (v % k) return 0;
1028
  v /= k;
1029
  a = gel(x,2); b = NULL;
1030
  if (!ispower(a, K, &b)) return gc_long(av,0);
1031
  if (d)
1032
  {
1033
    GEN p = characteristic(x);
1034
    a = leading_coeff(x);
1035
    if (!ispower(a, K, &b)) return gc_long(av,0);
1036
    x = RgX_normalize(x);
1037
    if (signe(p) && cmpii(p,K) <= 0)
1038
      pari_err_IMPL("ispower(general t_POL) in small characteristic");
1039
    y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1040
    if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
1041
  }
1042
  else
1043
    y = pol_1(varn(x));
1044
  if (pt)
1045
  {
1046
    if (!gequal1(a))
1047
    {
1048
      if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1049
      y = gmul(b,y);
1050
    }
1051
    if (v) y = RgX_shift_shallow(y, v);
1052
    *pt = gerepilecopy(av, y);
1053
  }
1054
  else set_avma(av);
1055
  return 1;
1056
}
1057

1058
long
1059
Z_ispowerall(GEN x, ulong k, GEN *pt)
1060
{
1061
  long s = signe(x);
1062
  ulong mask;
1063
  if (!s) { if (pt) *pt = gen_0; return 1; }
1064
  if (s > 0) {
1065
    if (k == 2) return Z_issquareall(x, pt);
1066
    if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1067
    if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1068
    if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1069
    return is_kth_power(x, k, pt);
1070
  }
1071
  if (!odd(k)) return 0;
1072
  if (Z_ispowerall(absi_shallow(x), k, pt))
1073
  {
1074
    if (pt) *pt = negi(*pt);
1075
    return 1;
1076
  };
1077
  return 0;
1078
}
1079

1080
/* is x a K-th power mod p ? Assume p prime. */
1081
int
1082
Fp_ispower(GEN x, GEN K, GEN p)
1083
{
1084
  pari_sp av = avma;
1085
  GEN p_1;
1086
  x = modii(x, p);
1087
  if (!signe(x) || equali1(x)) return gc_bool(av,1);
1088
  /* implies p > 2 */
1089
  p_1 = subiu(p,1);
1090
  K = gcdii(K, p_1);
1091
  if (absequaliu(K, 2)) return gc_bool(av, kronecker(x,p) > 0);
1092
  x = Fp_pow(x, diviiexact(p_1,K), p);
1093
  return gc_bool(av, equali1(x));
1094
}
1095

1096
/* x unit defined modulo 2^e, e > 0, p prime */
1097
static int
1098
U2_issquare(GEN x, long e)
1099
{
1100
  long r = signe(x)>=0?mod8(x):8-mod8(x);
1101
  if (e==1) return 1;
1102
  if (e==2) return (r&3L) == 1;
1103
  return r == 1;
1104
}
1105
/* x unit defined modulo p^e, e > 0, p prime */
1106
static int
1107
Up_issquare(GEN x, GEN p, long e)
1108
{ return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1109

1110
long
1111
Zn_issquare(GEN d, GEN fn)
1112
{
1113
  long j, np;
1114
  if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1115
  if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1116
  /* integer factorization */
1117
  np = nbrows(fn);
1118
  for (j = 1; j <= np; ++j)
1119
  {
1120
    GEN  r, p = gcoeff(fn, j, 1);
1121
    long e = itos(gcoeff(fn, j, 2));
1122
    long v = Z_pvalrem(d,p,&r);
1123
    if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
1124
  }
1125
  return 1;
1126
}
1127

1128
/* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
1129
GEN
1130
Zn_quad_roots(GEN N, GEN B, GEN C)
1131
{
1132
  pari_sp av = avma;
1133
  GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
1134
  long l, i, j, ct;
1135

1136
  if ((fa = check_arith_non0(N,"Zn_quad_roots")))
1137
  {
1138
    N = typ(N) == t_VEC? gel(N,1): factorback(N);
1139
    fa = clean_Z_factor(fa);
1140
  }
1141
  N = absi_shallow(N);
1142
  N4 = shifti(N,2);
1143
  D = modii(subii(sqri(B), shifti(C,2)), N4);
1144
  if (!signe(D))
1145
  { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
1146
    if (!fa) fa = Z_factor(N);
1147
    P = gel(fa,1);
1148
    E = ZV_to_zv(gel(fa,2));
1149
    l = lg(P);
1150
    for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
1151
    Np = factorback2(P, E); /* x = -B mod N' */
1152
    B = shifti(B,-1);
1153
    return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
1154
  }
1155
  if (!fa)
1156
    fa = Z_factor(N4);
1157
  else  /* convert to factorization of N4 = 4*N */
1158
    fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
1159
  P = gel(fa,1); l = lg(P);
1160
  E = ZV_to_zv(gel(fa,2));
1161
  F = cgetg(l, t_VEC);
1162
  mF= cgetg(l, t_VEC); F0 = gen_0;
1163
  Q = cgetg(l, t_VEC); Q0 = gen_1;
1164
  for (i = j = 1, ct = 0; i < l; i++)
1165
  {
1166
    GEN p = gel(P,i), q, f, mf, D0;
1167
    long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
1168
    if (d <= 0)
1169
    {
1170
      q = powiu(p, (s+1)>>1);
1171
      Q0 = mulii(Q0, q); continue;
1172
    }
1173
    /* d > 0 */
1174
    if (odd(t)) return NULL;
1175
    t2 = t >> 1;
1176
    if (i > 1)
1177
    { /* p > 2 */
1178
      if (kronecker(D0, p) == -1) return NULL;
1179
      q = powiu(p, s - t2);
1180
      f = Zp_sqrt(D0, p, d);
1181
      if (!f) return NULL; /* p was not actually prime... */
1182
      if (t2) f = mulii(powiu(p,t2), f);
1183
      mf = Fp_neg(f, q);
1184
    }
1185
    else
1186
    { /* p = 2 */
1187
      if (d <= 3)
1188
      {
1189
        if (d == 3 && Mod8(D0) != 1) return NULL;
1190
        if (d == 2 && Mod4(D0) != 1) return NULL;
1191
        Q0 = int2n(1+t2); F0 = NULL; continue;
1192
      }
1193
      if (Mod8(D0) != 1) return NULL;
1194
      q = int2n(d - 1 + t2);
1195
      f = Z2_sqrt(D0, d);
1196
      if (t2) f = shifti(f, t2);
1197
      mf = Fp_neg(f, q);
1198
    }
1199
    gel(Q,j) = q;
1200
    gel(F,j) = f;
1201
    gel(mF,j)= mf; j++;
1202
  }
1203
  setlg(Q,j);
1204
  setlg(F,j);
1205
  setlg(mF,j);
1206
  if (is_pm1(Q0)) A = leafcopy(F);
1207
  else
1208
  { /* append the fixed congruence (F0 mod Q0) */
1209
    if (!F0) F0 = shifti(Q0,-1);
1210
    A = shallowconcat(F, F0);
1211
    Q = shallowconcat(Q, Q0);
1212
  }
1213
  ct = 1 << (j-1);
1214
  T = ZV_producttree(Q);
1215
  R = ZV_chinesetree(Q,T);
1216
  Np = gmael(T, lg(T)-1, 1);
1217
  B = modii(B, Np);
1218
  if (!signe(B)) B = NULL;
1219
  Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
1220
  w = cgetg(3, t_VEC);
1221
  gel(w,1) = icopy(Np);
1222
  gel(w,2) = v = cgetg(ct+1, t_VEC);
1223
  l = lg(F);
1224
  for (j = 1; j <= ct; j++)
1225
  {
1226
    pari_sp av2 = avma;
1227
    long m = j - 1;
1228
    GEN u;
1229
    for (i = 1; i < l; i++)
1230
    {
1231
      gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
1232
      m >>= 1;
1233
    }
1234
    u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
1235
    if (B) u = subii(u,B);
1236
    gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
1237
  }
1238
  return gerepileupto(av, w);
1239
}
1240

1241
static long
1242
Qp_ispower(GEN x, GEN K, GEN *pt)
1243
{
1244
  pari_sp av = avma;
1245
  GEN z = Qp_sqrtn(x, K, NULL);
1246
  if (!z) return gc_long(av,0);
1247
  if (pt) *pt = z;
1248
  return 1;
1249
}
1250

1251
long
1252
ispower(GEN x, GEN K, GEN *pt)
1253
{
1254
  GEN z;
1255

1256
  if (!K) return gisanypower(x, pt);
1257
  if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1258
  if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1259
  if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1260
  switch(typ(x)) {
1261
    case t_INT:
1262
      if (lgefint(K) != 3) return 0;
1263
      return Z_ispowerall(x, itou(K), pt);
1264
    case t_FRAC:
1265
    {
1266
      GEN a = gel(x,1), b = gel(x,2);
1267
      ulong k;
1268
      if (lgefint(K) != 3) return 0;
1269
      k = itou(K);
1270
      if (pt) {
1271
        z = cgetg(3, t_FRAC);
1272
        if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1273
          *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1274
        }
1275
        set_avma((pari_sp)(z + 3)); return 0;
1276
      }
1277
      return Z_ispower(a, k) && Z_ispower(b, k);
1278
    }
1279
    case t_INTMOD:
1280
      return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1281
    case t_FFELT:
1282
      return FF_ispower(x, K, pt);
1283

1284
    case t_PADIC:
1285
      return Qp_ispower(x, K, pt);
1286
    case t_POLMOD:
1287
      return polmodispower(x, K, pt);
1288
    case t_POL:
1289
      return polispower(x, K, pt);
1290
    case t_RFRAC:
1291
      return rfracispower(x, K, pt);
1292
    case t_REAL:
1293
      if (signe(x) < 0 && !mpodd(K)) return 0;
1294
    case t_COMPLEX:
1295
      if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1296
      return 1;
1297

1298
    case t_SER:
1299
      if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
1300
        return 0;
1301
      if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1302
      return 1;
1303
  }
1304
  pari_err_TYPE("ispower",x);
1305
  return 0; /* LCOV_EXCL_LINE */
1306
}
1307

1308
long
1309
gisanypower(GEN x, GEN *pty)
1310
{
1311
  long tx = typ(x);
1312
  ulong k, h;
1313
  if (tx == t_INT) return Z_isanypower(x, pty);
1314
  if (tx == t_FRAC)
1315
  {
1316
    pari_sp av = avma;
1317
    GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1318
    long i, j, p, e;
1319
    int sw = (abscmpii(a, b) > 0);
1320

1321
    if (sw) swap(a, b);
1322
    k = Z_isanypower(a, pty? &a: NULL);
1323
    if (!k)
1324
    { /* a = -1,1 or not a pure power */
1325
      if (!is_pm1(a)) return gc_long(av,0);
1326
      if (signe(a) < 0) b = negi(b);
1327
      k = Z_isanypower(b, pty? &b: NULL);
1328
      if (!k || !pty) return gc_long(av,k);
1329
      *pty = gerepileupto(av, ginv(b));
1330
      return k;
1331
    }
1332
    fa = factoru(k);
1333
    P = gel(fa,1);
1334
    E = gel(fa,2); h = k;
1335
    for (i = lg(P) - 1; i > 0; i--)
1336
    {
1337
      p = P[i];
1338
      e = E[i];
1339
      for (j = 0; j < e; j++)
1340
        if (!is_kth_power(b, p, &b)) break;
1341
      if (j < e) k /= upowuu(p, e - j);
1342
    }
1343
    if (k == 1) return gc_long(av,0);
1344
    if (!pty) return gc_long(av,k);
1345
    if (k != h) a = powiu(a, h/k);
1346
    *pty = gerepilecopy(av, mkfrac(a, b));
1347
    return k;
1348
  }
1349
  pari_err_TYPE("gisanypower", x);
1350
  return 0; /* LCOV_EXCL_LINE */
1351
}
1352

1353
/* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1354
 * No need to optimize for 2,3,5,7 powers (done before) */
1355
static long
1356
split_exponent(ulong e, GEN *x)
1357
{
1358
  GEN fa, P, E;
1359
  long i, j, l, k = 1;
1360
  if (e == 1) return 1;
1361
  fa = factoru(e);
1362
  P = gel(fa,1);
1363
  E = gel(fa,2); l = lg(P);
1364
  for (i = 1; i < l; i++)
1365
  {
1366
    ulong p = P[i];
1367
    for (j = 0; j < E[i]; j++)
1368
    {
1369
      GEN y;
1370
      if (!is_kth_power(*x, p, &y)) break;
1371
      k *= p; *x = y;
1372
    }
1373
  }
1374
  return k;
1375
}
1376

1377
static long
1378
Z_isanypower_nosmalldiv(GEN *px)
1379
{ /* any prime divisor of x is > 102 */
1380
  const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1381
  const double LOG103 = 4.6347; /* lower bound for log(103) */
1382
  forprime_t T;
1383
  ulong mask = 7, e2;
1384
  long k, ex;
1385
  GEN y, x = *px;
1386

1387
  k = 1;
1388
  while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1389
  while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1390
  e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1391
  if (u_forprime_init(&T, 11, e2))
1392
  {
1393
    GEN logx = NULL;
1394
    const ulong Q = 30011; /* prime */
1395
    ulong p, xmodQ;
1396
    double dlogx = 0;
1397
    /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1398
     * for large p the modular checks are no longer competitively fast */
1399
    while ( (ex = is_pth_power(x, &y, &T, 30)) )
1400
    {
1401
      k *= ex; x = y;
1402
      e2 = (ulong)((expi(x) + 1) / LOG2_103);
1403
      u_forprime_restrict(&T, e2);
1404
    }
1405
    if (DEBUGLEVEL>4)
1406
      err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
1407
    xmodQ = umodiu(x, Q);
1408
    /* test Q | x, just in case */
1409
    if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1410
    /* x^(1/p) < 2^31 */
1411
    p = T.p;
1412
    if (p <= e2)
1413
    {
1414
      logx = logr_abs( itor(x, DEFAULTPREC) );
1415
      dlogx = rtodbl(logx);
1416
      e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1417
    }
1418
    while (p && p <= e2)
1419
    { /* is x a p-th power ? By computing y = round(x^(1/p)).
1420
       * Check whether y^p = x, first mod Q, then exactly. */
1421
      pari_sp av = avma;
1422
      long e;
1423
      GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1424
      ulong ymodQ = umodiu(y,Q);
1425
      if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1426
                   || !equalii(powiu(y, p), x)) set_avma(av);
1427
      else
1428
      {
1429
        k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1430
        e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1431
        u_forprime_restrict(&T, e2);
1432
        continue; /* if success, retry same p */
1433
      }
1434
      p = u_forprime_next(&T);
1435
    }
1436
  }
1437
  *px = x; return k;
1438
}
1439

1440
static ulong tinyprimes[] = {
1441
  2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1442
  73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1443
  157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1444
};
1445

1446
/* disregard the sign of x, caller will take care of x < 0 */
1447
static long
1448
Z_isanypower_aux(GEN x, GEN *pty)
1449
{
1450
  long ex, v, i, l, k;
1451
  GEN y, P, E;
1452
  ulong mask, e = 0;
1453

1454
  if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1455

1456
  if (signe(x) < 0) x = negi(x);
1457
  k = l = 1;
1458
  P = cgetg(26 + 1, t_VECSMALL);
1459
  E = cgetg(26 + 1, t_VECSMALL);
1460
  /* trial division */
1461
  for(i = 0; i < 26; i++)
1462
  {
1463
    ulong p = tinyprimes[i];
1464
    int stop;
1465
    v = Z_lvalrem_stop(&x, p, &stop);
1466
    if (v)
1467
    {
1468
      P[l] = p;
1469
      E[l] = v; l++;
1470
      e = ugcd(e, v); if (e == 1) goto END;
1471
    }
1472
    if (stop) {
1473
      if (is_pm1(x)) k = e;
1474
      goto END;
1475
    }
1476
  }
1477

1478
  if (e)
1479
  { /* Bingo. Result divides e */
1480
    long v3, v5, v7;
1481
    ulong e2 = e;
1482
    v = u_lvalrem(e2, 2, &e2);
1483
    if (v)
1484
    {
1485
      for (i = 0; i < v; i++)
1486
      {
1487
        if (!Z_issquareall(x, &y)) break;
1488
        k <<= 1; x = y;
1489
      }
1490
    }
1491
    mask = 0;
1492
    v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1493
    v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1494
    v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1495
    while ( (ex = is_357_power(x, &y, &mask)) ) {
1496
      x = y;
1497
      switch(ex)
1498
      {
1499
        case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1500
        case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1501
        case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1502
      }
1503
    }
1504
    k *= split_exponent(e2, &x);
1505
  }
1506
  else
1507
    k = Z_isanypower_nosmalldiv(&x);
1508
END:
1509
  if (pty && k != 1)
1510
  {
1511
    if (e)
1512
    { /* add missing small factors */
1513
      y = powuu(P[1], E[1] / k);
1514
      for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1515
      x = equali1(x)? y: mulii(x,y);
1516
    }
1517
    *pty = x;
1518
  }
1519
  return k == 1? 0: k;
1520
}
1521

1522
long
1523
Z_isanypower(GEN x, GEN *pty)
1524
{
1525
  pari_sp av = avma;
1526
  long k = Z_isanypower_aux(x, pty);
1527
  if (!k) return gc_long(av,0);
1528
  if (signe(x) < 0)
1529
  {
1530
    long v = vals(k);
1531
    if (v)
1532
    {
1533
      GEN y;
1534
      k >>= v;
1535
      if (k == 1) return gc_long(av,0);
1536
      if (!pty) return gc_long(av,k);
1537
      y = *pty;
1538
      y = powiu(y, 1<<v);
1539
      togglesign(y);
1540
      *pty = gerepileuptoint(av, y);
1541
      return k;
1542
    }
1543
    if (pty) togglesign_safe(pty);
1544
  }
1545
  if (pty) *pty = gerepilecopy(av, *pty); else set_avma(av);
1546
  return k;
1547
}
1548

1549
/* Faster than */
1550
/*   !cmpii(n, int2n(vali(n))) */
1551
/*   !cmpis(shifti(n, -vali(n)), 1) */
1552
/*   expi(n) == vali(n) */
1553
/*   hamming(n) == 1 */
1554
/* even for single-word values, and much faster for multiword values. */
1555
/* If all you have is a word, you can just use n & !(n & (n-1)). */
1556
long
1557
Z_ispow2(GEN n)
1558
{
1559
  GEN xp;
1560
  long i, lx;
1561
  ulong u;
1562
  if (signe(n) != 1) return 0;
1563
  xp = int_LSW(n);
1564
  lx = lgefint(n);
1565
  u = *xp;
1566
  for (i = 3; i < lx; ++i)
1567
  {
1568
    if (u) return 0;
1569
    xp = int_nextW(xp);
1570
    u = *xp;
1571
  }
1572
  return !(u & (u-1));
1573
}
1574

1575
static long
1576
isprimepower_i(GEN n, GEN *pt, long flag)
1577
{
1578
  pari_sp av = avma;
1579
  long i, v;
1580

1581
  if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1582
  if (signe(n) <= 0) return 0;
1583

1584
  if (lgefint(n) == 3)
1585
  {
1586
    ulong p;
1587
    v = uisprimepower(n[2], &p);
1588
    if (v)
1589
    {
1590
      if (pt) *pt = utoipos(p);
1591
      return v;
1592
    }
1593
    return 0;
1594
  }
1595
  for (i = 0; i < 26; i++)
1596
  {
1597
    ulong p = tinyprimes[i];
1598
    v = Z_lvalrem(n, p, &n);
1599
    if (v)
1600
    {
1601
      set_avma(av);
1602
      if (!is_pm1(n)) return 0;
1603
      if (pt) *pt = utoipos(p);
1604
      return v;
1605
    }
1606
  }
1607
  /* p | n => p >= 103 */
1608
  v = Z_isanypower_nosmalldiv(&n); /* expensive */
1609
  if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
1610
  if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
1611
  return v;
1612
}
1613
long
1614
isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1615
long
1616
ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1617

1618
long
1619
uisprimepower(ulong n, ulong *pp)
1620
{ /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1621
   * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1622
  const ulong CUTOFF = 200UL;
1623
  const long TINYCUTOFF = 46;  /* tinyprimes[45] = 199 */
1624
  const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1625
#ifdef LONG_IS_64BIT
1626
  /* primes preceeding the appropriate root of ULONG_MAX. */
1627
  const ulong ROOT9 = 137;
1628
  const ulong ROOT8 = 251;
1629
  const ulong ROOT7 = 563;
1630
  const ulong ROOT5 = 7129;
1631
  const ulong ROOT4 = 65521;
1632
#else
1633
  const ulong ROOT9 = 11;
1634
  const ulong ROOT8 = 13;
1635
  const ulong ROOT7 = 23;
1636
  const ulong ROOT5 = 83;
1637
  const ulong ROOT4 = 251;
1638
#endif
1639
  ulong mask;
1640
  long v, i;
1641
  int e;
1642
  if (n < 2) return 0;
1643
  if (!odd(n)) {
1644
    if (n & (n-1)) return 0;
1645
    *pp = 2; return vals(n);
1646
  }
1647
  if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1648
  for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1649
  {
1650
    ulong p = tinyprimes[i];
1651
    if (n % p == 0)
1652
    {
1653
      v = u_lvalrem(n, p, &n);
1654
      if (n == 1) { *pp = p; return v; }
1655
      return 0;
1656
    }
1657
  }
1658
  /* p | n => p >= CUTOFF */
1659

1660
  if (n < CUTOFF3)
1661
  {
1662
    if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1663
    if (uissquareall(n, &n)) { *pp = n; return 2; }
1664
    return 0;
1665
  }
1666

1667
  /* Check for squares, fourth powers, and eighth powers as appropriate. */
1668
  v = 1;
1669
  if (uissquareall(n, &n)) {
1670
    v <<= 1;
1671
    if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1672
      v <<= 1;
1673
      if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1674
    }
1675
  }
1676

1677
  if (CUTOFF > ROOT5) mask = 1;
1678
  else
1679
  {
1680
    const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1681
    if (n < CUTOFF5) mask = 1; else mask = 3;
1682
    if (CUTOFF <= ROOT7)
1683
    {
1684
      const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1685
      if (n >= CUTOFF7) mask = 7;
1686
    }
1687
  }
1688

1689
  if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1690
  if ((e = uis_357_power(n, &n, &mask))) v *= e;
1691

1692
  if (uisprime_101(n)) { *pp = n; return v; }
1693
  return 0;
1694
}
1695

1696
/*********************************************************************/
1697
/**                                                                 **/
1698
/**                        KRONECKER SYMBOL                         **/
1699
/**                                                                 **/
1700
/*********************************************************************/
1701
/* t = 3,5 mod 8 ?  (= 2 not a square mod t) */
1702
static int
1703
ome(long t)
1704
{
1705
  switch(t & 7)
1706
  {
1707
    case 3:
1708
    case 5: return 1;
1709
    default: return 0;
1710
  }
1711
}
1712
/* t a t_INT, is t = 3,5 mod 8 ? */
1713
static int
1714
gome(GEN t)
1715
{ return signe(t)? ome( mod2BIL(t) ): 0; }
1716

1717
/* assume y odd, return kronecker(x,y) * s */
1718
static long
1719
krouu_s(ulong x, ulong y, long s)
1720
{
1721
  ulong x1 = x, y1 = y, z;
1722
  while (x1)
1723
  {
1724
    long r = vals(x1);
1725
    if (r)
1726
    {
1727
      if (odd(r) && ome(y1)) s = -s;
1728
      x1 >>= r;
1729
    }
1730
    if (x1 & y1 & 2) s = -s;
1731
    z = y1 % x1; y1 = x1; x1 = z;
1732
  }
1733
  return (y1 == 1)? s: 0;
1734
}
1735

1736
long
1737
kronecker(GEN x, GEN y)
1738
{
1739
  pari_sp av = avma;
1740
  long s = 1, r;
1741
  ulong xu;
1742

1743
  if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1744
  if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1745
  switch (signe(y))
1746
  {
1747
    case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1748
    case 0: return is_pm1(x);
1749
  }
1750
  r = vali(y);
1751
  if (r)
1752
  {
1753
    if (!mpodd(x)) return gc_long(av,0);
1754
    if (odd(r) && gome(x)) s = -s;
1755
    y = shifti(y,-r);
1756
  }
1757
  x = modii(x,y);
1758
  while (lgefint(x) > 3) /* x < y */
1759
  {
1760
    GEN z;
1761
    r = vali(x);
1762
    if (r)
1763
    {
1764
      if (odd(r) && gome(y)) s = -s;
1765
      x = shifti(x,-r);
1766
    }
1767
    /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1768
    if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1769
    z = remii(y,x); y = x; x = z;
1770
    if (gc_needed(av,2))
1771
    {
1772
      if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1773
      gerepileall(av, 2, &x, &y);
1774
    }
1775
  }
1776
  xu = itou(x);
1777
  if (!xu) return is_pm1(y)? s: 0;
1778
  r = vals(xu);
1779
  if (r)
1780
  {
1781
    if (odd(r) && gome(y)) s = -s;
1782
    xu >>= r;
1783
  }
1784
  /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1785
  if (xu & mod2BIL(y) & 2) s = -s;
1786
  return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
1787
}
1788

1789
long
1790
krois(GEN x, long y)
1791
{
1792
  ulong yu;
1793
  long s = 1;
1794

1795
  if (y <= 0)
1796
  {
1797
    if (y == 0) return is_pm1(x);
1798
    yu = (ulong)-y; if (signe(x) < 0) s = -1;
1799
  }
1800
  else
1801
    yu = (ulong)y;
1802
  if (!odd(yu))
1803
  {
1804
    long r;
1805
    if (!mpodd(x)) return 0;
1806
    r = vals(yu); yu >>= r;
1807
    if (odd(r) && gome(x)) s = -s;
1808
  }
1809
  return krouu_s(umodiu(x, yu), yu, s);
1810
}
1811
/* assume y != 0 */
1812
long
1813
kroiu(GEN x, ulong y)
1814
{
1815
  long r;
1816
  if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1817
  if (!mpodd(x)) return 0;
1818
  r = vals(y); y >>= r;
1819
  return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
1820
}
1821

1822
/* assume y > 0, odd, return s * kronecker(x,y) */
1823
static long
1824
krouodd(ulong x, GEN y, long s)
1825
{
1826
  long r;
1827
  if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1828
  if (!x) return 0; /* y != 1 */
1829
  r = vals(x);
1830
  if (r)
1831
  {
1832
    if (odd(r) && gome(y)) s = -s;
1833
    x >>= r;
1834
  }
1835
  /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1836
  if (x & mod2BIL(y) & 2) s = -s;
1837
  return krouu_s(umodiu(y,x), x, s);
1838
}
1839

1840
long
1841
krosi(long x, GEN y)
1842
{
1843
  const pari_sp av = avma;
1844
  long s = 1, r;
1845
  switch (signe(y))
1846
  {
1847
    case -1: y = negi(y); if (x < 0) s = -1; break;
1848
    case 0: return (x==1 || x==-1);
1849
  }
1850
  r = vali(y);
1851
  if (r)
1852
  {
1853
    if (!odd(x)) return gc_long(av,0);
1854
    if (odd(r) && ome(x)) s = -s;
1855
    y = shifti(y,-r);
1856
  }
1857
  if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1858
  return gc_long(av, krouodd((ulong)x, y, s));
1859
}
1860

1861
long
1862
kroui(ulong x, GEN y)
1863
{
1864
  const pari_sp av = avma;
1865
  long s = 1, r;
1866
  switch (signe(y))
1867
  {
1868
    case -1: y = negi(y); break;
1869
    case 0: return x==1UL;
1870
  }
1871
  r = vali(y);
1872
  if (r)
1873
  {
1874
    if (!odd(x)) return gc_long(av,0);
1875
    if (odd(r) && ome(x)) s = -s;
1876
    y = shifti(y,-r);
1877
  }
1878
  return gc_long(av,  krouodd(x, y, s));
1879
}
1880

1881
long
1882
kross(long x, long y)
1883
{
1884
  ulong yu;
1885
  long s = 1;
1886

1887
  if (y <= 0)
1888
  {
1889
    if (y == 0) return (labs(x)==1);
1890
    yu = (ulong)-y; if (x < 0) s = -1;
1891
  }
1892
  else
1893
    yu = (ulong)y;
1894
  if (!odd(yu))
1895
  {
1896
    long r;
1897
    if (!odd(x)) return 0;
1898
    r = vals(yu); yu >>= r;
1899
    if (odd(r) && ome(x)) s = -s;
1900
  }
1901
  x %= (long)yu; if (x < 0) x += yu;
1902
  return krouu_s((ulong)x, yu, s);
1903
}
1904

1905
long
1906
krouu(ulong x, ulong y)
1907
{
1908
  long r;
1909
  if (odd(y)) return krouu_s(x, y, 1);
1910
  if (!odd(x)) return 0;
1911
  r = vals(y); y >>= r;
1912
  return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
1913
}
1914

1915
/*********************************************************************/
1916
/**                                                                 **/
1917
/**                          HILBERT SYMBOL                         **/
1918
/**                                                                 **/
1919
/*********************************************************************/
1920
/* x,y are t_INT or t_REAL */
1921
static long
1922
mphilbertoo(GEN x, GEN y)
1923
{
1924
  long sx = signe(x), sy = signe(y);
1925
  if (!sx || !sy) return 0;
1926
  return (sx < 0 && sy < 0)? -1: 1;
1927
}
1928

1929
long
1930
hilbertii(GEN x, GEN y, GEN p)
1931
{
1932
  pari_sp av;
1933
  long oddvx, oddvy, z;
1934

1935
  if (!p) return mphilbertoo(x,y);
1936
  if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1937
  if (!signe(x) || !signe(y)) return 0;
1938
  av = avma;
1939
  oddvx = odd(Z_pvalrem(x,p,&x));
1940
  oddvy = odd(Z_pvalrem(y,p,&y));
1941
  /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1942
  if (absequaliu(p, 2))
1943
  {
1944
    z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
1945
    if (oddvx && gome(y)) z = -z;
1946
    if (oddvy && gome(x)) z = -z;
1947
  }
1948
  else
1949
  {
1950
    z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
1951
    if (oddvx && kronecker(y,p) < 0) z = -z;
1952
    if (oddvy && kronecker(x,p) < 0) z = -z;
1953
  }
1954
  return gc_long(av, z);
1955
}
1956

1957
static void
1958
err_prec(void) { pari_err_PREC("hilbert"); }
1959
static void
1960
err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1961
static void
1962
err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1963

1964
/* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1965
 * Return lift(x) provided it's p-adic accuracy is large enough to decide
1966
 * hilbert()'s value [ problem at p = 2 ] */
1967
static GEN
1968
lift_intmod(GEN x, GEN *pp)
1969
{
1970
  GEN p = *pp, N = gel(x,1);
1971
  x = gel(x,2);
1972
  if (!p)
1973
  {
1974
    *pp = p = N;
1975
    switch(itos_or_0(p))
1976
    {
1977
      case 2:
1978
      case 4: err_prec();
1979
    }
1980
    return x;
1981
  }
1982
  if (!signe(p)) err_oo(N);
1983
  if (absequaliu(p,2))
1984
  { if (vali(N) <= 2) err_prec(); }
1985
  else
1986
  { if (!dvdii(N,p)) err_p(N,p); }
1987
  if (!signe(x)) err_prec();
1988
  return x;
1989
}
1990
/* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1991
 * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1992
 * to decide hilbert()'s value [ problem at p = 2 ]*/
1993
static GEN
1994
lift_padic(GEN x, GEN *pp)
1995
{
1996
  GEN p = *pp, q = gel(x,2), y = gel(x,4);
1997
  if (!p) *pp = p = q;
1998
  else if (!equalii(p,q)) err_p(p, q);
1999
  if (absequaliu(p,2) && precp(x) <= 2) err_prec();
2000
  if (!signe(y)) err_prec();
2001
  return odd(valp(x))? mulii(p,y): y;
2002
}
2003

2004
long
2005
hilbert(GEN x, GEN y, GEN p)
2006
{
2007
  pari_sp av = avma;
2008
  long tx = typ(x), ty = typ(y);
2009

2010
  if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
2011
  if (tx == t_REAL)
2012
  {
2013
    if (p && signe(p)) err_oo(p);
2014
    switch (ty)
2015
    {
2016
      case t_INT:
2017
      case t_REAL: return mphilbertoo(x,y);
2018
      case t_FRAC: return mphilbertoo(x,gel(y,1));
2019
      default: pari_err_TYPE2("hilbert",x,y);
2020
    }
2021
  }
2022
  if (ty == t_REAL)
2023
  {
2024
    if (p && signe(p)) err_oo(p);
2025
    switch (tx)
2026
    {
2027
      case t_INT:
2028
      case t_REAL: return mphilbertoo(x,y);
2029
      case t_FRAC: return mphilbertoo(gel(x,1),y);
2030
      default: pari_err_TYPE2("hilbert",x,y);
2031
    }
2032
  }
2033
  if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
2034
  if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
2035

2036
  if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
2037
  if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
2038

2039
  if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
2040
  if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
2041

2042
  if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
2043
  if (p && !signe(p)) p = NULL;
2044
  return gc_long(av, hilbertii(x,y,p));
2045
}
2046

2047
/*******************************************************************/
2048
/*                                                                 */
2049
/*                       SQUARE ROOT MODULO p                      */
2050
/*                                                                 */
2051
/*******************************************************************/
2052
static void
2053
checkp(ulong q, ulong p)
2054
{ if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
2055
/* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
2056
static ulong
2057
nonsquare1_Fl(ulong p)
2058
{
2059
  forprime_t S;
2060
  ulong q;
2061
  if ((p & 7UL) != 1) return 2UL;
2062
  q = p % 3; if (q == 2) return 3UL;
2063
  checkp(q, p);
2064
  q = p % 5; if (q == 2 || q == 3) return 5UL;
2065
  checkp(q, p);
2066
  q = p % 7; if (q != 4 && q >= 3) return 7UL;
2067
  checkp(q, p);
2068
  u_forprime_init(&S, 11, p);
2069
  while ((q = u_forprime_next(&S)))
2070
  {
2071
    long i = krouu(q, p);
2072
    if (i < 0) return q;
2073
    checkp(q, p);
2074
  }
2075
  checkp(0, p);
2076
  return 0; /*LCOV_EXCL_LINE*/
2077
}
2078
/* p > 2 a prime */
2079
ulong
2080
nonsquare_Fl(ulong p)
2081
{ return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
2082

2083
ulong
2084
Fl_2gener_pre(ulong p, ulong pi)
2085
{
2086
  ulong p1 = p-1;
2087
  long e = vals(p1);
2088
  if (e == 1) return p1;
2089
  return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
2090
}
2091

2092
/* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
2093
ulong
2094
Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
2095
{
2096
  long i, e, k;
2097
  ulong p1, q, v, w;
2098

2099
  if (!a) return 0;
2100
  p1 = p - 1; e = vals(p1);
2101
  if (e == 0) /* p = 2 */
2102
  {
2103
    if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
2104
    return ((a & 1) == 0)? 0: 1;
2105
  }
2106
  if (e == 1)
2107
  {
2108
    v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
2109
    if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
2110
    p1 = p - v; if (v > p1) v = p1;
2111
    return v;
2112
  }
2113
  q = p1 >> e; /* q = (p-1)/2^oo is odd */
2114
  p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
2115
  if (!p1) return 0;
2116
  v = Fl_mul_pre(a, p1, p, pi);
2117
  w = Fl_mul_pre(v, p1, p, pi);
2118
  if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
2119
  while (w != 1)
2120
  { /* a*w = v^2, y primitive 2^e-th root of 1
2121
       a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2122
    p1 = Fl_sqr_pre(w,p,pi);
2123
    for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
2124
    if (k == e) return ~0UL;
2125
    /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2126
    p1 = y;
2127
    for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
2128
    y = Fl_sqr_pre(p1, p, pi); e = k;
2129
    w = Fl_mul_pre(y, w, p, pi);
2130
    v = Fl_mul_pre(v, p1, p, pi);
2131
  }
2132
  p1 = p - v; if (v > p1) v = p1;
2133
  return v;
2134
}
2135

2136
ulong
2137
Fl_sqrt(ulong a, ulong p)
2138
{
2139
  ulong pi = get_Fl_red(p);
2140
  return Fl_sqrt_pre_i(a, 0, p, pi);
2141
}
2142

2143
ulong
2144
Fl_sqrt_pre(ulong a, ulong p, ulong pi)
2145
{
2146
  return Fl_sqrt_pre_i(a, 0, p, pi);
2147
}
2148

2149
static ulong
2150
Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
2151
{
2152
  ulong x, y, m;
2153
  ulong le1 = upowuu(l, e-1);
2154
  for (x = 2; ; x++)
2155
  {
2156
    y = Fl_powu_pre(x, r, p, pi);
2157
    if (y==1) continue;
2158
    m = Fl_powu_pre(y, le1, p, pi);
2159
    if (m != 1) break;
2160
  }
2161
  *pt_m = m;
2162
  return y;
2163
}
2164

2165
/* solve x^l = a , l prime in G of order q.
2166
 *
2167
 * q =  (l^e)*r, e >= 1, (r,l) = 1
2168
 * y generates the l-Sylow of G
2169
 * m = y^(l^(e-1)) != 1 */
2170
static ulong
2171
Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
2172
{
2173
  ulong p1, v, w, z, dl;
2174
  ulong u2;
2175
  if (a==0) return a;
2176
  u2 = Fl_inv(l%r, r);
2177
  v = Fl_powu_pre(a, u2, p, pi);
2178
  w = Fl_powu_pre(v, l, p, pi);
2179
  w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
2180
  if (w==1) return v;
2181
  if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2182
  while (w!=1)
2183
  {
2184
    ulong k = 0;
2185
    p1 = w;
2186
    do
2187
    {
2188
      z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2189
      if (++k == e) return ULONG_MAX;
2190
    } while (p1!=1);
2191
    dl = Fl_log_pre(z, m, l, p, pi);
2192
    dl = Fl_neg(dl, l);
2193
    p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2194
    m = Fl_powu_pre(m, dl, p, pi);
2195
    e = k;
2196
    v = Fl_mul_pre(p1,v,p,pi);
2197
    y = Fl_powu_pre(p1,l,p,pi);
2198
    w = Fl_mul_pre(y,w,p,pi);
2199
  }
2200
  return v;
2201
}
2202

2203
static ulong
2204
Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
2205
{
2206
  ulong r, e = u_lvalrem(p-1, l, &r);
2207
  return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
2208
}
2209

2210
ulong
2211
Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2212
{
2213
  return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2214
}
2215

2216
ulong
2217
Fl_sqrtl(ulong a, ulong l, ulong p)
2218
{
2219
  ulong pi = get_Fl_red(p);
2220
  return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2221
}
2222

2223
ulong
2224
Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
2225
{
2226
  ulong m, q = p-1, z;
2227
  ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
2228
  if (a==0)
2229
  {
2230
    if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
2231
    if (zetan) *zetan = 1UL;
2232
    return 0;
2233
  }
2234
  if (n==1)
2235
  {
2236
    if (zetan) *zetan = 1;
2237
    return n < 0? Fl_inv(a,p): a;
2238
  }
2239
  if (n==2)
2240
  {
2241
    if (zetan) *zetan = p-1;
2242
    return Fl_sqrt_pre_i(a, 0, p, pi);
2243
  }
2244
  if (a == 1 && !zetan) return a;
2245
  m = ugcd(nn, q);
2246
  z = 1;
2247
  if (m!=1)
2248
  {
2249
    GEN F = factoru(m);
2250
    long i, j, e;
2251
    ulong r, zeta, y, l;
2252
    for (i = nbrows(F); i; i--)
2253
    {
2254
      l = ucoeff(F,i,1);
2255
      j = ucoeff(F,i,2);
2256
      e = u_lvalrem(q,l, &r);
2257
      y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
2258
      if (zetan)
2259
        z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
2260
      if (a!=1)
2261
        do
2262
        {
2263
          a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
2264
          if (a==ULONG_MAX) return ULONG_MAX;
2265
        } while (--j);
2266
    }
2267
  }
2268
  if (m != nn)
2269
  {
2270
    ulong qm = q/m, nm = nn/m;
2271
    a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
2272
  }
2273
  if (n < 0) a = Fl_inv(a, p);
2274
  if (zetan) *zetan = z;
2275
  return a;
2276
}
2277

2278
ulong
2279
Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
2280
{
2281
  ulong pi = get_Fl_red(p);
2282
  return Fl_sqrtn_pre(a, n, p, pi, zetan);
2283
}
2284

2285
/* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2286
 * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2287
 * and in average is about 2 or 2.5 times worse. But try both algorithms for
2288
 * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2289
 *
2290
 * If X^2 := t^2 - a  is not a square in F_p (so X is in F_p^2), then
2291
 *   (t+X)^(p+1) = (t-X)(t+X) = a,   hence  sqrt(a) = (t+X)^((p+1)/2)  in F_p^2.
2292
 * If (a|p)=1, then sqrt(a) is in F_p.
2293
 * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
2294

2295
/* compute y^2, y = y[1] + y[2] X */
2296
static GEN
2297
sqrt_Cipolla_sqr(void *data, GEN y)
2298
{
2299
  GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2300
  GEN u2 = sqri(u), v2 = sqri(v);
2301
  v = subii(sqri(addii(v,u)), addii(u2,v2));
2302
  u = addii(u2, mulii(v2,n));
2303
  /* NOT mkvec2: must be gerepileupto-able */
2304
  retmkvec2(modii(u,p), modii(v,p));
2305
}
2306
/* compute (t+X) y^2 */
2307
static GEN
2308
sqrt_Cipolla_msqr(void *data, GEN y)
2309
{
2310
  GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2311
  ulong t = gt[2];
2312
  GEN d = addii(u, mului(t,v)), d2= sqri(d);
2313
  GEN b = remii(mulii(a,v), p);
2314
  u = subii(mului(t,d2), mulii(b,addii(u,d)));
2315
  v = subii(d2, mulii(b,v));
2316
  /* NOT mkvec2: must be gerepileupto-able */
2317
  retmkvec2(modii(u,p), modii(v,p));
2318
}
2319
/* assume a reduced mod p [ otherwise correct but inefficient ] */
2320
static GEN
2321
sqrt_Cipolla(GEN a, GEN p)
2322
{
2323
  pari_sp av1;
2324
  GEN u, v, n, y, pov2;
2325
  ulong t;
2326

2327
  if (kronecker(a, p) < 0) return NULL;
2328
  pov2 = shifti(p,-1);
2329
  if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2330

2331
  av1 = avma;
2332
  for(t=1; ; t++)
2333
  {
2334
    n = subsi((long)(t*t), a);
2335
    if (kronecker(n, p) < 0) break;
2336
    set_avma(av1);
2337
  }
2338

2339
  /* compute (t+X)^((p-1)/2) =: u+vX */
2340
  u = utoipos(t);
2341
  y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2342
                         sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2343
  /* Now u+vX = (t+X)^((p-1)/2); thus
2344
   *   (u+vX)(t+X) = sqrt(a) + 0 X
2345
   * Whence,
2346
   *   sqrt(a) = (u+vt)t - v*a
2347
   *   0       = (u+vt)
2348
   * Thus a square root is v*a */
2349

2350
  v = Fp_mul(gel(y, 2), a, p);
2351
  if (cmpii(v,pov2) > 0) v = subii(p,v);
2352
  return v;
2353
}
2354

2355
/* Return NULL if p is found to be composite */
2356
static GEN
2357
Fp_2gener_all(long e, GEN p)
2358
{
2359
  GEN y, m;
2360
  long k;
2361
  GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
2362
  if (e==0 && !equaliu(p,2)) return NULL;
2363
  for (k=2; ; k++)
2364
  {
2365
    long i = krosi(k, p);
2366
    if (i >= 0)
2367
    {
2368
      if (i) continue;
2369
      return NULL;
2370
    }
2371
    y = m = Fp_pow(utoi(k), q, p);
2372
    for (i=1; i<e; i++)
2373
      if (equali1(m = Fp_sqr(m, p))) break;
2374
    if (i == e) break; /* success */
2375
  }
2376
  return y;
2377
}
2378

2379
/* Return NULL if p is found to be composite */
2380
GEN
2381
Fp_2gener(GEN p)
2382
{ return Fp_2gener_all(vali(subis(p,1)),p); }
2383

2384
/* smallest square root */
2385
static GEN
2386
choose_sqrt(GEN v, GEN p)
2387
{
2388
  pari_sp av = avma;
2389
  GEN q = subii(p,v);
2390
  if (cmpii(v,q) > 0) v = q; else set_avma(av);
2391
  return v;
2392
}
2393
/* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2394
GEN
2395
Fp_sqrt_i(GEN a, GEN y, GEN p)
2396
{
2397
  pari_sp av = avma;
2398
  long i, k, e;
2399
  GEN p1, q, v, w;
2400

2401
  if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2402
  if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2403
  if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2404
  if (lgefint(p) == 3)
2405
  {
2406
    ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2407
    if (u == ~0UL) return NULL;
2408
    return utoi(u);
2409
  }
2410

2411
  a = modii(a, p); if (!signe(a)) { set_avma(av); return gen_0; }
2412
  p1 = subiu(p,1); e = vali(p1);
2413
  if (e <= 2)
2414
  { /* direct formulas more efficient */
2415
    pari_sp av2;
2416
    if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
2417
    if (e == 1)
2418
    {
2419
      q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
2420
      v = Fp_pow(a, q, p);
2421
    }
2422
    else
2423
    { /* Atkin's formula */
2424
      GEN i, a2 = shifti(a,1);
2425
      if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
2426
      q = shifti(p1, -3); /* (p-5)/8 */
2427
      v = Fp_pow(a2, q, p);
2428
      i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
2429
      v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
2430
    }
2431
    av2 = avma;
2432
    /* must check equality in case (a/p) = -1 or p not prime */
2433
    e = equalii(Fp_sqr(v,p), a); set_avma(av2);
2434
    return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
2435
  }
2436
  /* On average, Cipolla is better than Tonelli/Shanks if and only if
2437
   * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
2438
  if (e*(e-1) > 20 + 8 * expi(p))
2439
  {
2440
    v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
2441
    return gerepileuptoint(av,v);
2442
  }
2443
  if (!y)
2444
  {
2445
    y = Fp_2gener_all(e, p);
2446
    if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
2447
  }
2448
  q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2449
  p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
2450
  v = Fp_mul(a, p1, p);
2451
  w = Fp_mul(v, p1, p);
2452
  while (!equali1(w))
2453
  { /* a*w = v^2, y primitive 2^e-th root of 1
2454
       a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2455
    p1 = Fp_sqr(w,p);
2456
    for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
2457
    if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
2458
    /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2459
    p1 = y;
2460
    for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
2461
    y = Fp_sqr(p1, p); e = k;
2462
    w = Fp_mul(y, w, p);
2463
    v = Fp_mul(v, p1, p);
2464
    if (gc_needed(av,1))
2465
    {
2466
      if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2467
      gerepileall(av,3, &y,&w,&v);
2468
    }
2469
  }
2470
  return gerepileuptoint(av, choose_sqrt(v,p));
2471
}
2472

2473
GEN
2474
Fp_sqrt(GEN a, GEN p)
2475
{
2476
  return Fp_sqrt_i(a, NULL, p);
2477
}
2478

2479
/*********************************************************************/
2480
/**                                                                 **/
2481
/**                        GCD & BEZOUT                             **/
2482
/**                                                                 **/
2483
/*********************************************************************/
2484

2485
GEN
2486
lcmii(GEN x, GEN y)
2487
{
2488
  pari_sp av;
2489
  GEN a, b;
2490
  if (!signe(x) || !signe(y)) return gen_0;
2491
  av = avma; a = gcdii(x,y);
2492
  if (absequalii(a,y)) { set_avma(av); return absi(x); }
2493
  if (!equali1(a)) y = diviiexact(y,a);
2494
  b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2495
}
2496

2497
/* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
2498
 * set *pd = gcd(x,N) */
2499
GEN
2500
Fp_invgen(GEN x, GEN N, GEN *pd)
2501
{
2502
  GEN d, d0, e, v;
2503
  if (lgefint(N) == 3)
2504
  {
2505
    ulong dd, NN = N[2], xx = umodiu(x,NN);
2506
    if (!xx) { *pd = N; return gen_0; }
2507
    xx = Fl_invgen(xx, NN, &dd);
2508
    *pd = utoi(dd); return utoi(xx);
2509
  }
2510
  *pd = d = bezout(x, N, &v, NULL);
2511
  if (equali1(d)) return v;
2512
  /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
2513
  e = diviiexact(N,d);
2514
  d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
2515
  if (equali1(d0)) return v;
2516
  if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
2517
  return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
2518
}
2519

2520
/*********************************************************************/
2521
/**                                                                 **/
2522
/**                      CHINESE REMAINDERS                         **/
2523
/**                                                                 **/
2524
/*********************************************************************/
2525

2526
/* Chinese Remainder Theorem.  x and y must have the same type (integermod,
2527
 * polymod, or polynomial/vector/matrix recursively constructed with these
2528
 * as coefficients). Creates (with the same type) a z in the same residue
2529
 * class as x and the same residue class as y, if it is possible.
2530
 *
2531
 * We also allow (during recursion) two identical objects even if they are
2532
 * not integermod or polymod. For example:
2533
 *
2534
 * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
2535
 * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
2536
 * ? chinese(x, y)
2537
 * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
2538

2539
static GEN
2540
gen_chinese(GEN x, GEN(*f)(GEN,GEN))
2541
{
2542
  GEN z = gassoc_proto(f,x,NULL);
2543
  if (z == gen_1) retmkintmod(gen_0,gen_1);
2544
  return z;
2545
}
2546

2547
/* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
2548
 * call chinese: makes Mod(0,1) a better "neutral" element */
2549
static GEN
2550
chinese_intpol(GEN x,GEN y)
2551
{
2552
  pari_sp av = avma;
2553
  GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
2554
  return gerepileupto(av, chinese(z, y));
2555
}
2556

2557
GEN
2558
chinese1(GEN x) { return gen_chinese(x,chinese); }
2559

2560
GEN
2561
chinese(GEN x, GEN y)
2562
{
2563
  pari_sp av;
2564
  long tx = typ(x), ty;
2565
  GEN z,p1,p2,d,u,v;
2566

2567
  if (!y) return chinese1(x);
2568
  if (gidentical(x,y)) return gcopy(x);
2569
  ty = typ(y);
2570
  if (tx == ty) switch(tx)
2571
  {
2572
    case t_POLMOD:
2573
    {
2574
      GEN A = gel(x,1), B = gel(y,1);
2575
      GEN a = gel(x,2), b = gel(y,2);
2576
      if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2577
      if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
2578
      av = avma;
2579
      d = RgX_extgcd(A,B,&u,&v);
2580
      p2 = gsub(b, a);
2581
      if (!gequal0(gmod(p2, d))) break;
2582
      p1 = gdiv(A,d);
2583
      p2 = gadd(a, gmul(gmul(u,p1), p2));
2584

2585
      z = cgetg(3, t_POLMOD);
2586
      gel(z,1) = gmul(p1,B);
2587
      gel(z,2) = gmod(p2,gel(z,1));
2588
      return gerepileupto(av, z);
2589
    }
2590
    case t_INTMOD:
2591
    {
2592
      GEN A = gel(x,1), B = gel(y,1);
2593
      GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2594
      z = cgetg(3,t_INTMOD);
2595
      Z_chinese_pre(A, B, &C, &U, &d);
2596
      c = Z_chinese_post(a, b, C, U, d);
2597
      if (!c) pari_err_OP("chinese", x,y);
2598
      set_avma((pari_sp)z);
2599
      gel(z,1) = icopy(C);
2600
      gel(z,2) = icopy(c); return z;
2601
    }
2602
    case t_POL:
2603
    {
2604
      long i, lx = lg(x), ly = lg(y);
2605
      if (varn(x) != varn(y)) break;
2606
      if (lx < ly) { swap(x,y); lswap(lx,ly); }
2607
      z = cgetg(lx, t_POL); z[1] = x[1];
2608
      for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2609
      for (   ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2610
      return z;
2611
    }
2612

2613
    case t_VEC: case t_COL: case t_MAT:
2614
    {
2615
      long i, lx;
2616
      z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2617
      for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2618
      return z;
2619
    }
2620
  }
2621
  if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
2622
  if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
2623
  pari_err_OP("chinese",x,y);
2624
  return NULL; /* LCOV_EXCL_LINE */
2625
}
2626

2627
/* init chinese(Mod(.,A), Mod(.,B)) */
2628
void
2629
Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2630
{
2631
  GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2632
  GEN t = diviiexact(A,d);
2633
  *pU = mulii(u, t);
2634
  *pC = mulii(t, B);
2635
  if (pd) *pd = d;
2636
}
2637
/* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2638
 * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2639
 * If d not NULL, check whether a = b mod d. */
2640
GEN
2641
Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2642
{
2643
  GEN b_a;
2644
  if (!signe(a))
2645
  {
2646
    if (d && !dvdii(b, d)) return NULL;
2647
    return Fp_mul(b, U, C);
2648
  }
2649
  b_a = subii(b,a);
2650
  if (d && !dvdii(b_a, d)) return NULL;
2651
  return modii(addii(a, mulii(U, b_a)), C);
2652
}
2653
static ulong
2654
u_chinese_post(ulong a, ulong b, ulong C, ulong U)
2655
{
2656
  if (!a) return Fl_mul(b, U, C);
2657
  return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
2658
}
2659

2660
GEN
2661
Z_chinese(GEN a, GEN b, GEN A, GEN B)
2662
{
2663
  pari_sp av = avma;
2664
  GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2665
  return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2666
}
2667
GEN
2668
Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2669
{
2670
  GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2671
  return Z_chinese_post(a,b, *pC, U, NULL);
2672
}
2673

2674
/* return lift(chinese(a mod A, b mod B))
2675
 * assume(A,B)=1, a,b,A,B integers and C = A*B */
2676
GEN
2677
Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2678
{
2679
  pari_sp av = avma;
2680
  GEN U = mulii(Fp_inv(A,B), A);
2681
  return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2682
}
2683
ulong
2684
u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
2685
{ return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
2686

2687
/* chinese1 for coprime moduli in Z */
2688
static GEN
2689
chinese1_coprime_Z_aux(GEN x, GEN y)
2690
{
2691
  GEN z = cgetg(3, t_INTMOD);
2692
  GEN A = gel(x,1), a = gel(x, 2);
2693
  GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2694
  pari_sp av = avma;
2695
  GEN U = mulii(Fp_inv(A,B), A);
2696
  gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2697
  gel(z,1) = C; return z;
2698
}
2699
GEN
2700
chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
2701

2702
/*********************************************************************/
2703
/**                                                                 **/
2704
/**                    MODULAR EXPONENTIATION                       **/
2705
/**                                                                 **/
2706
/*********************************************************************/
2707

2708
/* xa, ya = t_VECSMALL */
2709
GEN
2710
ZV_producttree(GEN xa)
2711
{
2712
  long n = lg(xa)-1;
2713
  long m = n==1 ? 1: expu(n-1)+1;
2714
  GEN T = cgetg(m+1, t_VEC), t;
2715
  long i, j, k;
2716
  t = cgetg(((n+1)>>1)+1, t_VEC);
2717
  if (typ(xa)==t_VECSMALL)
2718
  {
2719
    for (j=1, k=1; k<n; j++, k+=2)
2720
      gel(t, j) = muluu(xa[k], xa[k+1]);
2721
    if (k==n) gel(t, j) = utoi(xa[k]);
2722
  } else {
2723
    for (j=1, k=1; k<n; j++, k+=2)
2724
      gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2725
    if (k==n) gel(t, j) = icopy(gel(xa,k));
2726
  }
2727
  gel(T,1) = t;
2728
  for (i=2; i<=m; i++)
2729
  {
2730
    GEN u = gel(T, i-1);
2731
    long n = lg(u)-1;
2732
    t = cgetg(((n+1)>>1)+1, t_VEC);
2733
    for (j=1, k=1; k<n; j++, k+=2)
2734
      gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2735
    if (k==n) gel(t, j) = gel(u, k);
2736
    gel(T, i) = t;
2737
  }
2738
  return T;
2739
}
2740

2741
/* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
2742
GEN
2743
Z_ZV_mod_tree(GEN A, GEN P, GEN T)
2744
{
2745
  long i,j,k;
2746
  long m = lg(T)-1, n = lg(P)-1;
2747
  GEN t;
2748
  GEN Tp = cgetg(m+1, t_VEC);
2749
  gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
2750
  for (i=m-1; i>=1; i--)
2751
  {
2752
    GEN u = gel(T, i);
2753
    GEN v = gel(Tp, i+1);
2754
    long n = lg(u)-1;
2755
    t = cgetg(n+1, t_VEC);
2756
    for (j=1, k=1; k<n; j++, k+=2)
2757
    {
2758
      gel(t, k)   = modii(gel(v, j), gel(u, k));
2759
      gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2760
    }
2761
    if (k==n) gel(t, k) = gel(v, j);
2762
    gel(Tp, i) = t;
2763
  }
2764
  {
2765
    GEN u = gel(T, i+1);
2766
    GEN v = gel(Tp, i+1);
2767
    long l = lg(u)-1;
2768
    if (typ(P)==t_VECSMALL)
2769
    {
2770
      GEN R = cgetg(n+1, t_VECSMALL);
2771
      for (j=1, k=1; j<=l; j++, k+=2)
2772
      {
2773
        uel(R,k) = umodiu(gel(v, j), P[k]);
2774
        if (k < n)
2775
          uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
2776
      }
2777
      return R;
2778
    }
2779
    else
2780
    {
2781
      GEN R = cgetg(n+1, t_VEC);
2782
      for (j=1, k=1; j<=l; j++, k+=2)
2783
      {
2784
        gel(R,k) = modii(gel(v, j), gel(P,k));
2785
        if (k < n)
2786
          gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
2787
      }
2788
      return R;
2789
    }
2790
  }
2791
}
2792

2793
/* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
2794
GEN
2795
ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
2796
{
2797
  long m = lg(T)-1, n = lg(A)-1;
2798
  long i,j,k;
2799
  GEN Tp = cgetg(m+1, t_VEC);
2800
  GEN M = gel(T, 1);
2801
  GEN t = cgetg(lg(M), t_VEC);
2802
  if (typ(P)==t_VECSMALL)
2803
  {
2804
    for (j=1, k=1; k<n; j++, k+=2)
2805
    {
2806
      pari_sp av = avma;
2807
      GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
2808
      GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
2809
      gel(t, j) = gerepileuptoint(av, tj);
2810
    }
2811
    if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
2812
  } else
2813
  {
2814
    for (j=1, k=1; k<n; j++, k+=2)
2815
    {
2816
      pari_sp av = avma;
2817
      GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
2818
      GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
2819
      gel(t, j) = gerepileuptoint(av, tj);
2820
    }
2821
    if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
2822
  }
2823
  gel(Tp, 1) = t;
2824
  for (i=2; i<=m; i++)
2825
  {
2826
    GEN u = gel(T, i-1), M = gel(T, i);
2827
    GEN t = cgetg(lg(M), t_VEC);
2828
    GEN v = gel(Tp, i-1);
2829
    long n = lg(v)-1;
2830
    for (j=1, k=1; k<n; j++, k+=2)
2831
    {
2832
      pari_sp av = avma;
2833
      gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2834
            mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2835
    }
2836
    if (k==n) gel(t, j) = gel(v, k);
2837
    gel(Tp, i) = t;
2838
  }
2839
  return gmael(Tp,m,1);
2840
}
2841

2842
static GEN
2843
ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2844
{
2845
  long i, l = lg(gel(vA,1)), n = lg(P);
2846
  GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
2847
  for (i=1; i < l; i++)
2848
  {
2849
    pari_sp av = avma;
2850
    GEN c, A = cgetg(n, typ(P));
2851
    long j;
2852
    for (j=1; j < n; j++) A[j] = mael(vA,j,i);
2853
    c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2854
    gel(V,i) = gerepileuptoint(av, c);
2855
  }
2856
  return V;
2857
}
2858

2859
static GEN
2860
nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2861
{
2862
  long i, j, l, n = lg(P);
2863
  GEN mod = gmael(T, lg(T)-1, 1), V, w;
2864
  w = cgetg(n, t_VECSMALL);
2865
  for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
2866
  l = vecsmall_max(w);
2867
  V = cgetg(l, t_POL);
2868
  V[1] = mael(vA,1,1);
2869
  for (i=2; i < l; i++)
2870
  {
2871
    pari_sp av = avma;
2872
    GEN c, A = cgetg(n, typ(P));
2873
    if (typ(P)==t_VECSMALL)
2874
      for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
2875
    else
2876
      for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
2877
    c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2878
    gel(V,i) = gerepileuptoint(av, c);
2879
  }
2880
  return ZX_renormalize(V, l);
2881
}
2882

2883
static GEN
2884
nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2885
{
2886
  long i, j, l = lg(gel(vA,1)), n = lg(P);
2887
  GEN A = cgetg(n, t_VEC);
2888
  GEN V = cgetg(l, t_COL);
2889
  for (i=1; i < l; i++)
2890
  {
2891
    for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2892
    gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
2893
  }
2894
  return V;
2895
}
2896

2897
static GEN
2898
polint_chinese(GEN worker, GEN mA, GEN P)
2899
{
2900
  long cnt, pending, n, i, j, l = lg(gel(mA,1));
2901
  struct pari_mt pt;
2902
  GEN done, va, M, A;
2903
  pari_timer ti;
2904

2905
  if (l == 1) return cgetg(1, t_MAT);
2906
  cnt = pending = 0;
2907
  n = lg(P);
2908
  A = cgetg(n, t_VEC);
2909
  va = mkvec(A);
2910
  M = cgetg(l, t_MAT);
2911
  if (DEBUGLEVEL>4) timer_start(&ti);
2912
  if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
2913
  mt_queue_start_lim(&pt, worker, l-1);
2914
  for (i=1; i<l || pending; i++)
2915
  {
2916
    long workid;
2917
    for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
2918
    mt_queue_submit(&pt, i, i<l? va: NULL);
2919
    done = mt_queue_get(&pt, &workid, &pending);
2920
    if (done)
2921
    {
2922
      gel(M,workid) = done;
2923
      if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
2924
    }
2925
  }
2926
  if (DEBUGLEVEL>5) err_printf("\n");
2927
  if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
2928
  mt_queue_end(&pt);
2929
  return M;
2930
}
2931

2932
GEN
2933
nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2934
{
2935
  return nxCV_polint_center_tree(vA, P, T, R, m2);
2936
}
2937

2938
static GEN
2939
nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2940
{
2941
  long i, j, l = lg(gel(vA,1)), n = lg(P);
2942
  GEN A = cgetg(n, t_VEC);
2943
  GEN V = cgetg(l, t_MAT);
2944
  for (i=1; i < l; i++)
2945
  {
2946
    for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2947
    gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
2948
  }
2949
  return V;
2950
}
2951

2952
static GEN
2953
nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2954
{
2955
  GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
2956
  return polint_chinese(worker, mA, P);
2957
}
2958

2959
static GEN
2960
nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2961
{
2962
  long i, j, l = lg(gel(vA,1)), n = lg(P);
2963
  GEN A = cgetg(n, t_VEC);
2964
  GEN V = cgetg(l, t_MAT);
2965
  for (i=1; i < l; i++)
2966
  {
2967
    for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2968
    gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
2969
  }
2970
  return V;
2971
}
2972

2973
GEN
2974
nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2975
{
2976
  return ncV_polint_center_tree(vA, P, T, R, m2);
2977
}
2978

2979
static GEN
2980
nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2981
{
2982
  GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
2983
  return polint_chinese(worker, mA, P);
2984
}
2985

2986
/* return [A mod P[i], i=1..#P] */
2987
GEN
2988
Z_ZV_mod(GEN A, GEN P)
2989
{
2990
  pari_sp av = avma;
2991
  return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2992
}
2993
/* P a t_VECSMALL */
2994
GEN
2995
Z_nv_mod(GEN A, GEN P)
2996
{
2997
  pari_sp av = avma;
2998
  return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2999
}
3000
/* B a ZX, T = ZV_producttree(P) */
3001
GEN
3002
ZX_nv_mod_tree(GEN B, GEN A, GEN T)
3003
{
3004
  pari_sp av;
3005
  long i, j, l = lg(B), n = lg(A)-1;
3006
  GEN V = cgetg(n+1, t_VEC);
3007
  for (j=1; j <= n; j++)
3008
  {
3009
    gel(V, j) = cgetg(l, t_VECSMALL);
3010
    mael(V, j, 1) = B[1]&VARNBITS;
3011
  }
3012
  av = avma;
3013
  for (i=2; i < l; i++)
3014
  {
3015
    GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
3016
    for (j=1; j <= n; j++)
3017
      mael(V, j, i) = v[j];
3018
    set_avma(av);
3019
  }
3020
  for (j=1; j <= n; j++)
3021
    (void) Flx_renormalize(gel(V, j), l);
3022
  return V;
3023
}
3024

3025
static GEN
3026
to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
3027

3028
GEN
3029
ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
3030
{
3031
  pari_sp av = avma;
3032
  long i, j, l = lg(P), n = lg(xa)-1;
3033
  GEN V = cgetg(n+1, t_VEC);
3034
  for (j=1; j <= n; j++)
3035
  {
3036
    gel(V, j) = cgetg(l, t_POL);
3037
    mael(V, j, 1) = P[1]&VARNBITS;
3038
  }
3039
  for (i=2; i < l; i++)
3040
  {
3041
    GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
3042
    for (j=1; j <= n; j++)
3043
      gmael(V, j, i) = gel(v,j);
3044
  }
3045
  for (j=1; j <= n; j++)
3046
    (void) FlxX_renormalize(gel(V, j), l);
3047
  return gerepilecopy(av, V);
3048
}
3049

3050
GEN
3051
ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
3052
{
3053
  pari_sp av = avma;
3054
  long i, j, l = lg(C), n = lg(xa)-1;
3055
  GEN V = cgetg(n+1, t_VEC);
3056
  for (j = 1; j <= n; j++)
3057
    gel(V, j) = cgetg(l, t_COL);
3058
  for (i = 1; i < l; i++)
3059
  {
3060
    GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
3061
    for (j = 1; j <= n; j++)
3062
      gmael(V, j, i) = gel(v,j);
3063
  }
3064
  return gerepilecopy(av, V);
3065
}
3066

3067
GEN
3068
ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
3069
{
3070
  pari_sp av = avma;
3071
  long i, j, l = lg(M), n = lg(xa)-1;
3072
  GEN V = cgetg(n+1, t_VEC);
3073
  for (j=1; j <= n; j++)
3074
    gel(V, j) = cgetg(l, t_MAT);
3075
  for (i=1; i < l; i++)
3076
  {
3077
    GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
3078
    for (j=1; j <= n; j++)
3079
      gmael(V, j, i) = gel(v,j);
3080
  }
3081
  return gerepilecopy(av, V);
3082
}
3083

3084
GEN
3085
ZV_nv_mod_tree(GEN B, GEN A, GEN T)
3086
{
3087
  pari_sp av;
3088
  long i, j, l = lg(B), n = lg(A)-1;
3089
  GEN V = cgetg(n+1, t_VEC);
3090
  for (j=1; j <= n; j++)
3091
    gel(V, j) = cgetg(l, t_VECSMALL);
3092
  av = avma;
3093
  for (i=1; i < l; i++)
3094
  {
3095
    GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
3096
    for (j=1; j <= n; j++)
3097
      mael(V, j, i) = v[j];
3098
    set_avma(av);
3099
  }
3100
  return V;
3101
}
3102

3103
GEN
3104
ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
3105
{
3106
  pari_sp av = avma;
3107
  long i, j, l = lg(M), n = lg(xa)-1;
3108
  GEN V = cgetg(n+1, t_VEC);
3109
  for (j=1; j <= n; j++)
3110
    gel(V, j) = cgetg(l, t_MAT);
3111
  for (i=1; i < l; i++)
3112
  {
3113
    GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
3114
    for (j=1; j <= n; j++)
3115
      gmael(V, j, i) = gel(v,j);
3116
  }
3117
  return gerepilecopy(av, V);
3118
}
3119

3120
static GEN
3121
ZV_sqr(GEN z)
3122
{
3123
  long i,l = lg(z);
3124
  GEN x = cgetg(l, t_VEC);
3125
  if (typ(z)==t_VECSMALL)
3126
    for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
3127
  else
3128
    for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
3129
  return x;
3130
}
3131

3132
static GEN
3133
ZT_sqr(GEN x)
3134
{
3135
  if (typ(x) == t_INT)
3136
    return sqri(x);
3137
  pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
3138
}
3139

3140
static GEN
3141
ZV_invdivexact(GEN y, GEN x)
3142
{
3143
  long i, l = lg(y);
3144
  GEN z = cgetg(l,t_VEC);
3145
  if (typ(x)==t_VECSMALL)
3146
    for (i=1; i<l; i++)
3147
    {
3148
      pari_sp av = avma;
3149
      ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
3150
      set_avma(av);
3151
      gel(z,i) = utoi(a);
3152
    }
3153
  else
3154
    for (i=1; i<l; i++)
3155
      gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
3156
  return z;
3157
}
3158

3159
/* P t_VECSMALL or t_VEC of t_INT  */
3160
GEN
3161
ZV_chinesetree(GEN P, GEN T)
3162
{
3163
  GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
3164
  GEN mod = gmael(T,lg(T)-1,1);
3165
  return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
3166
}
3167

3168
static GEN
3169
gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
3170
{
3171
  if (!pt_mod)
3172
    return gerepileupto(av, a);
3173
  else
3174
  {
3175
    GEN mod = gmael(T, lg(T)-1, 1);
3176
    gerepileall(av, 2, &a, &mod);
3177
    *pt_mod = mod;
3178
    return a;
3179
  }
3180
}
3181

3182
GEN
3183
ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3184
{
3185
  pari_sp av = avma;
3186
  GEN T = ZV_producttree(P);
3187
  GEN R = ZV_chinesetree(P, T);
3188
  GEN a = ZV_chinese_tree(A, P, T, R);
3189
  GEN mod = gmael(T, lg(T)-1, 1);
3190
  GEN ca = Fp_center(a, mod, shifti(mod,-1));
3191
  return gc_chinese(av, T, ca, pt_mod);
3192
}
3193

3194
GEN
3195
ZV_chinese(GEN A, GEN P, GEN *pt_mod)
3196
{
3197
  pari_sp av = avma;
3198
  GEN T = ZV_producttree(P);
3199
  GEN R = ZV_chinesetree(P, T);
3200
  GEN a = ZV_chinese_tree(A, P, T, R);
3201
  return gc_chinese(av, T, a, pt_mod);
3202
}
3203

3204
GEN
3205
nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3206
{
3207
  pari_sp av = avma;
3208
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3209
  GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3210
  return gerepileupto(av, a);
3211
}
3212

3213
GEN
3214
nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3215
{
3216
  pari_sp av = avma;
3217
  GEN T = ZV_producttree(P);
3218
  GEN R = ZV_chinesetree(P, T);
3219
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3220
  GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3221
  return gc_chinese(av, T, a, pt_mod);
3222
}
3223

3224
GEN
3225
ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3226
{
3227
  pari_sp av = avma;
3228
  GEN T = ZV_producttree(P);
3229
  GEN R = ZV_chinesetree(P, T);
3230
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3231
  GEN a = ncV_polint_center_tree(A, P, T, R, m2);
3232
  return gc_chinese(av, T, a, pt_mod);
3233
}
3234

3235
GEN
3236
ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3237
{
3238
  pari_sp av = avma;
3239
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3240
  GEN a = ncV_polint_center_tree(A, P, T, R, m2);
3241
  return gerepileupto(av, a);
3242
}
3243

3244
GEN
3245
nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3246
{
3247
  pari_sp av = avma;
3248
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3249
  GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3250
  return gerepileupto(av, a);
3251
}
3252

3253
GEN
3254
nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3255
{
3256
  pari_sp av = avma;
3257
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3258
  GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
3259
  return gerepileupto(av, a);
3260
}
3261

3262
GEN
3263
nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3264
{
3265
  pari_sp av = avma;
3266
  GEN T = ZV_producttree(P);
3267
  GEN R = ZV_chinesetree(P, T);
3268
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3269
  GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3270
  return gc_chinese(av, T, a, pt_mod);
3271
}
3272

3273
GEN
3274
nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3275
{
3276
  pari_sp av = avma;
3277
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3278
  GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
3279
  return gerepileupto(av, a);
3280
}
3281

3282
GEN
3283
nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3284
{
3285
  pari_sp av = avma;
3286
  GEN T = ZV_producttree(P);
3287
  GEN R = ZV_chinesetree(P, T);
3288
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3289
  GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
3290
  return gc_chinese(av, T, a, pt_mod);
3291
}
3292

3293
GEN
3294
nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3295
{
3296
  pari_sp av = avma;
3297
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3298
  GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
3299
  return gerepileupto(av, a);
3300
}
3301

3302
GEN
3303
nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3304
{
3305
  pari_sp av = avma;
3306
  GEN T = ZV_producttree(P);
3307
  GEN R = ZV_chinesetree(P, T);
3308
  GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3309
  GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
3310
  return gc_chinese(av, T, a, pt_mod);
3311
}
3312

3313
/**********************************************************************
3314
 **                                                                  **
3315
 **                    Powering  over (Z/NZ)^*, small N              **
3316
 **                                                                  **
3317
 **********************************************************************/
3318

3319
/* 2^n mod p; assume n > 1 */
3320
static ulong
3321
Fl_2powu_pre(ulong n, ulong p, ulong pi)
3322
{
3323
  ulong y = 2;
3324
  int j = 1+bfffo(n);
3325
  /* normalize, i.e set highest bit to 1 (we know n != 0) */
3326
  n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3327
  for (; j; n<<=1,j--)
3328
  {
3329
    y = Fl_sqr_pre(y,p,pi);
3330
    if (n & HIGHBIT) y = Fl_double(y, p);
3331
  }
3332
  return y;
3333
}
3334

3335
/* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
3336
static ulong
3337
Fl_2powu(ulong n, ulong p)
3338
{
3339
  ulong y = 2;
3340
  int j = 1+bfffo(n);
3341
  /* normalize, i.e set highest bit to 1 (we know n != 0) */
3342
  n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3343
  for (; j; n<<=1,j--)
3344
  {
3345
    y = (y*y) % p;
3346
    if (n & HIGHBIT) y = Fl_double(y, p);
3347
  }
3348
  return y;
3349
}
3350

3351
ulong
3352
Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
3353
{
3354
  ulong y, z, n;
3355
  if (n0 <= 1)
3356
  { /* frequent special cases */
3357
    if (n0 == 1) return x;
3358
    if (n0 == 0) return 1;
3359
  }
3360
  if (x <= 2)
3361
  {
3362
    if (x == 2) return Fl_2powu_pre(n0, p, pi);
3363
    return x; /* 0 or 1 */
3364
  }
3365
  y = 1; z = x; n = n0;
3366
  for(;;)
3367
  {
3368
    if (n&1) y = Fl_mul_pre(y,z,p,pi);
3369
    n>>=1; if (!n) return y;
3370
    z = Fl_sqr_pre(z,p,pi);
3371
  }
3372
}
3373

3374
ulong
3375
Fl_powu(ulong x, ulong n0, ulong p)
3376
{
3377
  ulong y, z, n;
3378
  if (n0 <= 2)
3379
  { /* frequent special cases */
3380
    if (n0 == 2) return Fl_sqr(x,p);
3381
    if (n0 == 1) return x;
3382
    if (n0 == 0) return 1;
3383
  }
3384
  if (x <= 1) return x; /* 0 or 1 */
3385
  if (p & HIGHMASK)
3386
    return Fl_powu_pre(x, n0, p, get_Fl_red(p));
3387
  if (x == 2) return Fl_2powu(n0, p);
3388
  y = 1; z = x; n = n0;
3389
  for(;;)
3390
  {
3391
    if (n&1) y = (y*z) % p;
3392
    n>>=1; if (!n) return y;
3393
    z = (z*z) % p;
3394
  }
3395
}
3396

3397
/* Reduce data dependency to maximize internal parallelism */
3398
GEN
3399
Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
3400
{
3401
  long i, k;
3402
  GEN powers = cgetg(n + 2, t_VECSMALL);
3403
  powers[1] = 1; if (n == 0) return powers;
3404
  powers[2] = x;
3405
  for (i = 3, k=2; i <= n; i+=2, k++)
3406
  {
3407
    powers[i] = Fl_sqr_pre(powers[k], p, pi);
3408
    powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
3409
  }
3410
  if (i==n+1)
3411
    powers[i] = Fl_sqr_pre(powers[k], p, pi);
3412
  return powers;
3413
}
3414

3415
GEN
3416
Fl_powers(ulong x, long n, ulong p)
3417
{
3418
  return Fl_powers_pre(x, n, p, get_Fl_red(p));
3419
}
3420

3421
/**********************************************************************
3422
 **                                                                  **
3423
 **                    Powering  over (Z/NZ)^*, large N              **
3424
 **                                                                  **
3425
 **********************************************************************/
3426

3427
static GEN
3428
Fp_dblsqr(GEN x, GEN N)
3429
{
3430
  GEN z = shifti(Fp_sqr(x, N), 1);
3431
  return cmpii(z, N) >= 0? subii(z, N): z;
3432
}
3433

3434
typedef struct muldata {
3435
  GEN (*sqr)(void * E, GEN x);
3436
  GEN (*mul)(void * E, GEN x, GEN y);
3437
  GEN (*mul2)(void * E, GEN x);
3438
} muldata;
3439

3440
/* modified Barrett reduction with one fold */
3441
/* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
3442

3443
static GEN
3444
Fp_invmBarrett(GEN p, long s)
3445
{
3446
  GEN R, Q = dvmdii(int2n(3*s),p,&R);
3447
  return mkvec2(Q,R);
3448
}
3449

3450
/* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
3451
 * a = r (mod N) */
3452
static GEN
3453
Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
3454
{
3455
  pari_sp av = avma;
3456
  GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
3457
  long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
3458
  GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
3459
  GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
3460
  GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
3461
  GEN r = subii(A, mulii(q, N));
3462
  GEN sr= subii(r,N);     /* 0 <= r < 4*N */
3463
  if (signe(sr)<0) return gerepileuptoint(av, r);
3464
  r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
3465
  if (signe(sr)<0) return gerepileuptoint(av, r);
3466
  r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
3467
  return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
3468
}
3469

3470
/* Montgomery reduction */
3471

3472
INLINE ulong
3473
init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
3474

3475
struct montred
3476
{
3477
  GEN N;
3478
  ulong inv;
3479
};
3480

3481
/* Montgomery reduction */
3482
static GEN
3483
_sqr_montred(void * E, GEN x)
3484
{
3485
  struct montred * D = (struct montred *) E;
3486
  return red_montgomery(sqri(x), D->N, D->inv);
3487
}
3488

3489
/* Montgomery reduction */
3490
static GEN
3491
_mul_montred(void * E, GEN x, GEN y)
3492
{
3493
  struct montred * D = (struct montred *) E;
3494
  return red_montgomery(mulii(x, y), D->N, D->inv);
3495
}
3496

3497
static GEN
3498
_mul2_montred(void * E, GEN x)
3499
{
3500
  struct montred * D = (struct montred *) E;
3501
  GEN z = shifti(_sqr_montred(E, x), 1);
3502
  long l = lgefint(D->N);
3503
  while (lgefint(z) > l) z = subii(z, D->N);
3504
  return z;
3505
}
3506

3507
static GEN
3508
_sqr_remii(void* N, GEN x)
3509
{ return remii(sqri(x), (GEN) N); }
3510

3511
static GEN
3512
_mul_remii(void* N, GEN x, GEN y)
3513
{ return remii(mulii(x, y), (GEN) N); }
3514

3515
static GEN
3516
_mul2_remii(void* N, GEN x)
3517
{ return Fp_dblsqr(x, (GEN) N); }
3518

3519
struct redbarrett
3520
{
3521
  GEN iM, N;
3522
  long s;
3523
};
3524

3525
static GEN
3526
_sqr_remiibar(void *E, GEN x)
3527
{
3528
  struct redbarrett * D = (struct redbarrett *) E;
3529
  return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
3530
}
3531

3532
static GEN
3533
_mul_remiibar(void *E, GEN x, GEN y)
3534
{
3535
  struct redbarrett * D = (struct redbarrett *) E;
3536
  return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
3537
}
3538

3539
static GEN
3540
_mul2_remiibar(void *E, GEN x)
3541
{
3542
  struct redbarrett * D = (struct redbarrett *) E;
3543
  return Fp_dblsqr(x, D->N);
3544
}
3545

3546
static long
3547
Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
3548
{
3549
  if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
3550
  {
3551
    struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
3552
    D->sqr = &_sqr_remiibar;
3553
    D->mul = &_mul_remiibar;
3554
    D->mul2 = &_mul2_remiibar;
3555
    E->N = N;
3556
    E->s = 1+(expi(N)>>1);
3557
    E->iM = Fp_invmBarrett(N, E->s);
3558
    *pt_E = (void*) E;
3559
    return 0;
3560
  }
3561
  else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
3562
  {
3563
    struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
3564
    *y = remii(shifti(*y, bit_accuracy(lN)), N);
3565
    D->sqr = &_sqr_montred;
3566
    D->mul = &_mul_montred;
3567
    D->mul2 = &_mul2_montred;
3568
    E->N = N;
3569
    E->inv = init_montdata(N);
3570
    *pt_E = (void*) E;
3571
    return 1;
3572
  }
3573
  else
3574
  {
3575
    D->sqr = &_sqr_remii;
3576
    D->mul = &_mul_remii;
3577
    D->mul2 = &_mul2_remii;
3578
    *pt_E = (void*) N;
3579
    return 0;
3580
  }
3581
}
3582

3583
GEN
3584
Fp_powu(GEN A, ulong k, GEN N)
3585
{
3586
  long lN = lgefint(N);
3587
  int base_is_2, use_montgomery;
3588
  muldata D;
3589
  void *E;
3590
  pari_sp av;
3591

3592
  if (lN == 3) {
3593
    ulong n = uel(N,2);
3594
    return utoi( Fl_powu(umodiu(A, n), k, n) );
3595
  }
3596
  if (k <= 2)
3597
  { /* frequent special cases */
3598
    if (k == 2) return Fp_sqr(A,N);
3599
    if (k == 1) return A;
3600
    if (k == 0) return gen_1;
3601
  }
3602
  av = avma; A = modii(A,N);
3603
  base_is_2 = 0;
3604
  if (lgefint(A) == 3) switch(A[2])
3605
  {
3606
    case 1: set_avma(av); return gen_1;
3607
    case 2:  base_is_2 = 1; break;
3608
  }
3609

3610
  /* TODO: Move this out of here and use for general modular computations */
3611
  use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
3612
  if (base_is_2)
3613
    A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
3614
  else
3615
    A = gen_powu_i(A, k, E, D.sqr, D.mul);
3616
  if (use_montgomery)
3617
  {
3618
    A = red_montgomery(A, N, ((struct montred *) E)->inv);
3619
    if (cmpii(A, N) >= 0) A = subii(A,N);
3620
  }
3621
  return gerepileuptoint(av, A);
3622
}
3623

3624
GEN
3625
Fp_pows(GEN A, long k, GEN N)
3626
{
3627
  if (lgefint(N) == 3) {
3628
    ulong n = N[2];
3629
    ulong a = umodiu(A, n);
3630
    if (k < 0) {
3631
      a = Fl_inv(a, n);
3632
      k = -k;
3633
    }
3634
    return utoi( Fl_powu(a, (ulong)k, n) );
3635
  }
3636
  if (k < 0) { A = Fp_inv(A, N); k = -k; };
3637
  return Fp_powu(A, (ulong)k, N);
3638
}
3639

3640
/* A^K mod N */
3641
GEN
3642
Fp_pow(GEN A, GEN K, GEN N)
3643
{
3644
  pari_sp av;
3645
  long s, lN = lgefint(N), sA, sy;
3646
  int base_is_2, use_montgomery;
3647
  GEN y;
3648
  muldata D;
3649
  void *E;
3650

3651
  s = signe(K);
3652
  if (!s) return dvdii(A,N)? gen_0: gen_1;
3653
  if (lN == 3 && lgefint(K) == 3)
3654
  {
3655
    ulong n = N[2], a = umodiu(A, n);
3656
    if (s < 0) a = Fl_inv(a, n);
3657
    if (a <= 1) return utoi(a); /* 0 or 1 */
3658
    return utoi(Fl_powu(a, uel(K,2), n));
3659
  }
3660

3661
  av = avma;
3662
  if (s < 0) y = Fp_inv(A,N);
3663
  else
3664
  {
3665
    y = modii(A,N);
3666
    if (!signe(y)) { set_avma(av); return gen_0; }
3667
  }
3668
  if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3669

3670
  base_is_2 = 0;
3671
  sy = abscmpii(y, shifti(N,-1)) > 0;
3672
  if (sy) y = subii(N,y);
3673
  sA = sy && mod2(K);
3674
  if (lgefint(y) == 3) switch(y[2])
3675
  {
3676
    case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
3677
    case 2:  base_is_2 = 1; break;
3678
  }
3679

3680
  /* TODO: Move this out of here and use for general modular computations */
3681
  use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
3682
  if (base_is_2)
3683
    y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
3684
  else
3685
    y = gen_pow_i(y, K, E, D.sqr, D.mul);
3686
  if (use_montgomery)
3687
  {
3688
    y = red_montgomery(y, N, ((struct montred *) E)->inv);
3689
    if (cmpii(y,N) >= 0) y = subii(y,N);
3690
  }
3691
  if (sA) y = subii(N, y);
3692
  return gerepileuptoint(av,y);
3693
}
3694

3695
static GEN
3696
_Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3697

3698
static GEN
3699
_Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3700

3701
static GEN
3702
_Fp_one(void *E) { (void) E; return gen_1; }
3703

3704
GEN
3705
Fp_pow_init(GEN x, GEN n, long k, GEN p)
3706
{
3707
  return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul);
3708
}
3709

3710
GEN
3711
Fp_pow_table(GEN R, GEN n, GEN p)
3712
{
3713
  return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul);
3714
}
3715

3716
GEN
3717
Fp_powers(GEN x, long n, GEN p)
3718
{
3719
  if (lgefint(p) == 3)
3720
    return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3721
  return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3722
}
3723

3724
GEN
3725
FpV_prod(GEN V, GEN p)
3726
{
3727
  return gen_product(V, (void *)p, &_Fp_mul);
3728
}
3729

3730
static GEN
3731
_Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3732

3733
static GEN
3734
_Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
3735

3736
static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3737

3738
static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3739
                                      equalii,equali1,Fp_easylog};
3740

3741
static GEN
3742
_Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3743

3744
static GEN
3745
_Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3746

3747
static GEN
3748
_Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3749

3750
static GEN
3751
_Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3752

3753
static GEN
3754
_Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3755

3756
static int
3757
_Fp_equal0(GEN x) { return signe(x)==0; }
3758

3759
static GEN
3760
_Fp_s(void *E, long x) { (void) E; return stoi(x); }
3761

3762
static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3763
                                        _Fp_inv,_Fp_equal0,_Fp_s};
3764

3765
const struct bb_field *get_Fp_field(void **E, GEN p)
3766
{
3767
  *E = (void*)p; return &Fp_field;
3768
}
3769

3770
/*********************************************************************/
3771
/**                                                                 **/
3772
/**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
3773
/**                                                                 **/
3774
/*********************************************************************/
3775
ulong
3776
Fl_order(ulong a, ulong o, ulong p)
3777
{
3778
  pari_sp av = avma;
3779
  GEN m, P, E;
3780
  long i;
3781
  if (a==1) return 1;
3782
  if (!o) o = p-1;
3783
  m = factoru(o);
3784
  P = gel(m,1);
3785
  E = gel(m,2);
3786
  for (i = lg(P)-1; i; i--)
3787
  {
3788
    ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3789
    if (y == 1) o = t;
3790
    else for (j = 1; j < e; j++)
3791
    {
3792
      y = Fl_powu(y, l, p);
3793
      if (y == 1) { o = t *  upowuu(l, j); break; }
3794
    }
3795
  }
3796
  return gc_ulong(av, o);
3797
}
3798

3799
/*Find the exact order of a assuming a^o==1*/
3800
GEN
3801
Fp_order(GEN a, GEN o, GEN p) {
3802
  if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
3803
  {
3804
    ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
3805
    return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3806
  }
3807
  return gen_order(a, o, (void*)p, &Fp_star);
3808
}
3809
GEN
3810
Fp_factored_order(GEN a, GEN o, GEN p)
3811
{ return gen_factored_order(a, o, (void*)p, &Fp_star); }
3812

3813
/* return order of a mod p^e, e > 0, pe = p^e */
3814
static GEN
3815
Zp_order(GEN a, GEN p, long e, GEN pe)
3816
{
3817
  GEN ap, op;
3818
  if (absequaliu(p, 2))
3819
  {
3820
    if (e == 1) return gen_1;
3821
    if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3822
    if (mod4(a) == 1)
3823
      op = gen_1;
3824
    else {
3825
      op = gen_2;
3826
      a = Fp_sqr(a, pe);
3827
    }
3828
  } else {
3829
    ap = (e == 1)? a: remii(a,p);
3830
    op = Fp_order(ap, subiu(p,1), p);
3831
    if (e == 1) return op;
3832
    a = Fp_pow(a, op, pe); /* 1 mod p */
3833
  }
3834
  if (equali1(a)) return op;
3835
  return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
3836
}
3837

3838
GEN
3839
znorder(GEN x, GEN o)
3840
{
3841
  pari_sp av = avma;
3842
  GEN b, a;
3843

3844
  if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3845
  b = gel(x,1); a = gel(x,2);
3846
  if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3847
  if (!o)
3848
  {
3849
    GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3850
    long i, l = lg(P);
3851
    o = gen_1;
3852
    for (i = 1; i < l; i++)
3853
    {
3854
      GEN p = gel(P,i);
3855
      long e = itos(gel(E,i));
3856

3857
      if (l == 2)
3858
        o = Zp_order(a, p, e, b);
3859
      else {
3860
        GEN pe = powiu(p,e);
3861
        o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3862
      }
3863
    }
3864
    return gerepileuptoint(av, o);
3865
  }
3866
  return Fp_order(a, o, b);
3867
}
3868
GEN
3869
order(GEN x) { return znorder(x, NULL); }
3870

3871
/*********************************************************************/
3872
/**                                                                 **/
3873
/**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
3874
/**                                                                 **/
3875
/*********************************************************************/
3876
static GEN
3877
Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3878
{
3879
  pari_sp av = avma;
3880
  GEN h1, h2, F, G;
3881
  if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
3882
  if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3883
  {
3884
    GEN M = cgetg(3, t_MAT);
3885
    gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3886
    gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3887
    return gerepileupto(av, M);
3888
  }
3889
  return gc_NULL(av);
3890
}
3891

3892
static GEN
3893
Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3894
{
3895
  GEN rel;
3896
  do
3897
  {
3898
    (*e)++; *g = Fp_mul(*g, b, p);
3899
    rel = Fp_log_halfgcd(bnd, C, *g, p);
3900
  } while (!rel);
3901
  return rel;
3902
}
3903

3904
struct Fp_log_rel
3905
{
3906
  GEN rel;
3907
  ulong prmax;
3908
  long nbrel, nbmax, nbgen;
3909
};
3910

3911
/* add u^e */
3912
static void
3913
addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
3914
{
3915
  pari_sp av = avma;
3916
  long off = r->prmax+1;
3917
  GEN F = cgetg(3, t_MAT);
3918
  gel(F,1) = vecsmall_append(gel(z,1), off+u);
3919
  gel(F,2) = vecsmall_append(gel(z,2), e);
3920
  gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3921
}
3922

3923
/* add u^-1 v^-1 */
3924
static void
3925
addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
3926
{
3927
  pari_sp av = avma;
3928
  long off = r->prmax+1;
3929
  GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3930
  GEN F = cgetg(3, t_MAT);
3931
  gel(F,1) = vecsmall_concat(gel(z,1), P);
3932
  gel(F,2) = vecsmall_concat(gel(z,2), E);
3933
  gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3934
}
3935

3936
/*
3937
Let p=C^2+c
3938
Solve h = (C+x)*(C+a)-p = 0 [mod l]
3939
h= -c+x*(C+a)+C*a = 0  [mod l]
3940
x = (c-C*a)/(C+a) [mod l]
3941
h = -c+C*(x+a)+a*x
3942
*/
3943

3944
GEN
3945
Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3946
{
3947
  pari_sp ltop = avma;
3948
  long th, n = lg(pi)-1;
3949
  long i, j;
3950
  GEN sieve = zero_zv(a+2)+1;
3951
  GEN L = cgetg(1+a+2, t_VEC);
3952
  pari_sp av = avma;
3953
  long rel = 1;
3954
  GEN z, h;
3955
  h = addis(C,a);
3956
  if ((z = Z_issmooth_fact(h, prmax)))
3957
  {
3958
    gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
3959
    av = avma;
3960
  }
3961
  for (i=1; i<=n; i++)
3962
  {
3963
    ulong li = pi[i], s = sz[i], al = a % li;
3964
    ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3965
    if (!iv) continue;
3966
    u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3967
    for(j = u; j<=a; j+=li)
3968
      sieve[j] += s;
3969
  }
3970
  if (a)
3971
  {
3972
    long e = expi(mulis(C,a));
3973
    th = e - expu(e) - 1;
3974
  } else th = -1;
3975
  for (j=0; j<a; j++)
3976
    if (sieve[j]>=th)
3977
    {
3978
      GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3979
      if ((z = Z_issmooth_fact(h, prmax)))
3980
      {
3981
        gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
3982
        av = avma;
3983
      } else set_avma(av);
3984
    }
3985
  /* j = a */
3986
  if (sieve[a]>=th)
3987
  {
3988
    GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3989
    if ((z = Z_issmooth_fact(h, prmax)))
3990
      gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
3991
  }
3992
  setlg(L, rel);
3993
  return gerepilecopy(ltop, L);
3994
}
3995

3996
static long
3997
Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3998
{
3999
  struct pari_mt pt;
4000
  long i, j, nb = 0;
4001
  GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
4002
               mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
4003
  long running, pending = 0;
4004
  GEN W = zerovec(r->nbgen);
4005
  mt_queue_start_lim(&pt, worker, r->nbgen);
4006
  for (i = 0; (running = (i < r->nbgen)) || pending; i++)
4007
  {
4008
    GEN done;
4009
    long idx;
4010
    mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
4011
    done = mt_queue_get(&pt, &idx, &pending);
4012
    if (!done || lg(done)==1) continue;
4013
    gel(W, idx+1) = done;
4014
    nb += lg(done)-1;
4015
    if (DEBUGLEVEL && (i&127)==0)
4016
      err_printf("%ld%% ",100*nb/r->nbmax);
4017
  }
4018
  mt_queue_end(&pt);
4019
  for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
4020
  {
4021
    long ll, m;
4022
    GEN L = gel(W,j);
4023
    if (isintzero(L)) continue;
4024
    ll = lg(L);
4025
    for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
4026
    {
4027
      GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
4028
      if (v[1] == 1)
4029
        addifsmooth1(r, h, v[2], v[3]);
4030
      else
4031
        addifsmooth2(r, h, v[2], v[3]);
4032
    }
4033
  }
4034
  return j;
4035
}
4036

4037
static GEN
4038
ECP_psi(GEN x, GEN y)
4039
{
4040
  long prec = realprec(x);
4041
  GEN lx = glog(x, prec), ly = glog(y, prec);
4042
  GEN u = gdiv(lx, ly);
4043
  return gpow(u, gneg(u),prec);
4044
}
4045

4046
struct computeG
4047
{
4048
  GEN C;
4049
  long bnd, nbi;
4050
};
4051

4052
static GEN
4053
_computeG(void *E, GEN gen)
4054
{
4055
  struct computeG * d = (struct computeG *) E;
4056
  GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
4057
  return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
4058
}
4059

4060
static long
4061
compute_nbgen(GEN C, long bnd, long nbi)
4062
{
4063
  struct computeG d;
4064
  d.C = shifti(C, 1);
4065
  d.bnd = bnd;
4066
  d.nbi = nbi;
4067
  return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
4068
}
4069

4070
static GEN
4071
_psi(void*E, GEN y)
4072
{
4073
  GEN lx = (GEN) E;
4074
  long prec = realprec(lx);
4075
  GEN ly = glog(y, prec);
4076
  GEN u = gdiv(lx, ly);
4077
  return gsub(gdiv(y ,ly), gpow(u, u, prec));
4078
}
4079

4080
static GEN
4081
opt_param(GEN x, long prec)
4082
{
4083
  return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
4084
}
4085

4086
static GEN
4087
check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
4088
{
4089
  pari_sp av = avma;
4090
  long lM = lg(M)-1, nbcol = lM;
4091
  long tbs = maxss(1, expu(nbg/expi(m)));
4092
  for (;;)
4093
  {
4094
    GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
4095
    GEN tab;
4096
    long i, f=0;
4097
    long l = lg(K), lm = lgefint(m);
4098
    GEN idx = diviiexact(subiu(p,1),m), g;
4099
    pari_timer ti;
4100
    if (DEBUGLEVEL) timer_start(&ti);
4101
    for(i=1; i<l; i++)
4102
      if (signe(gel(K,i)))
4103
        break;
4104
    g = Fp_pow(utoi(i), idx, p);
4105
    tab = Fp_pow_init(g, p, tbs, p);
4106
    K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
4107
    for(i=1; i<l; i++)
4108
    {
4109
      GEN k = gel(K,i);
4110
      GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
4111
      if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
4112
        gel(K,i) = cgetineg(lm);
4113
      else
4114
        f++;
4115
    }
4116
    if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
4117
    if(f > (nbg>>1)) return gerepileupto(av, K);
4118
    for(i=1; i<=nbcol; i++)
4119
    {
4120
      long a = 1+random_Fl(lM);
4121
      swap(gel(M,a),gel(M,i));
4122
    }
4123
    if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
4124
  }
4125
}
4126

4127
static GEN
4128
Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
4129
{
4130
  pari_sp av=avma;
4131
  GEN aa = gen_1;
4132
  long AV = 0;
4133
  for(;;)
4134
  {
4135
    GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
4136
    GEN F = gel(A,1), E = gel(A,2);
4137
    GEN Ao = gen_0;
4138
    long i, l = lg(F);
4139
    for(i=1; i<l; i++)
4140
    {
4141
      GEN Ki = gel(K,F[i]);
4142
      if (signe(Ki)<0) break;
4143
      Ao = addii(Ao, mulis(Ki, E[i]));
4144
    }
4145
    if (i==l) return Fp_divu(Ao, AV, m);
4146
    aa = gerepileuptoint(av, aa);
4147
  }
4148
}
4149

4150
static GEN
4151
Fp_log_index(GEN a, GEN b, GEN m, GEN p)
4152
{
4153
  pari_sp av = avma, av2;
4154
  long i, j, nbi, nbr = 0, nbrow, nbg;
4155
  GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
4156
  pari_timer ti;
4157
  struct Fp_log_rel r;
4158
  ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
4159
  ulong bnd = 4*bnds;
4160
  if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
4161

4162
  p_1 = subiu(p,1);
4163
  if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
4164
    m = diviiexact(p_1, Z_ppo(p_1, m));
4165
  pr = primes_upto_zv(bnd);
4166
  nbi = lg(pr)-1;
4167
  C = sqrtremi(p, &c);
4168
  av2 = avma;
4169
  for (i = 1; i <= nbi; ++i)
4170
  {
4171
    ulong lp = pr[i];
4172
    while (lp <= bnd)
4173
    {
4174
      nbr++;
4175
      lp *= pr[i];
4176
    }
4177
  }
4178
  pi = cgetg(nbr+1,t_VECSMALL);
4179
  Ci = cgetg(nbr+1,t_VECSMALL);
4180
  ci = cgetg(nbr+1,t_VECSMALL);
4181
  sz = cgetg(nbr+1,t_VECSMALL);
4182
  for (i = 1, j = 1; i <= nbi; ++i)
4183
  {
4184
    ulong lp = pr[i], sp = expu(2*lp-1);
4185
    while (lp <= bnd)
4186
    {
4187
      pi[j] = lp;
4188
      Ci[j] = umodiu(C, lp);
4189
      ci[j] = umodiu(c, lp);
4190
      sz[j] = sp;
4191
      lp *= pr[i];
4192
      j++;
4193
    }
4194
  }
4195
  r.nbrel = 0;
4196
  r.nbgen = compute_nbgen(C, bnd, nbi);
4197
  r.nbmax = 2*(nbi+r.nbgen);
4198
  r.rel = cgetg(r.nbmax+1,t_VEC);
4199
  r.prmax = pr[nbi];
4200
  if (DEBUGLEVEL)
4201
  {
4202
    err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
4203
    timer_start(&ti);
4204
  }
4205
  nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
4206
  nbrow = r.prmax + nbg;
4207
  if (DEBUGLEVEL)
4208
  {
4209
    err_printf("\n");
4210
    timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
4211
  }
4212
  setlg(r.rel,r.nbrel+1);
4213
  r.rel = gerepilecopy(av2, r.rel);
4214
  K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
4215
  if (DEBUGLEVEL) timer_start(&ti);
4216
  Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
4217
  if (DEBUGLEVEL) timer_printf(&ti," log element");
4218
  Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
4219
  if (DEBUGLEVEL) timer_printf(&ti," log generator");
4220
  d = gcdii(Ao,Bo);
4221
  l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
4222
  if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
4223
  return gerepileuptoint(av, l);
4224
}
4225

4226
static int
4227
Fp_log_use_index(long e, long p)
4228
{
4229
  return (e >= 27 && 20*(p+6)<=e*e);
4230
}
4231

4232
/* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
4233
static GEN
4234
Fp_easylog(void *E, GEN a, GEN g, GEN ord)
4235
{
4236
  pari_sp av = avma;
4237
  GEN p = (GEN)E;
4238
  /* assume a reduced mod p, p not necessarily prime */
4239
  if (equali1(a)) return gen_0;
4240
  /* p > 2 */
4241
  if (equalii(subiu(p,1), a))  /* -1 */
4242
  {
4243
    pari_sp av2;
4244
    GEN t;
4245
    ord = get_arith_Z(ord);
4246
    if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
4247
    t = shifti(ord,-1); /* only possible solution */
4248
    av2 = avma;
4249
    if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
4250
    set_avma(av2); return gerepileuptoint(av, t);
4251
  }
4252
  if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
4253
    return Fp_log_index(a, g, ord, p);
4254
  return gc_NULL(av); /* not easy */
4255
}
4256

4257
GEN
4258
Fp_log(GEN a, GEN g, GEN ord, GEN p)
4259
{
4260
  GEN v = get_arith_ZZM(ord);
4261
  GEN F = gmael(v,2,1);
4262
  long lF = lg(F)-1, lmax;
4263
  if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
4264
  lmax = expi(gel(F,lF));
4265
  if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
4266
    v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
4267
  return gen_PH_log(a,g,v,(void*)p,&Fp_star);
4268
}
4269

4270
static ulong
4271
Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
4272
{
4273
  ulong i, h=1;
4274
  for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
4275
    if(a==h) return i;
4276
  return ~0UL;
4277
}
4278

4279
static ulong
4280
Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4281
{
4282
  ulong i, h=1;
4283
  for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
4284
    if(a==h) return i;
4285
  return ~0UL;
4286
}
4287

4288
static ulong
4289
Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
4290
{
4291
  pari_sp av = avma;
4292
  GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
4293
  return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
4294
}
4295

4296
ulong
4297
Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4298
{
4299
  if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
4300
  return Fl_log_Fp(a, g, ord, p);
4301
}
4302

4303
ulong
4304
Fl_log(ulong a, ulong g, ulong ord, ulong p)
4305
{
4306
  if (ord <= 200)
4307
  return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
4308
                      : Fl_log_naive(a, g, ord, p);
4309
  return Fl_log_Fp(a, g, ord, p);
4310
}
4311

4312
/* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
4313
 * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
4314
static GEN
4315
znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
4316
{
4317
  long l = lg(P) - 1, e = E[l];
4318
  GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
4319
  GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
4320

4321
  if (l == 1) {
4322
    hpe = h;
4323
    gpe = g;
4324
  } else {
4325
    hpe = modii(h, pe);
4326
    gpe = modii(g, pe);
4327
  }
4328
  if (e == 1) {
4329
    hp = hpe;
4330
    gp = gpe;
4331
  } else {
4332
    hp = remii(hpe, p);
4333
    gp = remii(gpe, p);
4334
  }
4335
  if (hp == gen_0 || gp == gen_0) return NULL;
4336
  if (absequaliu(p, 2))
4337
  {
4338
    GEN n = int2n(e);
4339
    ogpe = Zp_order(gpe, gen_2, e, n);
4340
    a = Fp_log(hpe, gpe, ogpe, n);
4341
    if (typ(a) != t_INT) return NULL;
4342
  }
4343
  else
4344
  { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
4345
       is trivial */
4346
    /* [order(gp), factor(order(gp))] */
4347
    GEN v = Fp_factored_order(gp, subiu(p,1), p);
4348
    GEN ogp = gel(v,1);
4349
    if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
4350
    a = Fp_log(hp, gp, v, p);
4351
    if (typ(a) != t_INT) return NULL;
4352
    if (e == 1) ogpe = ogp;
4353
    else
4354
    { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
4355
      /* use p-adic log: O(log p + e) mul*/
4356
      long vpogpe, vpohpe;
4357

4358
      hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
4359
      gpe = Fp_pow(gpe, ogp, pe);
4360
      /* g,h = 1 mod p; compute b s.t. h = g^b */
4361

4362
      /* v_p(order g mod pe) */
4363
      vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
4364
      /* v_p(order h mod pe) */
4365
      vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
4366
      if (vpohpe > vpogpe) return NULL;
4367

4368
      ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
4369
      if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
4370
      b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
4371
      a = addii(a, mulii(ogp, padic_to_Q(b)));
4372
    }
4373
  }
4374
  /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
4375
  if (l == 1) return a;
4376

4377
  N = diviiexact(N, pe); /* make N coprime to p */
4378
  h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
4379
  g = Fp_pow(g, modii(ogpe, phi), N);
4380
  setlg(P, l); /* remove last element */
4381
  setlg(E, l);
4382
  b = znlog_rec(h, g, N, P, E, PHI);
4383
  if (!b) return NULL;
4384
  return addmulii(a, b, ogpe);
4385
}
4386

4387
static GEN
4388
get_PHI(GEN P, GEN E)
4389
{
4390
  long i, l = lg(P);
4391
  GEN PHI = cgetg(l, t_VEC);
4392
  gel(PHI,1) = gen_1;
4393
  for (i=1; i<l-1; i++)
4394
  {
4395
    GEN t, p = gel(P,i);
4396
    long e = E[i];
4397
    t = mulii(powiu(p, e-1), subiu(p,1));
4398
    if (i > 1) t = mulii(t, gel(PHI,i));
4399
    gel(PHI,i+1) = t;
4400
  }
4401
  return PHI;
4402
}
4403

4404
GEN
4405
znlog(GEN h, GEN g, GEN o)
4406
{
4407
  pari_sp av = avma;
4408
  GEN N, fa, P, E, x;
4409
  switch (typ(g))
4410
  {
4411
    case t_PADIC:
4412
    {
4413
      GEN p = gel(g,2);
4414
      long v = valp(g);
4415
      if (v < 0) pari_err_DIM("znlog");
4416
      if (v > 0) {
4417
        long k = gvaluation(h, p);
4418
        if (k % v) return cgetg(1,t_VEC);
4419
        k /= v;
4420
        if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
4421
        set_avma(av); return stoi(k);
4422
      }
4423
      N = gel(g,3);
4424
      g = Rg_to_Fp(g, N);
4425
      break;
4426
    }
4427
    case t_INTMOD:
4428
      N = gel(g,1);
4429
      g = gel(g,2); break;
4430
    default: pari_err_TYPE("znlog", g);
4431
      return NULL; /* LCOV_EXCL_LINE */
4432
  }
4433
  if (equali1(N)) { set_avma(av); return gen_0; }
4434
  h = Rg_to_Fp(h, N);
4435
  if (o) return gerepileupto(av, Fp_log(h, g, o, N));
4436
  fa = Z_factor(N);
4437
  P = gel(fa,1);
4438
  E = vec_to_vecsmall(gel(fa,2));
4439
  x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
4440
  if (!x) { set_avma(av); return cgetg(1,t_VEC); }
4441
  return gerepileuptoint(av, x);
4442
}
4443

4444
GEN
4445
Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
4446
{
4447
  if (lgefint(p)==3)
4448
  {
4449
    long nn = itos_or_0(n);
4450
    if (nn)
4451
    {
4452
      ulong pp = p[2];
4453
      ulong uz;
4454
      ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
4455
      if (r==ULONG_MAX) return NULL;
4456
      if (zeta) *zeta = utoi(uz);
4457
      return utoi(r);
4458
    }
4459
  }
4460
  a = modii(a,p);
4461
  if (!signe(a))
4462
  {
4463
    if (zeta) *zeta = gen_1;
4464
    if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
4465
    return gen_0;
4466
  }
4467
  if (absequaliu(n,2))
4468
  {
4469
    if (zeta) *zeta = subiu(p,1);
4470
    return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
4471
  }
4472
  return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
4473
}
4474

4475
/*********************************************************************/
4476
/**                                                                 **/
4477
/**                    FUNDAMENTAL DISCRIMINANTS                    **/
4478
/**                                                                 **/
4479
/*********************************************************************/
4480
static long
4481
fa_isfundamental(GEN F)
4482
{
4483
  GEN P = gel(F,1), E = gel(F,2);
4484
  long i, s, l = lg(P);
4485

4486
  if (l == 1) return 1;
4487
  s = signe(gel(P,1)); /* = signe(x) */
4488
  if (!s) return 0;
4489
  if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
4490
  if (l == 1) return 0;
4491
  if (!absequaliu(gel(P,1), 2))
4492
    i = 1; /* need x = 1 mod 4 */
4493
  else
4494
  {
4495
    i = 2;
4496
    switch(itou(gel(E,1)))
4497
    {
4498
      case 2: s = -s; break; /* need x/4 = 3 mod 4 */
4499
      case 3: s = 0; break; /* no condition mod 4 */
4500
      default: return 0;
4501
    }
4502
  }
4503
  for(; i < l; i++)
4504
  {
4505
    if (!equali1(gel(E,i))) return 0;
4506
    if (s && Mod4(gel(P,i)) == 3) s = -s;
4507
  }
4508
  return s >= 0;
4509
}
4510
long
4511
isfundamental(GEN x)
4512
{
4513
  if (typ(x) != t_INT)
4514
  {
4515
    pari_sp av = avma;
4516
    long v = fa_isfundamental(check_arith_all(x,"isfundamental"));
4517
    return gc_long(av,v);
4518
  }
4519
  return Z_isfundamental(x);
4520
}
4521

4522
/* x fundamental ? */
4523
long
4524
uposisfundamental(ulong x)
4525
{
4526
  ulong r = x & 15; /* x mod 16 */
4527
  if (!r) return 0;
4528
  switch(r & 3)
4529
  { /* x mod 4 */
4530
    case 0: return (r == 4)? 0: uissquarefree(x >> 2);
4531
    case 1: return uissquarefree(x);
4532
    default: return 0;
4533
  }
4534
}
4535
/* -x fundamental ? */
4536
long
4537
unegisfundamental(ulong x)
4538
{
4539
  ulong r = x & 15; /* x mod 16 */
4540
  if (!r) return 0;
4541
  switch(r & 3)
4542
  { /* x mod 4 */
4543
    case 0: return (r == 12)? 0: uissquarefree(x >> 2);
4544
    case 3: return uissquarefree(x);
4545
    default: return 0;
4546
  }
4547
}
4548
long
4549
sisfundamental(long x)
4550
{ return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
4551

4552
long
4553
Z_isfundamental(GEN x)
4554
{
4555
  long r;
4556
  switch(lgefint(x))
4557
  {
4558
    case 2: return 0;
4559
    case 3: return signe(x) < 0? unegisfundamental(x[2])
4560
                               : uposisfundamental(x[2]);
4561
  }
4562
  r = mod16(x);
4563
  if (!r) return 0;
4564
  if ((r & 3) == 0)
4565
  {
4566
    pari_sp av;
4567
    r >>= 2; /* |x|/4 mod 4 */
4568
    if (signe(x) < 0) r = 4-r;
4569
    if (r == 1) return 0;
4570
    av = avma;
4571
    r = Z_issquarefree( shifti(x,-2) );
4572
    return gc_long(av, r);
4573
  }
4574
  r &= 3; /* |x| mod 4 */
4575
  if (signe(x) < 0) r = 4-r;
4576
  return (r==1) ? Z_issquarefree(x) : 0;
4577
}
4578

4579
static GEN
4580
fa_quaddisc(GEN f)
4581
{
4582
  GEN P = gel(f,1), E = gel(f,2), s = gen_1;
4583
  long i, l = lg(P);
4584
  for (i = 1; i < l; i++) /* possibly including -1 */
4585
    if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
4586
  if (Mod4(s) > 1) s = shifti(s,2);
4587
  return s;
4588
}
4589

4590
GEN
4591
quaddisc(GEN x)
4592
{
4593
  const pari_sp av = avma;
4594
  if (is_rational_t(typ(x))) x = factor(x);
4595
  else x = check_arith_all(x,"quaddisc");
4596
  return gerepileuptoint(av, fa_quaddisc(x));
4597
}
4598

4599
/*********************************************************************/
4600
/**                                                                 **/
4601
/**                              FACTORIAL                          **/
4602
/**                                                                 **/
4603
/*********************************************************************/
4604
GEN
4605
mulu_interval_step(ulong a, ulong b, ulong step)
4606
{
4607
  pari_sp av = avma;
4608
  ulong k, l, N, n;
4609
  long lx;
4610
  GEN x;
4611

4612
  if (!a) return gen_0;
4613
  if (step == 1) return mulu_interval(a, b);
4614
  n = 1 + (b-a) / step;
4615
  b -= (b-a) % step;
4616
  if (n < 61)
4617
  {
4618
    if (n == 1) return utoipos(a);
4619
    x = muluu(a,a+step); if (n == 2) return x;
4620
    for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
4621
    return gerepileuptoint(av, x);
4622
  }
4623
  /* step | b-a */
4624
  lx = 1; x = cgetg(2 + n/2, t_VEC);
4625
  N = b + a;
4626
  for (k = a;; k += step)
4627
  {
4628
    l = N - k; if (l <= k) break;
4629
    gel(x,lx++) = muluu(k,l);
4630
  }
4631
  if (l == k) gel(x,lx++) = utoipos(k);
4632
  setlg(x, lx);
4633
  return gerepileuptoint(av, ZV_prod(x));
4634
}
4635
/* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
4636
 * first is slower ... ] */
4637
GEN
4638
mulu_interval(ulong a, ulong b)
4639
{
4640
  pari_sp av = avma;
4641
  ulong k, l, N, n;
4642
  long lx;
4643
  GEN x;
4644

4645
  if (!a) return gen_0;
4646
  n = b - a + 1;
4647
  if (n < 61)
4648
  {
4649
    if (n == 1) return utoipos(a);
4650
    x = muluu(a,a+1); if (n == 2) return x;
4651
    for (k=a+2; k<=b; k++) x = mului(k,x);
4652
    return gerepileuptoint(av, x);
4653
  }
4654
  lx = 1; x = cgetg(2 + n/2, t_VEC);
4655
  N = b + a;
4656
  for (k = a;; k++)
4657
  {
4658
    l = N - k; if (l <= k) break;
4659
    gel(x,lx++) = muluu(k,l);
4660
  }
4661
  if (l == k) gel(x,lx++) = utoipos(k);
4662
  setlg(x, lx);
4663
  return gerepileuptoint(av, ZV_prod(x));
4664
}
4665
GEN
4666
muls_interval(long a, long b)
4667
{
4668
  pari_sp av = avma;
4669
  long lx, k, l, N, n = b - a + 1;
4670
  GEN x;
4671

4672
  if (a <= 0 && b >= 0) return gen_0;
4673
  if (n < 61)
4674
  {
4675
    x = stoi(a);
4676
    for (k=a+1; k<=b; k++) x = mulsi(k,x);
4677
    return gerepileuptoint(av, x);
4678
  }
4679
  lx = 1; x = cgetg(2 + n/2, t_VEC);
4680
  N = b + a;
4681
  for (k = a;; k++)
4682
  {
4683
    l = N - k; if (l <= k) break;
4684
    gel(x,lx++) = mulss(k,l);
4685
  }
4686
  if (l == k) gel(x,lx++) = stoi(k);
4687
  setlg(x, lx);
4688
  return gerepileuptoint(av, ZV_prod(x));
4689
}
4690

4691
GEN
4692
mpfact(long n)
4693
{
4694
  pari_sp av = avma;
4695
  GEN a, v;
4696
  long k;
4697
  if (n <= 12) switch(n)
4698
  {
4699
    case 0: case 1: return gen_1;
4700
    case 2: return gen_2;
4701
    case 3: return utoipos(6);
4702
    case 4: return utoipos(24);
4703
    case 5: return utoipos(120);
4704
    case 6: return utoipos(720);
4705
    case 7: return utoipos(5040);
4706
    case 8: return utoipos(40320);
4707
    case 9: return utoipos(362880);
4708
    case 10:return utoipos(3628800);
4709
    case 11:return utoipos(39916800);
4710
    case 12:return utoipos(479001600);
4711
    default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
4712
  }
4713
  v = cgetg(expu(n) + 2, t_VEC);
4714
  for (k = 1;; k++)
4715
  {
4716
    long m = n >> (k-1), l;
4717
    if (m <= 2) break;
4718
    l = (1 + (n >> k)) | 1;
4719
    /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
4720
    a = mulu_interval_step(l, m, 2);
4721
    gel(v,k) = k == 1? a: powiu(a, k);
4722
  }
4723
  a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
4724
  a = shifti(a, factorial_lval(n, 2));
4725
  return gerepileuptoint(av, a);
4726
}
4727

4728
ulong
4729
factorial_Fl(long n, ulong p)
4730
{
4731
  long k;
4732
  ulong v;
4733
  if (p <= (ulong)n) return 0;
4734
  v = Fl_powu(2, factorial_lval(n, 2), p);
4735
  for (k = 1;; k++)
4736
  {
4737
    long m = n >> (k-1), l, i;
4738
    ulong a = 1;
4739
    if (m <= 2) break;
4740
    l = (1 + (n >> k)) | 1;
4741
    /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
4742
    for (i=l; i<=m; i+=2)
4743
      a = Fl_mul(a, i, p);
4744
    v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
4745
  }
4746
  return v;
4747
}
4748

4749
GEN
4750
factorial_Fp(long n, GEN p)
4751
{
4752
  pari_sp av = avma;
4753
  long k;
4754
  GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
4755
  for (k = 1;; k++)
4756
  {
4757
    long m = n >> (k-1), l, i;
4758
    GEN a = gen_1;
4759
    if (m <= 2) break;
4760
    l = (1 + (n >> k)) | 1;
4761
    /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
4762
    for (i=l; i<=m; i+=2)
4763
      a = Fp_mulu(a, i, p);
4764
    v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
4765
    v = gerepileuptoint(av, v);
4766
  }
4767
  return v;
4768
}
4769

4770
/*******************************************************************/
4771
/**                                                               **/
4772
/**                      LUCAS & FIBONACCI                        **/
4773
/**                                                               **/
4774
/*******************************************************************/
4775
static void
4776
lucas(ulong n, GEN *a, GEN *b)
4777
{
4778
  GEN z, t, zt;
4779
  if (!n) { *a = gen_2; *b = gen_1; return; }
4780
  lucas(n >> 1, &z, &t); zt = mulii(z, t);
4781
  switch(n & 3) {
4782
    case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
4783
    case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
4784
    case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
4785
    case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
4786
  }
4787
}
4788

4789
GEN
4790
fibo(long n)
4791
{
4792
  pari_sp av = avma;
4793
  GEN a, b;
4794
  if (!n) return gen_0;
4795
  lucas((ulong)(labs(n)-1), &a, &b);
4796
  a = diviuexact(addii(shifti(a,1),b), 5);
4797
  if (n < 0 && !odd(n)) setsigne(a, -1);
4798
  return gerepileuptoint(av, a);
4799
}
4800

4801
/*******************************************************************/
4802
/*                                                                 */
4803
/*                      CONTINUED FRACTIONS                        */
4804
/*                                                                 */
4805
/*******************************************************************/
4806
static GEN
4807
icopy_lg(GEN x, long l)
4808
{
4809
  long lx = lgefint(x);
4810
  GEN y;
4811

4812
  if (lx >= l) return icopy(x);
4813
  y = cgeti(l); affii(x, y); return y;
4814
}
4815

4816
/* continued fraction of a/b. If y != NULL, stop when partial quotients
4817
 * differ from y */
4818
static GEN
4819
Qsfcont(GEN a, GEN b, GEN y, ulong k)
4820
{
4821
  GEN  z, c;
4822
  ulong i, l, ly = lgefint(b);
4823

4824
  /* times 1 / log2( (1+sqrt(5)) / 2 )  */
4825
  l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
4826
  if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
4827
  if (l > LGBITS) l = LGBITS;
4828

4829
  z = cgetg(l,t_VEC);
4830
  l--;
4831
  if (y) {
4832
    pari_sp av = avma;
4833
    if (l >= (ulong)lg(y)) l = lg(y)-1;
4834
    for (i = 1; i <= l; i++)
4835
    {
4836
      GEN q = gel(y,i);
4837
      gel(z,i) = q;
4838
      c = b; if (!gequal1(q)) c = mulii(q, b);
4839
      c = subii(a, c);
4840
      if (signe(c) < 0)
4841
      { /* partial quotient too large */
4842
        c = addii(c, b);
4843
        if (signe(c) >= 0) i++; /* by 1 */
4844
        break;
4845
      }
4846
      if (cmpii(c, b) >= 0)
4847
      { /* partial quotient too small */
4848
        c = subii(c, b);
4849
        if (cmpii(c, b) < 0) {
4850
          /* by 1. If next quotient is 1 in y, add 1 */
4851
          if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
4852
          i++;
4853
        }
4854
        break;
4855
      }
4856
      if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
4857
      a = b; b = c;
4858
    }
4859
  } else {
4860
    a = icopy_lg(a, ly);
4861
    b = icopy(b);
4862
    for (i = 1; i <= l; i++)
4863
    {
4864
      gel(z,i) = truedvmdii(a,b,&c);
4865
      if (c == gen_0) { i++; break; }
4866
      affii(c, a); cgiv(c); c = a;
4867
      a = b; b = c;
4868
    }
4869
  }
4870
  i--;
4871
  if (i > 1 && gequal1(gel(z,i)))
4872
  {
4873
    cgiv(gel(z,i)); --i;
4874
    gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
4875
  }
4876
  setlg(z,i+1); return z;
4877
}
4878

4879
static GEN
4880
sersfcont(GEN a, GEN b, long k)
4881
{
4882
  long i, l = typ(a) == t_POL? lg(a): 3;
4883
  GEN y, c;
4884
  if (lg(b) > l) l = lg(b);
4885
  if (k > 0 && l > k+1) l = k+1;
4886
  y = cgetg(l,t_VEC);
4887
  for (i=1; i<l; i++)
4888
  {
4889
    gel(y,i) = poldivrem(a,b,&c);
4890
    if (gequal0(c)) { i++; break; }
4891
    a = b; b = c;
4892
  }
4893
  setlg(y, i); return y;
4894
}
4895

4896
GEN
4897
gboundcf(GEN x, long k)
4898
{
4899
  pari_sp av;
4900
  long tx = typ(x), e;
4901
  GEN y, a, b, c;
4902

4903
  if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
4904
  if (is_scalar_t(tx))
4905
  {
4906
    if (gequal0(x)) return mkvec(gen_0);
4907
    switch(tx)
4908
    {
4909
      case t_INT: return mkveccopy(x);
4910
      case t_REAL:
4911
        av = avma;
4912
        c = mantissa_real(x,&e);
4913
        if (e < 0) pari_err_PREC("gboundcf");
4914
        y = int2n(e);
4915
        a = Qsfcont(c,y, NULL, k);
4916
        b = addsi(signe(x), c);
4917
        return gerepilecopy(av, Qsfcont(b,y, a, k));
4918

4919
      case t_FRAC:
4920
        av = avma;
4921
        return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
4922
    }
4923
    pari_err_TYPE("gboundcf",x);
4924
  }
4925

4926
  switch(tx)
4927
  {
4928
    case t_POL: return mkveccopy(x);
4929
    case t_SER:
4930
      av = avma;
4931
      return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
4932
    case t_RFRAC:
4933
      av = avma;
4934
      return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
4935
  }
4936
  pari_err_TYPE("gboundcf",x);
4937
  return NULL; /* LCOV_EXCL_LINE */
4938
}
4939

4940
static GEN
4941
sfcont2(GEN b, GEN x, long k)
4942
{
4943
  pari_sp av = avma;
4944
  long lb = lg(b), tx = typ(x), i;
4945
  GEN y,p1;
4946

4947
  if (k)
4948
  {
4949
    if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
4950
    lb = k+1;
4951
  }
4952
  y = cgetg(lb,t_VEC);
4953
  if (lb==1) return y;
4954
  if (is_scalar_t(tx))
4955
  {
4956
    if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
4957
  }
4958
  else if (tx == t_SER) x = ser2rfrac_i(x);
4959

4960
  if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
4961
  for (i = 1;;)
4962
  {
4963
    if (tx == t_REAL)
4964
    {
4965
      long e = expo(x);
4966
      if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
4967
      gel(y,i) = floorr(x);
4968
      p1 = subri(x, gel(y,i));
4969
    }
4970
    else
4971
    {
4972
      gel(y,i) = gfloor(x);
4973
      p1 = gsub(x, gel(y,i));
4974
    }
4975
    if (++i >= lb) break;
4976
    if (gequal0(p1)) break;
4977
    x = gdiv(gel(b,i),p1);
4978
  }
4979
  setlg(y,i);
4980
  return gerepilecopy(av,y);
4981
}
4982

4983
GEN
4984
gcf(GEN x) { return gboundcf(x,0); }
4985
GEN
4986
gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4987
GEN
4988
contfrac0(GEN x, GEN b, long nmax)
4989
{
4990
  long tb;
4991

4992
  if (!b) return gboundcf(x,nmax);
4993
  tb = typ(b);
4994
  if (tb == t_INT) return gboundcf(x,itos(b));
4995
  if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4996
  if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4997
  return sfcont2(b,x,nmax);
4998
}
4999

5000
GEN
5001
contfracpnqn(GEN x, long n)
5002
{
5003
  pari_sp av = avma;
5004
  long i, lx = lg(x);
5005
  GEN M,A,B, p0,p1, q0,q1;
5006

5007
  if (lx == 1)
5008
  {
5009
    if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
5010
    if (n >= 0) return cgetg(1,t_MAT);
5011
    return matid(2);
5012
  }
5013
  switch(typ(x))
5014
  {
5015
    case t_VEC: case t_COL: A = x; B = NULL; break;
5016
    case t_MAT:
5017
      switch(lgcols(x))
5018
      {
5019
        case 2: A = row(x,1); B = NULL; break;
5020
        case 3: A = row(x,2); B = row(x,1); break;
5021
        default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
5022
                 return NULL; /*LCOV_EXCL_LINE*/
5023
      }
5024
      break;
5025
    default: pari_err_TYPE("pnqn",x);
5026
      return NULL; /*LCOV_EXCL_LINE*/
5027
  }
5028
  p1 = gel(A,1);
5029
  q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
5030
  if (n >= 0)
5031
  {
5032
    lx = minss(lx, n+2);
5033
    if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
5034
  }
5035
  else if (lx == 2)
5036
    return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
5037
  /* lx >= 3 */
5038
  p0 = gen_1;
5039
  q0 = gen_0; /* p[-1], q[-1] */
5040
  M = cgetg(lx, t_MAT);
5041
  gel(M,1) = mkcol2(p1,q1);
5042
  for (i=2; i<lx; i++)
5043
  {
5044
    GEN a = gel(A,i), p2,q2;
5045
    if (B) {
5046
      GEN b = gel(B,i);
5047
      p0 = gmul(b,p0);
5048
      q0 = gmul(b,q0);
5049
    }
5050
    p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
5051
    q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
5052
    gel(M,i) = mkcol2(p1,q1);
5053
  }
5054
  if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
5055
  return gerepilecopy(av, M);
5056
}
5057
GEN
5058
pnqn(GEN x) { return contfracpnqn(x,-1); }
5059
/* x = [a0, ..., an] from gboundcf, n >= 0;
5060
 * return [[p0, ..., pn], [q0,...,qn]] */
5061
GEN
5062
ZV_allpnqn(GEN x)
5063
{
5064
  long i, lx = lg(x);
5065
  GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
5066

5067
  gel(v,1) = P = cgetg(lx, t_VEC);
5068
  gel(v,2) = Q = cgetg(lx, t_VEC);
5069
  p0 = gen_1; q0 = gen_0;
5070
  gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
5071
  for (i=2; i<lx; i++)
5072
  {
5073
    GEN a = gel(x,i);
5074
    gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
5075
    gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
5076
  }
5077
  return v;
5078
}
5079

5080
/* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
5081
static GEN
5082
mod_to_frac(GEN x, GEN N, GEN B)
5083
{
5084
  GEN a, b, A;
5085
  if (B) A = divii(shifti(N, -1), B);
5086
  else
5087
  {
5088
    A = sqrti(shifti(N, -1));
5089
    B = A;
5090
  }
5091
  if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
5092
  return equali1(b)? a: mkfrac(a,b);
5093
}
5094

5095
static GEN
5096
mod_to_rfrac(GEN x, GEN N, long B)
5097
{
5098
  GEN a, b;
5099
  long A, d = degpol(N);
5100
  if (B >= 0) A = d-1 - B;
5101
  else
5102
  {
5103
    B = d >> 1;
5104
    A = odd(d)? B : B-1;
5105
  }
5106
  if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
5107
  if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
5108
  return gdiv(a,b);
5109
}
5110

5111
/* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
5112
 * of the continued fraction of x with b <= k maximal */
5113
static GEN
5114
bestappr_frac(GEN x, GEN k)
5115
{
5116
  pari_sp av;
5117
  GEN p0, p1, p, q0, q1, q, a, y;
5118

5119
  if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
5120
  av = avma; y = x;
5121
  p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
5122
  q1 = gen_0; q0 = gen_1;
5123
  x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
5124
  for(;;)
5125
  {
5126
    x = ginv(x); /* > 1 */
5127
    a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
5128
    if (cmpii(a,k) > 0)
5129
    { /* next partial quotient will overflow limits */
5130
      GEN n, d;
5131
      a = divii(subii(k, q1), q0);
5132
      p = addii(mulii(a,p0), p1); p1=p0; p0=p;
5133
      q = addii(mulii(a,q0), q1); q1=q0; q0=q;
5134
      /* compare |y-p0/q0|, |y-p1/q1| */
5135
      n = gel(y,1);
5136
      d = gel(y,2);
5137
      if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
5138
                   mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
5139
                   { p1 = p0; q1 = q0; }
5140
      break;
5141
    }
5142
    p = addii(mulii(a,p0), p1); p1=p0; p0=p;
5143
    q = addii(mulii(a,q0), q1); q1=q0; q0=q;
5144

5145
    if (cmpii(q0,k) > 0) break;
5146
    x = gsub(x,a); /* 0 <= x < 1 */
5147
    if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
5148

5149
  }
5150
  return gerepileupto(av, gdiv(p1,q1));
5151
}
5152
/* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
5153
 * of the continued fraction of x with b <= k maximal */
5154
static GEN
5155
bestappr_real(GEN x, GEN k)
5156
{
5157
  pari_sp av = avma;
5158
  GEN kr, p0, p1, p, q0, q1, q, a, y = x;
5159

5160
  p1 = gen_1; a = p0 = floorr(x);
5161
  q1 = gen_0; q0 = gen_1;
5162
  x = subri(x,a); /* 0 <= x < 1 */
5163
  if (!signe(x)) { cgiv(x); return a; }
5164
  kr = itor(k, realprec(x));
5165
  for(;;)
5166
  {
5167
    long d;
5168
    x = invr(x); /* > 1 */
5169
    if (cmprr(x,kr) > 0)
5170
    { /* next partial quotient will overflow limits */
5171
      a = divii(subii(k, q1), q0);
5172
      p = addii(mulii(a,p0), p1); p1=p0; p0=p;
5173
      q = addii(mulii(a,q0), q1); q1=q0; q0=q;
5174
      /* compare |y-p0/q0|, |y-p1/q1| */
5175
      if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
5176
                   mulir(q0, subri(mulir(q1,y), p1))) < 0)
5177
                   { p1 = p0; q1 = q0; }
5178
      break;
5179
    }
5180
    d = nbits2prec(expo(x) + 1);
5181
    if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
5182

5183
    a = truncr(x); /* truncr(x) will NOT raise e_PREC */
5184
    p = addii(mulii(a,p0), p1); p1=p0; p0=p;
5185
    q = addii(mulii(a,q0), q1); q1=q0; q0=q;
5186

5187
    if (cmpii(q0,k) > 0) break;
5188
    x = subri(x,a); /* 0 <= x < 1 */
5189
    if (!signe(x)) { p1 = p0; q1 = q0; break; }
5190
  }
5191
  if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
5192
  return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
5193
}
5194

5195
/* k t_INT or NULL */
5196
static GEN
5197
bestappr_Q(GEN x, GEN k)
5198
{
5199
  long lx, tx = typ(x), i;
5200
  GEN a, y;
5201

5202
  switch(tx)
5203
  {
5204
    case t_INT: return icopy(x);
5205
    case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
5206
    case t_REAL:
5207
      if (!signe(x)) return gen_0;
5208
      /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
5209
      i = bit_prec(x); if (i <= expo(x)) return NULL;
5210
      return bestappr_real(x, k? k: int2n(i));
5211

5212
    case t_INTMOD: {
5213
      pari_sp av = avma;
5214
      a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
5215
      return gerepilecopy(av, a);
5216
    }
5217
    case t_PADIC: {
5218
      pari_sp av = avma;
5219
      long v = valp(x);
5220
      a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
5221
      if (v) a = gmul(a, powis(gel(x,2), v));
5222
      return gerepilecopy(av, a);
5223
    }
5224

5225
    case t_COMPLEX: {
5226
      pari_sp av = avma;
5227
      y = cgetg(3, t_COMPLEX);
5228
      gel(y,2) = bestappr(gel(x,2), k);
5229
      gel(y,1) = bestappr(gel(x,1), k);
5230
      if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
5231
      return y;
5232
    }
5233
    case t_SER:
5234
      if (ser_isexactzero(x)) return gcopy(x);
5235
      /* fall through */
5236
    case t_POLMOD: case t_POL: case t_RFRAC:
5237
    case t_VEC: case t_COL: case t_MAT:
5238
      y = cgetg_copy(x, &lx);
5239
      if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
5240
      for (; i<lx; i++)
5241
      {
5242
        a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
5243
        gel(y,i) = a;
5244
      }
5245
      if (tx == t_POL) return normalizepol(y);
5246
      if (tx == t_SER) return normalize(y);
5247
      return y;
5248
  }
5249
  pari_err_TYPE("bestappr_Q",x);
5250
  return NULL; /* LCOV_EXCL_LINE */
5251
}
5252

5253
static GEN
5254
bestappr_ser(GEN x, long B)
5255
{
5256
  long dN, v = valp(x), lx = lg(x);
5257
  GEN t;
5258
  x = normalizepol(ser2pol_i(x, lx));
5259
  dN = lx-2;
5260
  if (v > 0)
5261
  {
5262
    x = RgX_shift_shallow(x, v);
5263
    dN += v;
5264
  }
5265
  else if (v < 0)
5266
  {
5267
    if (B >= 0) B = maxss(B+v, 0);
5268
  }
5269
  t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
5270
  if (!t) return NULL;
5271
  if (v < 0)
5272
  {
5273
    GEN a, b;
5274
    long vx;
5275
    if (typ(t) == t_POL) return RgX_mulXn(t, v);
5276
    /* t_RFRAC */
5277
    vx = varn(x);
5278
    a = gel(t,1);
5279
    b = gel(t,2);
5280
    v -= RgX_valrem(b, &b);
5281
    if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
5282
    if (v < 0) b = RgX_shift(b, -v);
5283
    else if (v > 0) {
5284
      if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
5285
      a = RgX_shift(a, v);
5286
    }
5287
    t = mkrfraccopy(a, b);
5288
  }
5289
  return t;
5290
}
5291
static GEN bestappr_RgX(GEN x, long B);
5292
/* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
5293
 * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
5294
static GEN
5295
bestappr_RgX(GEN x, long B)
5296
{
5297
  long i, lx, tx = typ(x);
5298
  GEN y, t;
5299
  switch(tx)
5300
  {
5301
    case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
5302
    case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
5303
      return gcopy(x);
5304

5305
    case t_RFRAC: {
5306
      pari_sp av = avma;
5307
      if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
5308
      x = rfrac_to_ser(x, 2*B+1);
5309
      t = bestappr_ser(x, B); if (!t) return NULL;
5310
      return gerepileupto(av, t);
5311
    }
5312
    case t_POLMOD: {
5313
      pari_sp av = avma;
5314
      t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
5315
      return gerepileupto(av, t);
5316
    }
5317
    case t_SER: {
5318
      pari_sp av = avma;
5319
      t = bestappr_ser(x, B); if (!t) return NULL;
5320
      return gerepileupto(av, t);
5321
    }
5322

5323
    case t_VEC: case t_COL: case t_MAT:
5324
      y = cgetg_copy(x, &lx);
5325
      if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
5326
      for (; i<lx; i++)
5327
      {
5328
        t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
5329
        gel(y,i) = t;
5330
      }
5331
      return y;
5332
  }
5333
  pari_err_TYPE("bestappr_RgX",x);
5334
  return NULL; /* LCOV_EXCL_LINE */
5335
}
5336

5337
/* allow k = NULL: maximal accuracy */
5338
GEN
5339
bestappr(GEN x, GEN k)
5340
{
5341
  pari_sp av = avma;
5342
  if (k) { /* replace by floor(k) */
5343
    switch(typ(k))
5344
    {
5345
      case t_INT:
5346
        break;
5347
      case t_REAL: case t_FRAC:
5348
        k = floor_safe(k); /* left on stack for efficiency */
5349
        if (!signe(k)) k = gen_1;
5350
        break;
5351
      default:
5352
        pari_err_TYPE("bestappr [bound type]", k);
5353
        break;
5354
    }
5355
  }
5356
  x = bestappr_Q(x, k);
5357
  if (!x) { set_avma(av); return cgetg(1,t_VEC); }
5358
  return x;
5359
}
5360
GEN
5361
bestapprPade(GEN x, long B)
5362
{
5363
  pari_sp av = avma;
5364
  GEN t = bestappr_RgX(x, B);
5365
  if (!t) { set_avma(av); return cgetg(1,t_VEC); }
5366
  return t;
5367
}
5368

5369
/***********************************************************************/
5370
/**                                                                   **/
5371
/**         FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS)         **/
5372
/**                                                                   **/
5373
/***********************************************************************/
5374

5375
static GEN
5376
get_quad(GEN f, GEN pol, long r)
5377
{
5378
  GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
5379
  return mkquad(pol, r? subii(p1,q1): p1, q1);
5380
}
5381

5382
/* replace f by f * [a,1; 1,0] */
5383
static void
5384
update_f(GEN f, GEN a)
5385
{
5386
  GEN p1;
5387
  p1 = gcoeff(f,1,1);
5388
  gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
5389
  gcoeff(f,1,2) = p1;
5390

5391
  p1 = gcoeff(f,2,1);
5392
  gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
5393
  gcoeff(f,2,2) = p1;
5394
}
5395

5396
/*
5397
 * fm is a vector of matrices and i an index
5398
 * the bits of i give the non-zero entries
5399
 * the product of the non zero entries is the
5400
 * actual result.
5401
 * if i odd, fm[1] is implicitely [fm[1],1;1,0]
5402
 */
5403

5404
static void
5405
update_fm(GEN f, GEN a, long i)
5406
{
5407
  if (!odd(i))
5408
    gel(f,1) = a;
5409
  else
5410
  {
5411
    long v = vals(i+1), k;
5412
    GEN b = gel(f,1);
5413
    GEN u = mkmat22(addiu(mulii(a, b), 1), b, a, gen_1);
5414
    gel(f,1) = gen_0;
5415
    for (k = 1; k < v; k++)
5416
    {
5417
      u = ZM2_mul(gel(f, k+1), u);
5418
      gel(f,k+1) = gen_0; /* for gerepileall */
5419
    }
5420
    gel(f,v+1) = u;
5421
  }
5422
}
5423

5424
static GEN
5425
prod_fm(GEN f, long i)
5426
{
5427
  long v = vals(i);
5428
  GEN u;
5429
  long k;
5430
  if (!v) u = mkmat22(gel(f,1),gen_1,gen_1,gen_0);
5431
  else u = gel(f,v+1);
5432
  v++;
5433
  for (i>>=v, k = v+1; i; i>>=1, k++)
5434
    if (odd(i))
5435
      u = ZM2_mul(gel(f,k), u);
5436
  return u;
5437
}
5438

5439
GEN
5440
quadunit(GEN x)
5441
{
5442
  pari_sp av = avma, av2;
5443
  GEN pol, y, a, u, v, sqd, f;
5444
  long r, i = 1;
5445

5446
  check_quaddisc_real(x, &r, "quadunit");
5447
  pol = quadpoly(x);
5448
  sqd = sqrti(x); av2 = avma;
5449
  a = shifti(addui(r,sqd),-1);
5450
  f = zerovec(2+(expi(x)>>1));
5451
  gel(f,1) = a;
5452
  u = stoi(r); v = gen_2;
5453
  for(;;)
5454
  {
5455
    GEN u1, v1;
5456
    u1 = subii(mulii(a,v),u);
5457
    v1 = divii(subii(x,sqri(u1)),v);
5458
    if ( equalii(v,v1) ) {
5459
      f = prod_fm(f,i);
5460
      y = get_quad(f,pol,r);
5461
      update_f(f,a);
5462
      y = gdiv(get_quad(f,pol,r), conj_i(y));
5463
      break;
5464
    }
5465
    a = divii(addii(sqd,u1), v1);
5466
    if ( equalii(u,u1) ) {
5467
      y = get_quad(prod_fm(f,i),pol,r);
5468
      y = gdiv(y, conj_i(y));
5469
      break;
5470
    }
5471
    update_fm(f,a,i++);
5472
    u = u1; v = v1;
5473
    if (gc_needed(av2,2))
5474
    {
5475
      if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
5476
      gerepileall(av2,4, &a,&f,&u,&v);
5477
    }
5478
  }
5479
  if (signe(gel(y,3)) < 0) y = gneg(y);
5480
  return gerepileupto(av, y);
5481
}
5482

5483
GEN
5484
quadunit0(GEN x, long v)
5485
{
5486
  GEN y = quadunit(x);
5487
  if (v==-1) v = fetch_user_var("w");
5488
  setvarn(gel(y,1), v);
5489
  return y;
5490
}
5491

5492
GEN
5493
quadregulator(GEN x, long prec)
5494
{
5495
  pari_sp av = avma, av2;
5496
  GEN R, rsqd, u, v, sqd;
5497
  long r, Rexpo;
5498

5499
  check_quaddisc_real(x, &r, "quadregulator");
5500
  sqd = sqrti(x);
5501
  rsqd = gsqrt(x,prec);
5502
  Rexpo = 0; R = real2n(1, prec); /* = 2 */
5503
  av2 = avma;
5504
  u = stoi(r); v = gen_2;
5505
  for(;;)
5506
  {
5507
    GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
5508
    GEN v1 = divii(subii(x,sqri(u1)),v);
5509
    if (equalii(v,v1))
5510
    {
5511
      R = sqrr(R); shiftr_inplace(R, -1);
5512
      R = mulrr(R, divri(addir(u1,rsqd),v));
5513
      break;
5514
    }
5515
    if (equalii(u,u1))
5516
    {
5517
      R = sqrr(R); shiftr_inplace(R, -1);
5518
      break;
5519
    }
5520
    R = mulrr(R, divri(addir(u1,rsqd),v));
5521
    Rexpo += expo(R); setexpo(R,0);
5522
    u = u1; v = v1;
5523
    if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
5524
    if (gc_needed(av2,2))
5525
    {
5526
      if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
5527
      gerepileall(av2,3, &R,&u,&v);
5528
    }
5529
  }
5530
  R = logr_abs(divri(R,v));
5531
  if (Rexpo)
5532
  {
5533
    GEN t = mulsr(Rexpo, mplog2(prec));
5534
    shiftr_inplace(t, 1);
5535
    R = addrr(R,t);
5536
  }
5537
  return gerepileuptoleaf(av, R);
5538
}
5539

5540
/*************************************************************************/
5541
/**                                                                     **/
5542
/**                            CLASS NUMBER                             **/
5543
/**                                                                     **/
5544
/*************************************************************************/
5545

5546
int
5547
qfb_equal1(GEN f) { return equali1(gel(f,1)); }
5548

5549
static GEN qfi_pow(void *E, GEN f, GEN n)
5550
{ return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
5551
static GEN qfi_comp(void *E, GEN f, GEN g)
5552
{ return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
5553
static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
5554
                                         gidentical,qfb_equal1,NULL};
5555

5556
GEN
5557
qfi_order(GEN q, GEN o)
5558
{ return gen_order(q, o, NULL, &qfi_group); }
5559

5560
GEN
5561
qfi_log(GEN a, GEN g, GEN o)
5562
{ return gen_PH_log(a, g, o, NULL, &qfi_group); }
5563

5564
GEN
5565
qfi_Shanks(GEN a, GEN g, long n)
5566
{
5567
  pari_sp av = avma;
5568
  GEN T, X;
5569
  long rt_n, c;
5570

5571
  a = qfbred_i(a);
5572
  g = qfbred_i(g);
5573

5574
  rt_n = sqrt((double)n);
5575
  c = n / rt_n;
5576
  c = (c * rt_n < n + 1) ? c + 1 : c;
5577

5578
  T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
5579
  X = gen_Shanks(T, a, c, NULL, &qfi_group);
5580
  return X? gerepileuptoint(av, X): gc_NULL(av);
5581
}
5582

5583
GEN
5584
qfbclassno0(GEN x,long flag)
5585
{
5586
  switch(flag)
5587
  {
5588
    case 0: return map_proto_G(classno,x);
5589
    case 1: return map_proto_G(classno2,x);
5590
    default: pari_err_FLAG("qfbclassno");
5591
  }
5592
  return NULL; /* LCOV_EXCL_LINE */
5593
}
5594

5595
/* f^h = 1, return order(f). Set *pfao to its factorization */
5596
static GEN
5597
find_order(void *E, GEN f, GEN h, GEN *pfao)
5598
{
5599
  GEN v = gen_factored_order(f, h, E, &qfi_group);
5600
  *pfao = gel(v,2); return gel(v,1);
5601
}
5602

5603
static int
5604
ok_q(GEN q, GEN h, GEN d2, long r2)
5605
{
5606
  if (d2)
5607
  {
5608
    if (r2 <= 2 && !mpodd(q)) return 0;
5609
    return is_pm1(Z_ppo(q,d2));
5610
  }
5611
  else
5612
  {
5613
    if (r2 <= 1 && !mpodd(q)) return 0;
5614
    return is_pm1(Z_ppo(q,h));
5615
  }
5616
}
5617

5618
/* a,b given by their factorizations. Return factorization of lcm(a,b).
5619
 * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
5620
static GEN
5621
split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
5622
{
5623
  GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
5624
  GEN A = gen_1, B = gen_1;
5625
  long i, l = lg(P);
5626
  GEN E = cgetg(l, t_COL);
5627
  for (i=1; i<l; i++)
5628
  {
5629
    GEN p = gel(P,i);
5630
    long va = Z_pval(a,p);
5631
    long vb = Z_pval(b,p);
5632
    if (va < vb)
5633
    {
5634
      B = mulii(B,powiu(p,vb));
5635
      gel(E,i) = utoi(vb);
5636
    }
5637
    else
5638
    {
5639
      A = mulii(A,powiu(p,va));
5640
      gel(E,i) = utoi(va);
5641
    }
5642
  }
5643
  *pA = A;
5644
  *pB = B; return mkmat2(P,E);
5645
}
5646

5647
/* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
5648
static void
5649
update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
5650
{
5651
  GEN A,B, g1 = *pg1, d1 = *pd1;
5652
  *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
5653
  *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)),  qfbpow_i(f, diviiexact(o,B)));
5654
  *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
5655
}
5656

5657
/* Write x = Df^2, where D = fundamental discriminant,
5658
 * P^E = factorisation of conductor f, with E[i] >= 0 */
5659
static void
5660
corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
5661
{
5662
  long s = signe(x), l, i;
5663
  GEN fa = absZ_factor(x);
5664
  GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
5665

5666
  l = lg(P); d = gen_1;
5667
  for (i=1; i<l; i++)
5668
  {
5669
    if (E[i] & 1) d = mulii(d, gel(P,i));
5670
    E[i] >>= 1;
5671
  }
5672
  if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
5673
  *ptD = (s < 0)? negi(d): d;
5674
  *ptP = P;
5675
  *ptE = E;
5676
}
5677

5678
static GEN
5679
conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
5680
{
5681
  long l, i, s = signe(x);
5682
  GEN E, H, D, P, reg;
5683

5684
  corediscfact(x, xmod4, &D, &P, &E);
5685
  H = gen_1; l = lg(P);
5686
  /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
5687
  for (i=1; i<l; i++)
5688
  {
5689
    long e = E[i];
5690
    if (e)
5691
    {
5692
      GEN p = gel(P,i);
5693
      H = mulii(H, subis(p, kronecker(D,p)));
5694
      if (e >= 2) H = mulii(H, powiu(p,e-1));
5695
    }
5696
  }
5697

5698
  /* divide by [ O_K^* : O^* ] */
5699
  if (s < 0)
5700
  {
5701
    reg = NULL;
5702
    switch(itou_or_0(D))
5703
    {
5704
      case 4: H = shifti(H,-1); break;
5705
      case 3: H = divis(H,3); break;
5706
    }
5707
  } else {
5708
    reg = quadregulator(D,DEFAULTPREC);
5709
    if (!equalii(x,D))
5710
      H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
5711
  }
5712
  if (ptreg) *ptreg = reg;
5713
  *ptD = D; return H;
5714
}
5715

5716
static long
5717
two_rank(GEN x)
5718
{
5719
  GEN p = gel(absZ_factor(x),1);
5720
  long l = lg(p)-1;
5721
#if 0 /* positive disc not needed */
5722
  if (signe(x) > 0)
5723
  {
5724
    long i;
5725
    for (i=1; i<=l; i++)
5726
      if (mod4(gel(p,i)) == 3) { l--; break; }
5727
  }
5728
#endif
5729
  return l-1;
5730
}
5731

5732
static GEN
5733
sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
5734
/* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
5735
static GEN
5736
get_forms(GEN D, GEN *pL)
5737
{
5738
  const long MAXFORM = 20;
5739
  GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
5740
  GEN forms = vectrunc_init(MAXFORM+1);
5741
  long s, nforms = 0;
5742
  ulong p;
5743
  forprime_t S;
5744
  L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
5745
  s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
5746
  if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
5747
  if      (s < 10)   s = 200;
5748
  else if (s < 20)   s = 1000;
5749
  else if (s < 5000) s = 5000;
5750
  u_forprime_init(&S, 2, s);
5751
  while ( (p = u_forprime_next(&S)) )
5752
  {
5753
    long d, k = kroiu(D,p);
5754
    pari_sp av2;
5755
    if (!k) continue;
5756
    if (k > 0)
5757
    {
5758
      if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
5759
      d = p-1;
5760
    }
5761
    else
5762
      d = p+1;
5763
    av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
5764
  }
5765
  *pL = L; return forms;
5766
}
5767

5768
/* h ~ #G, return o = order of f, set fao = its factorization */
5769
static  GEN
5770
Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
5771
{
5772
  long s = minss(itos(sqrti(h)), 10000);
5773
  GEN T = gen_Shanks_init(f, s, E, &qfi_group);
5774
  GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
5775
  return find_order(E, f, addiu(v,1), pfao);
5776
}
5777

5778
/* if g = 1 in  G/<f> ? */
5779
static int
5780
equal1(void *E, GEN T, ulong N, GEN g)
5781
{ return !!gen_Shanks(T, g, N, E, &qfi_group); }
5782

5783
/* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
5784
 * FIXME: should be gen_order, but equal1 has the wrong prototype */
5785
static GEN
5786
relative_order(void *E, GEN a, GEN o, ulong N,  GEN T)
5787
{
5788
  pari_sp av = avma;
5789
  long i, l;
5790
  GEN m;
5791

5792
  m = get_arith_ZZM(o);
5793
  if (!m) pari_err_TYPE("gen_order [missing order]",a);
5794
  o = gel(m,1);
5795
  m = gel(m,2); l = lgcols(m);
5796
  for (i = l-1; i; i--)
5797
  {
5798
    GEN t, y, p = gcoeff(m,i,1);
5799
    long j, e = itos(gcoeff(m,i,2));
5800
    if (l == 2) {
5801
      t = gen_1;
5802
      y = a;
5803
    } else {
5804
      t = diviiexact(o, powiu(p,e));
5805
      y = powgi(a, t);
5806
    }
5807
    if (equal1(E, T,N,y)) o = t;
5808
    else {
5809
      for (j = 1; j < e; j++)
5810
      {
5811
        y = powgi(y, p);
5812
        if (equal1(E, T,N,y)) break;
5813
      }
5814
      if (j < e) {
5815
        if (j > 1) p = powiu(p, j);
5816
        o = mulii(t, p);
5817
      }
5818
    }
5819
  }
5820
  return gerepilecopy(av, o);
5821
}
5822

5823
/* h(x) for x<0 using Baby Step/Giant Step.
5824
 * Assumes G is not too far from being cyclic.
5825
 *
5826
 * Compute G^2 instead of G so as to kill most of the noncyclicity */
5827
GEN
5828
classno(GEN x)
5829
{
5830
  pari_sp av = avma;
5831
  long r2, k, s, i, l;
5832
  GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
5833
  void *E;
5834

5835
  if (signe(x) >= 0) return classno2(x);
5836

5837
  check_quaddisc(x, &s, &k, "classno");
5838
  if (abscmpiu(x,12) <= 0) return gen_1;
5839

5840
  Hf = conductor_part(x, k, &D, NULL);
5841
  if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
5842
  forms =  get_forms(D, &L);
5843
  r2 = two_rank(D);
5844
  hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
5845

5846
  l = lg(forms);
5847
  order_bound = const_vec(l-1, NULL);
5848
  E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
5849
  g1 = gel(forms,1);
5850
  gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
5851
  q = diviiround(hin, d1); /* approximate order of G/<g1> */
5852
  d2 = NULL; /* not computed yet */
5853
  if (is_pm1(q)) goto END;
5854
  for (i=2; i < l; i++)
5855
  {
5856
    GEN o, fao, a, F, fd, f = gel(forms,i);
5857
    fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
5858
    F = qfbpow_i(fd, q);
5859
    a = gel(F,1);
5860
    o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
5861
    /* f^(d1 q) = 1 */
5862
    fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
5863
    o = find_order(E, f, fao, &fao);
5864
    gel(order_bound,i) = o;
5865
    /* o = order of f, fao = factor(o) */
5866
    update_g1(&g1,&d1,&fad1, f,o,fao);
5867
    q = diviiround(hin, d1);
5868
    if (is_pm1(q)) goto END;
5869
  }
5870
  /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
5871
  if (expi(q) > 3)
5872
  { /* q large: compute d2, 2nd elt divisor */
5873
    ulong N, n = 2*itou(sqrti(d1));
5874
    GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
5875
    d2 = gen_1;
5876
    N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
5877
    for (i = 1; i < l; i++)
5878
    {
5879
      GEN d, f = gel(forms,i), B = gel(order_bound,i);
5880
      if (!B) B = find_order(E, f, fad1, /*junk*/&d);
5881
      f = qfbpow_i(f,d2);
5882
      if (equal1(E, T,N,f)) continue;
5883
      B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
5884
      /* f^B = 1 */
5885
      d = relative_order(E, f, B, N,T);
5886
      d2= mulii(d,d2);
5887
      D = mulii(d1,d2);
5888
      q = diviiround(hin,D);
5889
      if (is_pm1(q)) { d1 = D; goto END; }
5890
    }
5891
    /* very probably, d2 is the 2nd elementary divisor */
5892
    d1 = D; /* product of first two elt divisors */
5893
  }
5894
  /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
5895
   * 2-rank */
5896
  if (!ok_q(q,d1,d2,r2))
5897
  {
5898
    GEN q0 = q;
5899
    long d;
5900
    if (cmpii(mulii(q,d1), hin) < 0)
5901
    { /* try q = q0+1,-1,+2,-2 */
5902
      d = 1;
5903
      do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
5904
    }
5905
    else
5906
    { /* q0-1,+1,-2,+2  */
5907
      d = -1;
5908
      do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
5909
    }
5910
  }
5911
  d1 = mulii(d1,q);
5912

5913
END:
5914
  return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
5915
}
5916

5917
GEN
5918
quadclassno(GEN x)
5919
{
5920
  pari_sp av = avma;
5921
  GEN Hf, D;
5922
  long s, r;
5923
  check_quaddisc(x, &s, &r, "quadclassno");
5924
  if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5925
  Hf = conductor_part(x, r, &D, NULL);
5926
  return gerepileuptoint(av, mulii(Hf, gel(quadclassunit0(D,0,NULL,0),1)));
5927
}
5928

5929
/* use Euler products */
5930
GEN
5931
classno2(GEN x)
5932
{
5933
  pari_sp av = avma;
5934
  const long prec = DEFAULTPREC;
5935
  long n, i, r, s;
5936
  GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
5937

5938
  check_quaddisc(x, &s, &r, "classno2");
5939
  if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5940

5941
  Hf = conductor_part(x, r, &D, &reg);
5942
  if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
5943

5944
  Pi = mppi(prec);
5945
  d = absi_shallow(D); dr = itor(d, prec);
5946
  logd = logr_abs(dr);
5947
  p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
5948
  if (s > 0)
5949
  {
5950
    GEN invlogd = invr(logd);
5951
    p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
5952
    if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
5953
  }
5954
  n = itos_or_0( mptrunc(p1) );
5955
  if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
5956

5957
  p4 = divri(Pi,d);
5958
  p7 = invr(sqrtr_abs(Pi));
5959
  half = real2n(-1, prec);
5960
  if (s > 0)
5961
  { /* i = 1, shortcut */
5962
    p1 = sqrtr_abs(dr);
5963
    p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5964
    S = addrr(mulrr(p1,p5), eint1(p4,prec));
5965
    for (i=2; i<=n; i++)
5966
    {
5967
      long k = kroiu(D,i); if (!k) continue;
5968
      p2 = mulir(sqru(i), p4);
5969
      p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5970
      p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
5971
      S = (k>0)? addrr(S,p5): subrr(S,p5);
5972
    }
5973
    S = shiftr(divrr(S,reg),-1);
5974
  }
5975
  else
5976
  { /* i = 1, shortcut */
5977
    p1 = gdiv(sqrtr_abs(dr), Pi);
5978
    p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5979
    S = addrr(p5, divrr(p1, mpexp(p4)));
5980
    for (i=2; i<=n; i++)
5981
    {
5982
      long k = kroiu(D,i); if (!k) continue;
5983
      p2 = mulir(sqru(i), p4);
5984
      p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5985
      p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
5986
      S = (k>0)? addrr(S,p5): subrr(S,p5);
5987
    }
5988
  }
5989
  return gerepileuptoint(av, mulii(Hf, roundr(S)));
5990
}
5991

5992
/* 1 + q + ... + q^v, v > 0 */
5993
static GEN
5994
geomsumu(ulong q, long v)
5995
{
5996
  GEN u = utoipos(1+q);
5997
  for (; v > 1; v--) u = addui(1, mului(q, u));
5998
  return u;
5999
}
6000
static GEN
6001
geomsum(GEN q, long v)
6002
{
6003
  GEN u;
6004
  if (lgefint(q) == 3) return geomsumu(q[2], v);
6005
  u = addiu(q,1);
6006
  for (; v > 1; v--) u = addui(1, mulii(q, u));
6007
  return u;
6008
}
6009

6010
static GEN
6011
hclassno6_large(GEN x)
6012
{
6013
  long i, l, s, xmod4;
6014
  GEN H = NULL, D, P, E;
6015

6016
  x = negi(x);
6017
  check_quaddisc(x, &s, &xmod4, "hclassno");
6018
  corediscfact(x, xmod4, &D, &P, &E);
6019
  l = lg(P);
6020
  if (l > 1 && lgefint(x) == 3)
6021
  { /* F != 1, second chance */
6022
    ulong h = hclassno6u_from_cache(x[2]);
6023
    if (h) H = utoipos(h);
6024
  }
6025
  if (!H)
6026
  {
6027
    GEN Q = quadclassunit0(D, 0, NULL, 0);
6028
    H = gel(Q,1);
6029
    switch(itou_or_0(D))
6030
    {
6031
      case 3: H = shifti(H,1);break;
6032
      case 4: H = muliu(H,3); break;
6033
      default:H = muliu(H,6); break;
6034
    }
6035
  }
6036
  /* H \prod_{p^e||f}  (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
6037
  for (i = 1; i < l; i++)
6038
  {
6039
    long e = E[i], s;
6040
    GEN p, t;
6041
    if (!e) continue;
6042
    p = gel(P,i); s = kronecker(D,p);
6043
    if (e == 1) t = addiu(p, 1-s);
6044
    else if (s == 1) t = powiu(p, e);
6045
    else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
6046
    H = mulii(H,t);
6047
  }
6048
  return H;
6049
}
6050

6051
/* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
6052
GEN
6053
hclassno6(GEN x)
6054
{
6055
  ulong d = itou_or_0(x);
6056
  if (d)
6057
  { /* create cache if d small */
6058
    ulong h = d < 500000 ? hclassno6u(d): hclassno6u_from_cache(d);
6059
    if (h) return utoipos(h);
6060
  }
6061
  return hclassno6_large(x);
6062
}
6063

6064
GEN
6065
hclassno(GEN x)
6066
{
6067
  long a, s;
6068
  if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
6069
  s = signe(x);
6070
  if (s < 0) return gen_0;
6071
  if (!s) return gdivgs(gen_1, -12);
6072
  a = mod4(x); if (a == 1 || a == 2) return gen_0;
6073
  return gdivgs(hclassno6(x), 6);
6074
}
6075

6076