1 /*
2 Copyright (C) 2020 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <https://www.gnu.org/licenses/>.
10 */
11
12 #include "n_poly.h"
13 #include "mpn_extras.h"
14
15
n_poly_mod_is_canonical(const n_poly_t A,nmod_t mod)16 int n_poly_mod_is_canonical(const n_poly_t A, nmod_t mod)
17 {
18 slong i;
19 if (A->length < 0)
20 return 0;
21 for (i = 0; i < A->length; i++)
22 {
23 if (A->coeffs[i] >= mod.n)
24 return 0;
25 if (A->coeffs[i] == 0 && i + 1 == A->length)
26 return 0;
27 }
28 return 1;
29 }
30
31
n_poly_mod_set_coeff_ui(n_poly_t poly,slong j,ulong c,nmod_t ctx)32 void n_poly_mod_set_coeff_ui(
33 n_poly_t poly,
34 slong j,
35 ulong c,
36 nmod_t ctx)
37 {
38 if (c >= ctx.n)
39 NMOD_RED(c, c, ctx);
40
41 n_poly_set_coeff(poly, j, c);
42 }
43
44
n_poly_mod_add_ui(n_poly_t res,const n_poly_t poly,ulong c,nmod_t ctx)45 void n_poly_mod_add_ui(n_poly_t res, const n_poly_t poly, ulong c, nmod_t ctx)
46 {
47 if (c >= ctx.n)
48 NMOD_RED(c, c, ctx);
49
50 if (poly->length < 1)
51 {
52 n_poly_set_ui(res, c);
53 }
54 else
55 {
56 n_poly_set(res, poly);
57 res->coeffs[0] = nmod_add(res->coeffs[0], c, ctx);
58 _n_poly_normalise(res);
59 }
60 }
61
62
n_poly_mod_div_root(n_poly_t Q,const n_poly_t A,mp_limb_t c,nmod_t ctx)63 mp_limb_t n_poly_mod_div_root(n_poly_t Q,
64 const n_poly_t A, mp_limb_t c, nmod_t ctx)
65 {
66 mp_limb_t rem;
67
68 slong len = A->length;
69
70 if (len < 2)
71 {
72 if (len == 1)
73 {
74 rem = A->coeffs[0];
75 n_poly_zero(Q);
76 return rem;
77 }
78
79 n_poly_zero(Q);
80 return 0;
81 }
82
83 n_poly_fit_length(Q, len - 1);
84 rem = _nmod_poly_div_root(Q->coeffs, A->coeffs, len, c, ctx);
85 Q->length = len - 1;
86 return rem;
87 }
88
n_poly_mod_pow(n_poly_t res,const n_poly_t poly,ulong e,nmod_t ctx)89 void n_poly_mod_pow(n_poly_t res, const n_poly_t poly, ulong e, nmod_t ctx)
90 {
91 const slong len = poly->length;
92 slong rlen;
93
94 if ((len < 2) | (e < UWORD(3)))
95 {
96 if (len == 0)
97 {
98 if (e == 0)
99 n_poly_one(res);
100 else
101 n_poly_zero(res);
102 }
103 else if (len == 1)
104 {
105 n_poly_set_ui(res,
106 n_powmod2_ui_preinv(poly->coeffs[0], e, ctx.n, ctx.ninv));
107 }
108 else if (e == 0)
109 {
110 n_poly_one(res);
111 }
112 else if (e == 1)
113 n_poly_set(res, poly);
114 else /* e == UWORD(2) */
115 n_poly_mod_mul(res, poly, poly, ctx);
116
117 return;
118 }
119
120 rlen = (slong) e * (len - 1) + 1;
121
122 if (res != poly)
123 {
124 n_poly_fit_length(res, rlen);
125 _nmod_poly_pow(res->coeffs, poly->coeffs, len, e, ctx);
126 }
127 else
128 {
129 n_poly_t t;
130 n_poly_init2(t, rlen);
131 _nmod_poly_pow(t->coeffs, poly->coeffs, len, e, ctx);
132 n_poly_swap(res, t);
133 n_poly_clear(t);
134 }
135
136 res->length = rlen;
137 _n_poly_normalise(res);
138 }
139
140
n_poly_mod_mul(n_poly_t res,const n_poly_t poly1,const n_poly_t poly2,nmod_t ctx)141 void n_poly_mod_mul(n_poly_t res, const n_poly_t poly1,
142 const n_poly_t poly2, nmod_t ctx)
143 {
144 slong len1, len2, len_out;
145
146 len1 = poly1->length;
147 len2 = poly2->length;
148
149 if (len1 == 0 || len2 == 0)
150 {
151 n_poly_zero(res);
152
153 return;
154 }
155
156 len_out = poly1->length + poly2->length - 1;
157
158 if (res == poly1 || res == poly2)
159 {
160 n_poly_t temp;
161
162 n_poly_init2(temp, len_out);
163
164 if (len1 >= len2)
165 _nmod_poly_mul(temp->coeffs, poly1->coeffs, len1,
166 poly2->coeffs, len2, ctx);
167 else
168 _nmod_poly_mul(temp->coeffs, poly2->coeffs, len2,
169 poly1->coeffs, len1, ctx);
170
171 n_poly_swap(temp, res);
172 n_poly_clear(temp);
173 }
174 else
175 {
176 n_poly_fit_length(res, len_out);
177
178 if (len1 >= len2)
179 _nmod_poly_mul(res->coeffs, poly1->coeffs, len1,
180 poly2->coeffs, len2, ctx);
181 else
182 _nmod_poly_mul(res->coeffs, poly2->coeffs, len2,
183 poly1->coeffs, len1, ctx);
184 }
185
186 res->length = len_out;
187 _n_poly_normalise(res);
188 }
189
n_poly_mod_mullow(n_poly_t res,const n_poly_t poly1,const n_poly_t poly2,slong trunc,nmod_t ctx)190 void n_poly_mod_mullow(
191 n_poly_t res,
192 const n_poly_t poly1,
193 const n_poly_t poly2,
194 slong trunc,
195 nmod_t ctx)
196 {
197 slong len1, len2, len_out;
198
199 len1 = poly1->length;
200 len2 = poly2->length;
201
202 len_out = poly1->length + poly2->length - 1;
203 if (trunc > len_out)
204 trunc = len_out;
205
206 if (len1 <= 0 || len2 <= 0 || trunc <= 0)
207 {
208 n_poly_zero(res);
209 return;
210 }
211
212 if (res == poly1 || res == poly2)
213 {
214 n_poly_t temp;
215
216 n_poly_init2(temp, trunc);
217
218 if (len1 >= len2)
219 _nmod_poly_mullow(temp->coeffs, poly1->coeffs, len1,
220 poly2->coeffs, len2, trunc, ctx);
221 else
222 _nmod_poly_mullow(temp->coeffs, poly2->coeffs, len2,
223 poly1->coeffs, len1, trunc, ctx);
224
225 n_poly_swap(temp, res);
226 n_poly_clear(temp);
227 }
228 else
229 {
230 n_poly_fit_length(res, trunc);
231
232 if (len1 >= len2)
233 _nmod_poly_mullow(res->coeffs, poly1->coeffs, len1,
234 poly2->coeffs, len2, trunc, ctx);
235 else
236 _nmod_poly_mullow(res->coeffs, poly2->coeffs, len2,
237 poly1->coeffs, len1, trunc, ctx);
238 }
239
240 res->length = trunc;
241 _n_poly_normalise(res);
242 }
243
244
n_poly_mod_mulmod(n_poly_t res,const n_poly_t poly1,const n_poly_t poly2,const n_poly_t f,nmod_t ctx)245 void n_poly_mod_mulmod(n_poly_t res, const n_poly_t poly1,
246 const n_poly_t poly2, const n_poly_t f, nmod_t ctx)
247 {
248 slong len1, len2, lenf;
249 mp_ptr fcoeffs;
250
251 lenf = f->length;
252 len1 = poly1->length;
253 len2 = poly2->length;
254
255 if (lenf == 0)
256 {
257 flint_printf("Exception (nmod_poly_mulmod). Divide by zero.\n");
258 flint_abort();
259 }
260
261 if (lenf == 1 || len1 == 0 || len2 == 0)
262 {
263 n_poly_zero(res);
264 return;
265 }
266
267 if (len1 + len2 - lenf > 0)
268 {
269 if (f == res)
270 {
271 fcoeffs = flint_malloc(sizeof(mp_limb_t) * lenf);
272 _nmod_vec_set(fcoeffs, f->coeffs, lenf);
273 }
274 else
275 fcoeffs = f->coeffs;
276
277 n_poly_fit_length(res, lenf - 1);
278 _nmod_poly_mulmod(res->coeffs, poly1->coeffs, len1,
279 poly2->coeffs, len2,
280 fcoeffs, lenf,
281 ctx);
282 if (f == res)
283 flint_free(fcoeffs);
284
285 res->length = lenf - 1;
286 _n_poly_normalise(res);
287 }
288 else
289 {
290 n_poly_mod_mul(res, poly1, poly2, ctx);
291 }
292 }
293
294
n_poly_mod_div(n_poly_t Q,const n_poly_t A,const n_poly_t B,nmod_t ctx)295 void n_poly_mod_div(n_poly_t Q, const n_poly_t A, const n_poly_t B, nmod_t ctx)
296 {
297 n_poly_t tQ;
298 mp_ptr q;
299 slong A_len, B_len;
300
301 B_len = B->length;
302
303 if (B_len == 0)
304 {
305 if (ctx.n == 1)
306 {
307 n_poly_set(Q, A);
308 return;
309 }
310 else
311 {
312 flint_printf("Exception (n_poly_mod_div). Division by zero.\n");
313 flint_abort();
314 }
315 }
316
317 A_len = A->length;
318
319 if (A_len < B_len)
320 {
321 n_poly_zero(Q);
322 return;
323 }
324
325 if (Q == A || Q == B)
326 {
327 n_poly_init2(tQ, A_len - B_len + 1);
328 q = tQ->coeffs;
329 }
330 else
331 {
332 n_poly_fit_length(Q, A_len - B_len + 1);
333 q = Q->coeffs;
334 }
335
336 _nmod_poly_div(q, A->coeffs, A_len, B->coeffs, B_len, ctx);
337
338 if (Q == A || Q == B)
339 {
340 n_poly_swap(tQ, Q);
341 n_poly_clear(tQ);
342 }
343
344 Q->length = A_len - B_len + 1;
345 }
346
n_poly_mod_rem(n_poly_t R,const n_poly_t A,const n_poly_t B,nmod_t ctx)347 void n_poly_mod_rem(n_poly_t R, const n_poly_t A, const n_poly_t B, nmod_t ctx)
348 {
349 const slong lenA = A->length, lenB = B->length;
350 n_poly_t tR;
351 mp_ptr r;
352
353 if (lenB == 0)
354 {
355 flint_printf("Exception (nmod_poly_rem). Division by zero.\n");
356 flint_abort();
357 }
358 if (lenA < lenB)
359 {
360 n_poly_set(R, A);
361 return;
362 }
363
364 if (R == A || R == B)
365 {
366 n_poly_init2(tR, lenB - 1);
367 r = tR->coeffs;
368 }
369 else
370 {
371 n_poly_fit_length(R, lenB - 1);
372 r = R->coeffs;
373 }
374
375 _nmod_poly_rem(r, A->coeffs, lenA, B->coeffs, lenB, ctx);
376
377 if (R == A || R == B)
378 {
379 n_poly_swap(R, tR);
380 n_poly_clear(tR);
381 }
382
383 R->length = lenB - 1;
384 _n_poly_normalise(R);
385 }
386
387
n_poly_mod_divrem(n_poly_t Q,n_poly_t R,const n_poly_t A,const n_poly_t B,nmod_t ctx)388 void n_poly_mod_divrem(n_poly_t Q, n_poly_t R,
389 const n_poly_t A, const n_poly_t B, nmod_t ctx)
390 {
391 const slong lenA = A->length, lenB = B->length;
392 n_poly_t tQ, tR;
393 mp_ptr q, r;
394
395 if (lenB == 0)
396 {
397 if (ctx.n == 1)
398 {
399 n_poly_set(Q, A);
400 n_poly_zero(R);
401 return;
402 }
403 else
404 {
405 flint_printf("Exception (n_poly_mod_divrem). Division by zero.");
406 flint_abort();
407 }
408 }
409
410 if (lenA < lenB)
411 {
412 n_poly_set(R, A);
413 n_poly_zero(Q);
414 return;
415 }
416
417 if (Q == A || Q == B)
418 {
419 n_poly_init2(tQ, lenA - lenB + 1);
420 q = tQ->coeffs;
421 }
422 else
423 {
424 n_poly_fit_length(Q, lenA - lenB + 1);
425 q = Q->coeffs;
426 }
427
428 if (R == A || R == B)
429 {
430 n_poly_fit_length(tR, lenB - 1);
431 r = tR->coeffs;
432 }
433 else
434 {
435 n_poly_fit_length(R, lenB - 1);
436 r = R->coeffs;
437 }
438
439 _nmod_poly_divrem(q, r, A->coeffs, lenA, B->coeffs, lenB, ctx);
440
441 if (Q == A || Q == B)
442 {
443 n_poly_swap(Q, tQ);
444 n_poly_clear(tQ);
445 }
446 if (R == A || R == B)
447 {
448 n_poly_swap(R, tR);
449 n_poly_clear(tR);
450 }
451
452 Q->length = lenA - lenB + 1;
453 R->length = lenB - 1;
454
455 _n_poly_normalise(R);
456 }
457
458
n_poly_mod_invmod(n_poly_t A,const n_poly_t B,const n_poly_t P,nmod_t ctx)459 int n_poly_mod_invmod(n_poly_t A, const n_poly_t B, const n_poly_t P, nmod_t ctx)
460 {
461 const slong lenB = B->length, lenP = P->length;
462 mp_limb_t * a;
463 n_poly_t tA;
464 int ans;
465
466 if (lenP < 2)
467 {
468 printf("Exception (nmod_poly_invmod). lenP < 2.\n");
469 flint_abort();
470 }
471 if (lenB == 0)
472 {
473 n_poly_zero(A);
474 return 0;
475 }
476 if (lenB >= lenP)
477 {
478 n_poly_t T;
479
480 n_poly_init(T);
481 n_poly_mod_rem(T, B, P, ctx);
482 ans = n_poly_mod_invmod(A, T, P, ctx);
483 n_poly_clear(T);
484 return ans;
485 }
486
487 if (A != B && A != P)
488 {
489 n_poly_fit_length(A, lenP - 1);
490 a = A->coeffs;
491 }
492 else
493 {
494 n_poly_init2(tA, lenP - 1);
495 a = tA->coeffs;
496 }
497
498 ans = _nmod_poly_invmod(a, B->coeffs, lenB, P->coeffs, lenP, ctx);
499
500 if (A == B || A == P)
501 {
502 n_poly_swap(A, tA);
503 n_poly_clear(tA);
504 }
505
506 A->length = lenP - 1;
507 _n_poly_normalise(A);
508 return ans;
509 }
510
511
n_poly_mod_gcd(n_poly_t G,const n_poly_t A,const n_poly_t B,nmod_t ctx)512 void n_poly_mod_gcd(n_poly_t G, const n_poly_t A, const n_poly_t B, nmod_t ctx)
513 {
514 if (A->length < B->length)
515 {
516 n_poly_mod_gcd(G, B, A, ctx);
517 }
518 else /* lenA >= lenB >= 0 */
519 {
520 slong lenA = A->length, lenB = B->length, lenG;
521 n_poly_t tG;
522 mp_ptr g;
523
524 if (lenA == 0) /* lenA = lenB = 0 */
525 {
526 n_poly_zero(G);
527 }
528 else if (lenB == 0) /* lenA > lenB = 0 */
529 {
530 n_poly_mod_make_monic(G, A, ctx);
531 }
532 else /* lenA >= lenB >= 1 */
533 {
534 if (G == A || G == B)
535 {
536 n_poly_init2(tG, FLINT_MIN(lenA, lenB));
537 g = tG->coeffs;
538 }
539 else
540 {
541 n_poly_fit_length(G, FLINT_MIN(lenA, lenB));
542 g = G->coeffs;
543 }
544
545 lenG = _nmod_poly_gcd(g, A->coeffs, lenA, B->coeffs, lenB, ctx);
546
547 if (G == A || G == B)
548 {
549 n_poly_swap(tG, G);
550 n_poly_clear(tG);
551 }
552 G->length = lenG;
553
554 if (G->length == 1)
555 G->coeffs[0] = 1;
556 else
557 n_poly_mod_make_monic(G, G, ctx);
558 }
559 }
560 }
561
n_poly_mod_xgcd(n_poly_t G,n_poly_t S,n_poly_t T,const n_poly_t A,const n_poly_t B,nmod_t ctx)562 void n_poly_mod_xgcd(
563 n_poly_t G,
564 n_poly_t S,
565 n_poly_t T,
566 const n_poly_t A,
567 const n_poly_t B,
568 nmod_t ctx)
569 {
570 if (A->length < B->length)
571 {
572 n_poly_mod_xgcd(G, T, S, B, A, ctx);
573 }
574 else /* lenA >= lenB >= 0 */
575 {
576 const slong lenA = A->length, lenB = B->length;
577 mp_limb_t inv;
578
579 if (lenA == 0) /* lenA = lenB = 0 */
580 {
581 n_poly_zero(G);
582 n_poly_zero(S);
583 n_poly_zero(T);
584 }
585 else if (lenB == 0) /* lenA > lenB = 0 */
586 {
587 inv = n_invmod(A->coeffs[lenA - 1], ctx.n);
588 _n_poly_mod_scalar_mul_nmod(G, A, inv, ctx);
589 n_poly_zero(T);
590 n_poly_set_coeff(S, 0, inv);
591 S->length = 1;
592 }
593 else if (lenB == 1) /* lenA >= lenB = 1 */
594 {
595 n_poly_fit_length(T, 1);
596 T->length = 1;
597 T->coeffs[0] = n_invmod(B->coeffs[0], ctx.n);
598 n_poly_one(G);
599 n_poly_zero(S);
600 }
601 else /* lenA >= lenB >= 2 */
602 {
603 mp_ptr g, s, t;
604 slong lenG;
605
606 if (G == A || G == B)
607 {
608 g = _nmod_vec_init(FLINT_MIN(lenA, lenB));
609 }
610 else
611 {
612 n_poly_fit_length(G, FLINT_MIN(lenA, lenB));
613 g = G->coeffs;
614 }
615 if (S == A || S == B)
616 {
617 s = _nmod_vec_init(lenB - 1);
618 }
619 else
620 {
621 n_poly_fit_length(S, lenB - 1);
622 s = S->coeffs;
623 }
624 if (T == A || T == B)
625 {
626 t = _nmod_vec_init(lenA - 1);
627 }
628 else
629 {
630 n_poly_fit_length(T, lenA - 1);
631 t = T->coeffs;
632 }
633
634 if (lenA >= lenB)
635 lenG = _nmod_poly_xgcd(g, s, t, A->coeffs, lenA,
636 B->coeffs, lenB, ctx);
637 else
638 lenG = _nmod_poly_xgcd(g, t, s, B->coeffs, lenB,
639 A->coeffs, lenA, ctx);
640
641 if (G == A || G == B)
642 {
643 flint_free(G->coeffs);
644 G->coeffs = g;
645 G->alloc = FLINT_MIN(lenA, lenB);
646 }
647 if (S == A || S == B)
648 {
649 flint_free(S->coeffs);
650 S->coeffs = s;
651 S->alloc = lenB - 1;
652 }
653 if (T == A || T == B)
654 {
655 flint_free(T->coeffs);
656 T->coeffs = t;
657 T->alloc = lenA - 1;
658 }
659
660 G->length = lenG;
661 S->length = FLINT_MAX(lenB - lenG, 1);
662 T->length = FLINT_MAX(lenA - lenG, 1);
663 MPN_NORM(S->coeffs, S->length);
664 MPN_NORM(T->coeffs, T->length);
665
666 if (G->coeffs[lenG - 1] != 1)
667 {
668 inv = nmod_inv(G->coeffs[lenG - 1], ctx);
669 _n_poly_mod_scalar_mul_nmod(G, G, inv, ctx);
670 _n_poly_mod_scalar_mul_nmod(S, S, inv, ctx);
671 _n_poly_mod_scalar_mul_nmod(T, T, inv, ctx);
672 }
673 }
674 }
675 }
676
677
n_poly_reverse(n_poly_t output,const n_poly_t input,slong m)678 void n_poly_reverse(n_poly_t output, const n_poly_t input, slong m)
679 {
680 n_poly_fit_length(output, m);
681 _nmod_poly_reverse(output->coeffs, input->coeffs, input->length, m);
682 output->length = m;
683 _n_poly_normalise(output);
684 }
685
686
n_poly_mod_mulmod_preinv(n_poly_t res,const n_poly_t poly1,const n_poly_t poly2,const n_poly_t f,const n_poly_t finv,nmod_t ctx)687 void n_poly_mod_mulmod_preinv(
688 n_poly_t res,
689 const n_poly_t poly1,
690 const n_poly_t poly2,
691 const n_poly_t f,
692 const n_poly_t finv,
693 nmod_t ctx)
694 {
695 slong len1, len2, lenf;
696 mp_ptr fcoeffs;
697
698 lenf = f->length;
699 len1 = poly1->length;
700 len2 = poly2->length;
701
702 if (lenf <= len1 || lenf <= len2)
703 {
704 flint_printf("n_poly_mod_mulmod_preinv: Input is larger than modulus.");
705 flint_abort();
706 }
707
708 if (lenf == 1 || len1 == 0 || len2 == 0)
709 {
710 n_poly_zero(res);
711 return;
712 }
713
714 if (len1 + len2 > lenf)
715 {
716 if (f == res)
717 {
718 fcoeffs = flint_malloc(sizeof(mp_limb_t) * lenf);
719 _nmod_vec_set(fcoeffs, f->coeffs, lenf);
720 }
721 else
722 {
723 fcoeffs = f->coeffs;
724 }
725
726 n_poly_fit_length(res, lenf - 1);
727 _nmod_poly_mulmod_preinv(res->coeffs, poly1->coeffs, len1,
728 poly2->coeffs, len2,
729 fcoeffs, lenf,
730 finv->coeffs, finv->length, ctx);
731 if (f == res)
732 flint_free(fcoeffs);
733
734 res->length = lenf - 1;
735 _n_poly_normalise(res);
736 }
737 else
738 {
739 n_poly_mod_mul(res, poly1, poly2, ctx);
740 }
741 }
742
n_poly_mod_inv_series(n_poly_t Qinv,const n_poly_t Q,slong n,nmod_t ctx)743 void n_poly_mod_inv_series(n_poly_t Qinv, const n_poly_t Q, slong n, nmod_t ctx)
744 {
745 slong Qlen = Q->length;
746
747 Qlen = FLINT_MIN(Qlen, n);
748
749 if (Qlen == 0)
750 {
751 flint_throw(FLINT_ERROR, "n_poly_mod_inv_series_newton: Division by zero.");
752 }
753
754 if (Qinv != Q)
755 {
756 n_poly_fit_length(Qinv, n);
757 _nmod_poly_inv_series_newton(Qinv->coeffs, Q->coeffs, Qlen, n, ctx);
758 }
759 else
760 {
761 n_poly_t t;
762 n_poly_init2(t, n);
763 _nmod_poly_inv_series_newton(t->coeffs, Q->coeffs, Qlen, n, ctx);
764 n_poly_swap(Qinv, t);
765 n_poly_clear(t);
766 }
767
768 Qinv->length = n;
769 _n_poly_normalise(Qinv);
770 }
771
772
n_poly_mod_div_series(n_poly_t Q,const n_poly_t A,const n_poly_t B,slong order,nmod_t ctx)773 void n_poly_mod_div_series(n_poly_t Q, const n_poly_t A, const n_poly_t B,
774 slong order, nmod_t ctx)
775 {
776 slong Blen = B->length;
777 slong Alen = A->length;
778
779 if (order < 1 || Blen == 0 || B->coeffs[0] == 0)
780 {
781 flint_printf("Exception (n_poly_div_series). Division by zero.\n");
782 flint_abort();
783 }
784
785 if (Alen == 0)
786 {
787 n_poly_zero(Q);
788 return;
789 }
790
791 if (Q != A && Q != B)
792 {
793 n_poly_fit_length(Q, order);
794 _nmod_poly_div_series(Q->coeffs, A->coeffs, Alen, B->coeffs, Blen, order, ctx);
795 }
796 else
797 {
798 n_poly_t t;
799 n_poly_init(t);
800 _nmod_poly_div_series(t->coeffs, A->coeffs, Alen, B->coeffs, Blen, order, ctx);
801 n_poly_swap(Q, t);
802 n_poly_clear(t);
803 }
804
805 Q->length = order;
806 _n_poly_normalise(Q);
807 }
808
n_poly_mod_scalar_mul_ui(n_poly_t A,const n_poly_t B,mp_limb_t c,nmod_t ctx)809 void n_poly_mod_scalar_mul_ui(n_poly_t A, const n_poly_t B, mp_limb_t c, nmod_t ctx)
810 {
811 if (c >= ctx.n)
812 {
813 NMOD_RED(c, c, ctx);
814 }
815
816 if (c == 0 || B->length < 1)
817 {
818 n_poly_zero(A);
819 return;
820 }
821
822 _n_poly_mod_scalar_mul_nmod(A, B, c, ctx);
823 _n_poly_normalise(A);
824 }
825
826 /* multiply A by (x^k + c) */
n_poly_mod_shift_left_scalar_addmul(n_poly_t A,slong k,mp_limb_t c,nmod_t ctx)827 void n_poly_mod_shift_left_scalar_addmul(n_poly_t A, slong k, mp_limb_t c,
828 nmod_t ctx)
829 {
830 mp_limb_t * Acoeffs;
831 slong i;
832 slong Alen = A->length;
833
834 n_poly_fit_length(A, Alen + k);
835
836 Acoeffs = A->coeffs;
837
838 flint_mpn_copyd(Acoeffs + k, Acoeffs, Alen);
839 flint_mpn_zero(Acoeffs, k);
840
841 for (i = 0; i < A->length; i++)
842 Acoeffs[i] = nmod_addmul(Acoeffs[i], c, Acoeffs[i + k], ctx);
843
844 A->length = Alen + k;
845 }
846
847 /* A = B + C*(d1*x+d0) */
n_poly_mod_addmul_linear(n_poly_t A,const n_poly_t B,const n_poly_t C,mp_limb_t d1,mp_limb_t d0,nmod_t ctx)848 void n_poly_mod_addmul_linear(
849 n_poly_t A,
850 const n_poly_t B,
851 const n_poly_t C,
852 mp_limb_t d1, mp_limb_t d0,
853 nmod_t ctx)
854 {
855 slong i;
856 mp_limb_t * Acoeffs, * Bcoeffs, * Ccoeffs;
857 slong Blen = B->length;
858 slong Clen = C->length;
859 slong Alen = FLINT_MAX(B->length, C->length + 1);
860
861 FLINT_ASSERT(A != B);
862 FLINT_ASSERT(A != C);
863
864 n_poly_fit_length(A, Alen);
865 Acoeffs = A->coeffs;
866 Bcoeffs = B->coeffs;
867 Ccoeffs = C->coeffs;
868
869 for (i = 0; i < Alen; i++)
870 {
871 ulong p1, p0, t0 = 0, t1 = 0, t2 = 0;
872
873 if (i < Blen)
874 {
875 t0 = Bcoeffs[i];
876 }
877 if (i < Clen)
878 {
879 umul_ppmm(p1, p0, Ccoeffs[i], d0);
880 add_ssaaaa(t1, t0, t1, t0, p1, p0);
881 }
882 if (0 < i && i - 1 < Clen)
883 {
884 umul_ppmm(p1, p0, Ccoeffs[i - 1], d1);
885 add_sssaaaaaa(t2, t1, t0, t2, t1, t0, 0, p1, p0);
886 }
887 NMOD_RED3(Acoeffs[i], t2, t1, t0, ctx);
888 }
889
890 A->length = Alen;
891 _n_poly_normalise(A);
892 }
893
894 /* A = B + C*d0 */
n_poly_mod_scalar_addmul_nmod(n_poly_t A,const n_poly_t B,const n_poly_t C,mp_limb_t d0,nmod_t ctx)895 void n_poly_mod_scalar_addmul_nmod(
896 n_poly_t A,
897 const n_poly_t B,
898 const n_poly_t C,
899 mp_limb_t d0,
900 nmod_t ctx)
901 {
902 slong i;
903 mp_limb_t t0, t1;
904 mp_limb_t * Acoeffs, * Bcoeffs, * Ccoeffs;
905 slong Blen = B->length;
906 slong Clen = C->length;
907 slong Alen = FLINT_MAX(B->length, C->length);
908
909 n_poly_fit_length(A, Alen);
910 Acoeffs = A->coeffs;
911 Bcoeffs = B->coeffs;
912 Ccoeffs = C->coeffs;
913
914 if (ctx.norm >= (FLINT_BITS + 1)/2)
915 {
916 for (i = 0; i + 2 <= FLINT_MIN(Blen, Clen); i += 2)
917 {
918 NMOD_RED2(t0, 0, Bcoeffs[i + 0] + d0*Ccoeffs[i + 0], ctx);
919 NMOD_RED2(t1, 0, Bcoeffs[i + 1] + d0*Ccoeffs[i + 1], ctx);
920 Acoeffs[i + 0] = t0;
921 Acoeffs[i + 1] = t1;
922 }
923 for ( ; i < FLINT_MIN(Blen, Clen); i++)
924 {
925 NMOD_RED2(Acoeffs[i], 0, Bcoeffs[i] + d0*Ccoeffs[i], ctx);
926 }
927
928 for ( ; i + 2 <= Clen; i += 2)
929 {
930 NMOD_RED2(t0, 0, d0*Ccoeffs[i + 0], ctx);
931 NMOD_RED2(t1, 0, d0*Ccoeffs[i + 1], ctx);
932 Acoeffs[i + 0] = t0;
933 Acoeffs[i + 1] = t1;
934 }
935 for ( ; i < Clen; i++)
936 {
937 NMOD_RED2(Acoeffs[i], 0, d0*Ccoeffs[i], ctx);
938 }
939 }
940 else
941 {
942 for (i = 0; i < FLINT_MIN(Blen, Clen); i++)
943 {
944 t0 = Bcoeffs[i];
945 NMOD_ADDMUL(t0, d0, Ccoeffs[i], ctx);
946 Acoeffs[i] = t0;
947 }
948
949 while (i < Clen)
950 {
951 Acoeffs[i] = nmod_mul(d0, Ccoeffs[i], ctx);
952 i++;
953 }
954 }
955
956 while (i < Blen)
957 {
958 Acoeffs[i] = Bcoeffs[i];
959 i++;
960 }
961
962 A->length = Alen;
963 _n_poly_normalise(A);
964 }
965
966
n_poly_mod_remove(n_poly_t f,const n_poly_t p,nmod_t ctx)967 ulong n_poly_mod_remove(n_poly_t f, const n_poly_t p, nmod_t ctx)
968 {
969 n_poly_t q, r;
970 ulong i = 0;
971
972 n_poly_init(q);
973 n_poly_init(r);
974
975 while (1)
976 {
977 if (f->length < p->length)
978 break;
979 n_poly_mod_divrem(q, r, f, p, ctx);
980 if (r->length == 0)
981 n_poly_swap(q, f);
982 else
983 break;
984 i++;
985 }
986
987 n_poly_clear(q);
988 n_poly_clear(r);
989
990 return i;
991 }
992
_n_poly_eval_pow(n_poly_t P,n_poly_t alphapow,int nlimbs,nmod_t ctx)993 mp_limb_t _n_poly_eval_pow(n_poly_t P, n_poly_t alphapow, int nlimbs, nmod_t ctx)
994 {
995 mp_limb_t * Pcoeffs = P->coeffs;
996 slong Plen = P->length;
997 mp_limb_t * alpha_powers = alphapow->coeffs;
998 mp_limb_t res;
999 slong k;
1000
1001 if (Plen > alphapow->length)
1002 {
1003 slong oldlength = alphapow->length;
1004 FLINT_ASSERT(2 <= oldlength);
1005 n_poly_fit_length(alphapow, Plen);
1006 alphapow->length = Plen;
1007 alpha_powers = alphapow->coeffs;
1008 for (k = oldlength; k < Plen; k++)
1009 alpha_powers[k] = nmod_mul(alpha_powers[k - 1], alpha_powers[1], ctx);
1010 }
1011
1012 NMOD_VEC_DOT(res, k, Plen, Pcoeffs[k], alpha_powers[k], ctx, nlimbs);
1013
1014 return res;
1015 }
1016
n_poly_mod_eval_pow(n_poly_t P,n_poly_t alphapow,nmod_t ctx)1017 mp_limb_t n_poly_mod_eval_pow(n_poly_t P, n_poly_t alphapow, nmod_t ctx)
1018 {
1019 int nlimbs = _nmod_vec_dot_bound_limbs(P->length, ctx);
1020 return _n_poly_eval_pow(P, alphapow, nlimbs, ctx);
1021 }
1022
n_poly_mod_eval2_pow(mp_limb_t * vp,mp_limb_t * vm,const n_poly_t P,n_poly_t alphapow,nmod_t ctx)1023 void n_poly_mod_eval2_pow(
1024 mp_limb_t * vp,
1025 mp_limb_t * vm,
1026 const n_poly_t P,
1027 n_poly_t alphapow,
1028 nmod_t ctx)
1029 {
1030 const mp_limb_t * Pcoeffs = P->coeffs;
1031 slong Plen = P->length;
1032 mp_limb_t * alpha_powers = alphapow->coeffs;
1033 mp_limb_t p1, p0, a0, a1, a2, q1, q0, b0, b1, b2;
1034 slong k;
1035
1036 a0 = a1 = a2 = 0;
1037 b0 = b1 = b2 = 0;
1038
1039 if (Plen > alphapow->length)
1040 {
1041 slong oldlength = alphapow->length;
1042 FLINT_ASSERT(2 <= oldlength);
1043 n_poly_fit_length(alphapow, Plen);
1044 for (k = oldlength; k < Plen; k++)
1045 {
1046 alphapow->coeffs[k] = nmod_mul(alphapow->coeffs[k - 1],
1047 alphapow->coeffs[1], ctx);
1048 }
1049 alphapow->length = Plen;
1050 alpha_powers = alphapow->coeffs;
1051 }
1052
1053 for (k = 0; k + 2 <= Plen; k += 2)
1054 {
1055 umul_ppmm(p1, p0, Pcoeffs[k + 0], alpha_powers[k + 0]);
1056 umul_ppmm(q1, q0, Pcoeffs[k + 1], alpha_powers[k + 1]);
1057 add_sssaaaaaa(a2, a1, a0, a2, a1, a0, 0, p1, p0);
1058 add_sssaaaaaa(b2, b1, b0, b2, b1, b0, 0, q1, q0);
1059 }
1060
1061 if (k < Plen)
1062 {
1063 umul_ppmm(p1, p0, Pcoeffs[k + 0], alpha_powers[k + 0]);
1064 add_sssaaaaaa(a2, a1, a0, a2, a1, a0, 0, p1, p0);
1065 k++;
1066 }
1067
1068 FLINT_ASSERT(k == Plen);
1069
1070 NMOD_RED3(p0, a2, a1, a0, ctx);
1071 NMOD_RED3(q0, b2, b1, b0, ctx);
1072
1073 *vp = nmod_add(p0, q0, ctx);
1074 *vm = nmod_sub(p0, q0, ctx);
1075 }
1076
n_poly_mod_eval_step2(n_poly_t Acur,const n_poly_t Ainc,nmod_t mod)1077 mp_limb_t n_poly_mod_eval_step2(
1078 n_poly_t Acur,
1079 const n_poly_t Ainc,
1080 nmod_t mod)
1081 {
1082 slong i, Alen = Acur->length;
1083 mp_limb_t * cur = Acur->coeffs;
1084 const mp_limb_t * inc = Ainc->coeffs;
1085 ulong t0, t1, t2, p0, p1;
1086
1087 FLINT_ASSERT(2*Alen == Ainc->length);
1088
1089 t2 = t1 = t0 = 0;
1090 for (i = 0; i < Alen; i++)
1091 {
1092 umul_ppmm(p1, p0, cur[i], inc[2*i + 0]);
1093 add_sssaaaaaa(t2, t1, t0, t2, t1, t0, 0, p1, p0);
1094 cur[i] = nmod_mul(cur[i], inc[2*i + 1], mod);
1095 }
1096 NMOD_RED3(t0, t2, t1, t0, mod);
1097 return t0;
1098 }
1099
1100