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