Skip to content

Commit 9946b40

Browse files
fredrik-johanssonfredrik-johansson
committed
Generalize gr_poly mullow algorithms to mulmid
1 parent 17e900d commit 9946b40

File tree

12 files changed

+488
-55
lines changed

12 files changed

+488
-55
lines changed

‎doc/source/gr_poly.rst‎

Lines changed: 13 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -183,30 +183,34 @@ Arithmetic
183183
by default, which currently always delegates to :func:`_gr_poly_mullow_classical`.
184184
This can be overridden by specific rings.
185185

186-
.. function:: int _gr_poly_mulmid_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
187-
int gr_poly_mulmid_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
188-
int _gr_poly_mulmid_generic(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
186+
.. function:: int _gr_poly_mulmid_generic(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
189187
int _gr_poly_mulmid(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
190188
int gr_poly_mulmid(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, gr_ctx_t ctx)
191189

192190
Multiplication algorithms
193191
-------------------------------------------------------------------------------
194192

195-
.. function:: int _gr_poly_mullow_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
193+
.. function:: int _gr_poly_mulmid_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
194+
int gr_poly_mulmid_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
195+
int _gr_poly_mullow_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
196196
int gr_poly_mullow_classical(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
197197

198198
Multiply using the classical (schoolbook) algorithm, performing
199199
a sequence of dot products.
200200

201-
.. function:: int _gr_poly_mullow_bivariate_KS(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
201+
.. function:: int _gr_poly_mulmid_bivariate_KS(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
202+
int gr_poly_mulmid_bivariate_KS(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
203+
int _gr_poly_mullow_bivariate_KS(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
202204
int gr_poly_mullow_bivariate_KS(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
203205

204206
Assuming that the coefficients are polynomials of type ``gr_poly``,
205207
reduce the bivariate product to a single univariate product using
206208
Kronecker substitution. Returns ``GR_UNABLE`` if the elements are not
207209
of type ``gr_poly``.
208210

209-
.. function:: int _gr_poly_mullow_complex_reorder(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx)
211+
.. function:: int _gr_poly_mulmid_complex_reorder(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx)
212+
int gr_poly_mulmid_complex_reorder(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx)
213+
int _gr_poly_mullow_complex_reorder(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx)
210214
int gr_poly_mullow_complex_reorder(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx)
211215

212216
Assuming that the coefficients of the polynomials of type ``ctx`` are
@@ -270,7 +274,6 @@ Multiplication algorithms
270274
the interpolation matrix is computed with a common denominator and this
271275
denominator is removed from the final result with an exact division.
272276

273-
274277
The suffix "serial" indicates that the evaluations are done one point at a time
275278
and that the output is constructed incrementally with the goal to minimize
276279
memory allocation. This implementation is therefore well suited for huge
@@ -282,6 +285,9 @@ Multiplication algorithms
282285

283286
The underscore method requires positive lengths and does not support aliasing.
284287

288+
Scalar arithmetic
289+
-------------------------------------------------------------------------------
290+
285291
.. function:: int gr_poly_add_scalar(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
286292
int gr_poly_add_ui(gr_poly_t res, const gr_poly_t poly, ulong c, gr_ctx_t ctx)
287293
int gr_poly_add_si(gr_poly_t res, const gr_poly_t poly, slong c, gr_ctx_t ctx)

‎src/gr/fmpz_poly.c‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -826,6 +826,16 @@ _gr_fmpz_poly_gr_poly_mullow(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr
826826
return _gr_poly_mullow_bivariate_KS(res, poly1, len1, poly2, len2, n, ctx);
827827
}
828828

829+
static int
830+
_gr_fmpz_poly_gr_poly_mulmid(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
831+
{
832+
if (len1 < MUL_KS_CUTOFF || len2 < MUL_KS_CUTOFF || nhi < MUL_KS_CUTOFF
833+
|| len1 + len2 - 1 - nlo < MUL_KS_CUTOFF || 2 * (nhi - nlo) < MUL_KS_CUTOFF)
834+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
835+
else
836+
return _gr_poly_mulmid_bivariate_KS(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
837+
}
838+
829839

830840
int _fmpz_poly_methods_initialized = 0;
831841

@@ -923,6 +933,7 @@ gr_method_tab_input _fmpz_poly_methods_input[] =
923933
{GR_METHOD_CANONICAL_ASSOCIATE, (gr_funcptr) _gr_fmpz_poly_canonical_associate},
924934
{GR_METHOD_FACTOR, (gr_funcptr) _gr_fmpz_poly_factor},
925935
{GR_METHOD_POLY_MULLOW, (gr_funcptr) _gr_fmpz_poly_gr_poly_mullow},
936+
{GR_METHOD_POLY_MULMID, (gr_funcptr) _gr_fmpz_poly_gr_poly_mulmid},
926937
{0, (gr_funcptr) NULL},
927938
};
928939

‎src/gr/fmpzi.c‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -943,6 +943,27 @@ _gr_fmpzi_poly_mullow(fmpzi_struct * res, const fmpzi_struct * poly1, slong len1
943943
}
944944
}
945945

946+
static int
947+
_gr_fmpzi_poly_mulmid(fmpzi_struct * res, const fmpzi_struct * poly1, slong len1, const fmpzi_struct * poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
948+
{
949+
if (len1 < MUL_REORDER_CUTOFF ||
950+
len2 < MUL_REORDER_CUTOFF ||
951+
FLINT_MIN(nhi, len1 + len2 - 1 - nlo) < MUL_REORDER_CUTOFF ||
952+
2 * (nhi - nlo) < MUL_REORDER_CUTOFF)
953+
{
954+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
955+
}
956+
else
957+
{
958+
gr_ctx_t fmpz_ctx;
959+
gr_ctx_init_fmpz(fmpz_ctx);
960+
/* Not currently using Karatsuba because it's often worse if the real
961+
and imaginary parts are even slightly unbalanced */
962+
return _gr_poly_mulmid_complex_reorder(res, poly1, len1, poly2, len2, nlo, nhi, 0, ctx, fmpz_ctx);
963+
}
964+
}
965+
966+
946967
/*
947968
int
948969
_gr_fmpzi_sgn(fmpzi_t res, const fmpzi_t x, const gr_ctx_t ctx)
@@ -1085,6 +1106,7 @@ gr_method_tab_input _fmpzi_methods_input[] =
10851106
{GR_METHOD_VEC_DOT_REV, (gr_funcptr) _gr_fmpzi_vec_dot_rev},
10861107
*/
10871108
{GR_METHOD_POLY_MULLOW, (gr_funcptr) _gr_fmpzi_poly_mullow},
1109+
{GR_METHOD_POLY_MULMID, (gr_funcptr) _gr_fmpzi_poly_mulmid},
10881110
/*
10891111
{GR_METHOD_MAT_MUL, (gr_funcptr) _gr_fmpzi_mat_mul},
10901112
{GR_METHOD_MAT_DET, (gr_funcptr) _gr_fmpzi_mat_det},

‎src/gr/polynomial.c‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -737,6 +737,21 @@ _polynomial_gr_poly_mullow(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr po
737737
return _gr_poly_mullow_bivariate_KS(res, poly1, len1, poly2, len2, n, ctx);
738738
}
739739

740+
static int
741+
_polynomial_gr_poly_mulmid(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
742+
{
743+
gr_ctx_struct * cctx = POLYNOMIAL_ELEM_CTX(ctx);
744+
745+
if (len1 < MUL_KS_CUTOFF || len2 < MUL_KS_CUTOFF || nhi < MUL_KS_CUTOFF
746+
|| len1 + len2 - 1 - nlo < MUL_KS_CUTOFF || 2 * (nhi - nlo) < MUL_KS_CUTOFF)
747+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
748+
749+
if (!want_KS(cctx, 0))
750+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
751+
752+
return _gr_poly_mulmid_bivariate_KS(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
753+
}
754+
740755

741756
int _gr_poly_methods_initialized = 0;
742757

@@ -829,6 +844,7 @@ gr_method_tab_input _gr_poly_methods_input[] =
829844

830845
{GR_METHOD_FACTOR, (gr_funcptr) polynomial_factor},
831846
{GR_METHOD_POLY_MULLOW, (gr_funcptr) _polynomial_gr_poly_mullow},
847+
{GR_METHOD_POLY_MULMID, (gr_funcptr) _polynomial_gr_poly_mulmid},
832848

833849
{0, (gr_funcptr) NULL},
834850
};

‎src/gr_poly.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,8 @@ GR_POLY_INLINE WARN_UNUSED_RESULT int _gr_poly_mullow(gr_ptr res, gr_srcptr poly
142142
WARN_UNUSED_RESULT int _gr_poly_mullow_generic(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx);
143143
WARN_UNUSED_RESULT int gr_poly_mullow(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx);
144144

145+
WARN_UNUSED_RESULT int _gr_poly_mulmid_complex_reorder(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx);
146+
WARN_UNUSED_RESULT int gr_poly_mulmid_complex_reorder(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx);
145147
WARN_UNUSED_RESULT int _gr_poly_mullow_complex_reorder(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx);
146148
WARN_UNUSED_RESULT int gr_poly_mullow_complex_reorder(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, int karatsuba, gr_ctx_t ctx, gr_ctx_t real_ctx);
147149

@@ -162,6 +164,9 @@ WARN_UNUSED_RESULT int gr_poly_mulmid(gr_poly_t res, const gr_poly_t poly1, cons
162164
WARN_UNUSED_RESULT int _gr_poly_mulmid_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
163165
WARN_UNUSED_RESULT int gr_poly_mulmid_classical(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, gr_ctx_t ctx);
164166

167+
WARN_UNUSED_RESULT int _gr_poly_mulmid_bivariate_KS(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
168+
WARN_UNUSED_RESULT int gr_poly_mulmid_bivariate_KS(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, gr_ctx_t ctx);
169+
165170
GR_POLY_INLINE WARN_UNUSED_RESULT int
166171
_gr_poly_mullow_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
167172
{

‎src/gr_poly/mullow_bivariate_KS.c‎

Lines changed: 83 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -15,17 +15,17 @@
1515
#include "gr_poly.h"
1616

1717
int
18-
_gr_poly_mullow_bivariate_KS(gr_ptr res,
18+
_gr_poly_mulmid_bivariate_KS(gr_ptr res,
1919
gr_srcptr poly1, slong len1,
20-
gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
20+
gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
2121
{
22-
slong i, l, max_len1, max_len2, inner_len;
22+
slong i, l, max_len1, max_len2, inner_len, n;
2323
gr_ctx_struct * cctx;
2424
gr_ctx_t tmp_ctx;
2525
const gr_poly_struct * ppoly1, * ppoly2;
2626
gr_poly_struct * pres;
2727
gr_ptr R, P1, P2, pp, tmp_zero;
28-
slong csz;
28+
slong sz, csz;
2929
int status = GR_SUCCESS;
3030
int squaring;
3131

@@ -42,10 +42,36 @@ _gr_poly_mullow_bivariate_KS(gr_ptr res,
4242
else
4343
return GR_UNABLE;
4444

45+
sz = ctx->sizeof_elem;
4546
csz = cctx->sizeof_elem;
4647

47-
len1 = FLINT_MIN(len1, n);
48-
len2 = FLINT_MIN(len2, n);
48+
len1 = FLINT_MIN(len1, nhi);
49+
len2 = FLINT_MIN(len2, nhi);
50+
51+
if (nlo != 0)
52+
{
53+
slong nlo2 = (len1 + len2 - 1) - nlo;
54+
55+
if (len1 > nlo2)
56+
{
57+
slong trunc = len1 - nlo2;
58+
poly1 = GR_ENTRY(poly1, trunc, sz);
59+
len1 -= trunc;
60+
nlo -= trunc;
61+
nhi -= trunc;
62+
}
63+
64+
if (len2 > nlo2)
65+
{
66+
slong trunc = len2 - nlo2;
67+
poly2 = GR_ENTRY(poly2, trunc, sz);
68+
len2 -= trunc;
69+
nlo -= trunc;
70+
nhi -= trunc;
71+
}
72+
}
73+
74+
n = nhi - nlo;
4975

5076
squaring = (poly1 == poly2) && (len1 == len2);
5177
ppoly1 = poly1;
@@ -81,12 +107,12 @@ _gr_poly_mullow_bivariate_KS(gr_ptr res,
81107

82108
GR_TMP_INIT_VEC(R, n * inner_len, cctx);
83109
GR_TMP_INIT_VEC(tmp_zero, inner_len, cctx);
84-
P1 = GR_TMP_ALLOC(len1 * inner_len * cctx->sizeof_elem);
110+
P1 = GR_TMP_ALLOC(len1 * inner_len * csz);
85111

86112
if (squaring)
87113
P2 = P1;
88114
else
89-
P2 = GR_TMP_ALLOC(len2 * inner_len * cctx->sizeof_elem);
115+
P2 = GR_TMP_ALLOC(len2 * inner_len * csz);
90116

91117
for (i = 0; i < len1; i++)
92118
{
@@ -105,7 +131,7 @@ _gr_poly_mullow_bivariate_KS(gr_ptr res,
105131
}
106132
}
107133

108-
status = _gr_poly_mullow(R, P1, len1 * inner_len, P2, len2 * inner_len, n * inner_len, cctx);
134+
status = _gr_poly_mulmid(R, P1, len1 * inner_len, P2, len2 * inner_len, nlo * inner_len, nhi * inner_len, cctx);
109135

110136
for (i = 0; i < n; i++)
111137
{
@@ -120,10 +146,56 @@ _gr_poly_mullow_bivariate_KS(gr_ptr res,
120146

121147
GR_TMP_CLEAR_VEC(R, n * inner_len, cctx);
122148
GR_TMP_CLEAR_VEC(tmp_zero, inner_len, cctx);
123-
GR_TMP_FREE(P1, len1 * inner_len * cctx->sizeof_elem);
149+
GR_TMP_FREE(P1, len1 * inner_len * csz);
124150
if (!squaring)
125-
GR_TMP_FREE(P2, len2 * inner_len * cctx->sizeof_elem);
151+
GR_TMP_FREE(P2, len2 * inner_len * csz);
152+
153+
return status;
154+
}
155+
156+
int
157+
_gr_poly_mullow_bivariate_KS(gr_ptr res,
158+
gr_srcptr poly1, slong len1,
159+
gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
160+
{
161+
return _gr_poly_mulmid_bivariate_KS(res, poly1, len1, poly2, len2, 0, n, ctx);
162+
}
163+
164+
int
165+
gr_poly_mulmid_bivariate_KS(gr_poly_t res, const gr_poly_t poly1,
166+
const gr_poly_t poly2,
167+
slong nlo, slong nhi, gr_ctx_t ctx)
168+
{
169+
slong len1 = poly1->length;
170+
slong len2 = poly2->length;
171+
int status;
172+
slong len;
126173

174+
FLINT_ASSERT(nlo >= 0);
175+
FLINT_ASSERT(nhi >= 0);
176+
177+
if (len1 == 0 || len2 == 0 || nlo >= FLINT_MIN(nhi, len1 + len2 - 1))
178+
return gr_poly_zero(res, ctx);
179+
180+
nhi = FLINT_MIN(nhi, len1 + len2 - 1);
181+
len = nhi - nlo;
182+
183+
if (res == poly1 || res == poly2)
184+
{
185+
gr_poly_t t;
186+
gr_poly_init2(t, len, ctx);
187+
status = _gr_poly_mulmid_bivariate_KS(t->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
188+
gr_poly_swap(res, t, ctx);
189+
gr_poly_clear(t, ctx);
190+
}
191+
else
192+
{
193+
gr_poly_fit_length(res, len, ctx);
194+
status = _gr_poly_mulmid_bivariate_KS(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
195+
}
196+
197+
_gr_poly_set_length(res, len, ctx);
198+
_gr_poly_normalise(res, ctx);
127199
return status;
128200
}
129201

0 commit comments

Comments
 (0)