Skip to content

Commit 17e900d

Browse files
fredrik-johanssonfredrik-johansson
committed
Add _gr_poly_mulmid_classical and use in mullow_classical, mulmid_generic
1 parent 56accbc commit 17e900d

File tree

6 files changed

+237
-45
lines changed

6 files changed

+237
-45
lines changed

‎doc/source/gr_poly.rst‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -183,8 +183,10 @@ 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(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_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)
187188
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)
189+
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)
188190
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)
189191

190192
Multiplication algorithms
@@ -1309,11 +1311,13 @@ polynomials up to length *maxn*. If *ctx* is set to ``NULL``, a random
13091311
ring is generated on each test iteration, otherwise the given ring is used.
13101312

13111313
.. function:: void _gr_poly_test_mullow(gr_method_poly_binary_trunc_op mullow_impl, gr_method_poly_binary_trunc_op mullow_ref, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
1314+
void _gr_poly_test_mulmid(gr_method_poly_binary_trunc_op mulmid_impl, gr_method_poly_binary_trunc_op mulmid_ref, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
13121315

1313-
Tests the given function ``mullow_impl`` for correctness as an implementation
1314-
of :func:`_gr_poly_mullow`.
1316+
Tests the given function ``mullow_impl`` (``mulmid_impl``) for correctness as an implementation
1317+
of :func:`_gr_poly_mullow` (:func:`_gr_poly_mulmid`).
13151318
A reference implementation to compare against can be provided as
1316-
``mullow_ref``; if ``NULL``, classical multiplication is used.
1319+
``mullow_ref`` (``mulmid_ref``); if ``NULL``, classical
1320+
multiplication is used.
13171321

13181322
.. function:: void _gr_poly_test_divrem(gr_method_poly_binary_binary_op divrem_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
13191323

‎src/gr_poly.h‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -142,9 +142,6 @@ 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_mullow_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx);
146-
WARN_UNUSED_RESULT 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);
147-
148145
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);
149146
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);
150147

@@ -162,6 +159,17 @@ GR_POLY_INLINE WARN_UNUSED_RESULT int _gr_poly_mulmid(gr_ptr res, gr_srcptr poly
162159
WARN_UNUSED_RESULT 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);
163160
WARN_UNUSED_RESULT 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);
164161

162+
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);
163+
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);
164+
165+
GR_POLY_INLINE WARN_UNUSED_RESULT int
166+
_gr_poly_mullow_classical(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
167+
{
168+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, 0, n, ctx);
169+
}
170+
171+
WARN_UNUSED_RESULT 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);
172+
165173
WARN_UNUSED_RESULT int _gr_poly_mul_karatsuba(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx);
166174
WARN_UNUSED_RESULT int gr_poly_mul_karatsuba(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx);
167175
WARN_UNUSED_RESULT int _gr_poly_mul_toom33(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx);

‎src/gr_poly/mullow_classical.c‎

Lines changed: 93 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -13,40 +13,70 @@
1313
#include "gr_poly.h"
1414

1515
int
16-
_gr_poly_mullow_classical(gr_ptr res,
16+
_gr_poly_mulmid_classical(gr_ptr res,
1717
gr_srcptr poly1, slong len1,
18-
gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
18+
gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
1919
{
20-
int status;
21-
len1 = FLINT_MIN(len1, n);
22-
len2 = FLINT_MIN(len2, n);
20+
int status = GR_SUCCESS;
21+
slong sz = ctx->sizeof_elem;
22+
23+
FLINT_ASSERT(len1 != 0);
24+
FLINT_ASSERT(len2 != 0);
25+
FLINT_ASSERT(nhi != 0);
26+
FLINT_ASSERT(nlo < nhi);
27+
FLINT_ASSERT(nlo >= 0);
28+
FLINT_ASSERT(nhi <= len1 + len2 - 1);
29+
30+
len1 = FLINT_MIN(len1, nhi);
31+
len2 = FLINT_MIN(len2, nhi);
2332

24-
if (n == 1)
33+
if (nlo != 0)
34+
{
35+
slong nlo2 = (len1 + len2 - 1) - nlo;
36+
37+
if (len1 > nlo2)
38+
{
39+
slong trunc = len1 - nlo2;
40+
poly1 = GR_ENTRY(poly1, trunc, sz);
41+
len1 -= trunc;
42+
nlo -= trunc;
43+
nhi -= trunc;
44+
}
45+
46+
if (len2 > nlo2)
47+
{
48+
slong trunc = len2 - nlo2;
49+
poly2 = GR_ENTRY(poly2, trunc, sz);
50+
len2 -= trunc;
51+
nlo -= trunc;
52+
nhi -= trunc;
53+
}
54+
}
55+
56+
if (nhi == 1)
2557
return gr_mul(res, poly1, poly2, ctx);
2658

2759
if (len1 == 1)
28-
return _gr_scalar_mul_vec(res, poly1, poly2, n, ctx);
60+
return _gr_scalar_mul_vec(res, poly1, GR_ENTRY(poly2, nlo, sz), nhi - nlo, ctx);
2961

3062
if (len2 == 1)
31-
return _gr_vec_mul_scalar(res, poly1, n, poly2, ctx);
63+
return _gr_vec_mul_scalar(res, GR_ENTRY(poly1, nlo, sz), nhi - nlo, poly2, ctx);
64+
65+
res = GR_ENTRY(res, -nlo, sz);
3266

3367
/* Squaring */
3468
/* fixme: only valid for commutative rings, but we also want squaring over
3569
ringlike structures, e.g. floats */
3670
if (poly1 == poly2 && len1 == len2 /* && gr_ctx_is_commutative_ring(ctx) == T_TRUE */)
3771
{
38-
slong i, start, stop, sz;
39-
sz = ctx->sizeof_elem;
40-
41-
status = GR_SUCCESS;
72+
slong i, start, stop;
4273

4374
/* todo: double, square, addmul */
4475

45-
status |= gr_sqr(res, poly1, ctx);
46-
status |= gr_mul(GR_ENTRY(res, 1, sz), poly1, GR_ENTRY(poly1, 1, sz), ctx);
47-
status |= gr_mul_two(GR_ENTRY(res, 1, sz), GR_ENTRY(res, 1, sz), ctx);
76+
if (nlo == 0)
77+
status |= gr_sqr(res, poly1, ctx);
4878

49-
for (i = 2; i < FLINT_MIN(n, 2 * len1 - 3); i++)
79+
for (i = FLINT_MAX(1, nlo); i < FLINT_MIN(nhi, 2 * len1 - 2); i++)
5080
{
5181
start = FLINT_MAX(0, i - len1 + 1);
5282
stop = FLINT_MIN(len1 - 1, (i + 1) / 2 - 1);
@@ -58,36 +88,72 @@ _gr_poly_mullow_classical(gr_ptr res,
5888
status |= gr_addmul(GR_ENTRY(res, i, sz), GR_ENTRY(poly1, i / 2, sz), GR_ENTRY(poly1, i / 2, sz), ctx);
5989
}
6090

61-
if (len1 > 2 && n >= 2 * len1 - 2)
62-
{
63-
status |= gr_mul(GR_ENTRY(res, 2 * len1 - 3, sz), GR_ENTRY(poly1, len1 - 1, sz), GR_ENTRY(poly1, len1 - 2, sz), ctx);
64-
status |= gr_mul_two(GR_ENTRY(res, 2 * len1 - 3, sz), GR_ENTRY(res, 2 * len1 - 3, sz), ctx);
65-
}
66-
67-
if (n >= 2 * len1 - 1)
91+
if (nhi >= 2 * len1 - 1)
6892
status |= gr_sqr(GR_ENTRY(res, 2 * len1 - 2, sz), GR_ENTRY(poly1, len1 - 1, sz), ctx);
6993

7094
return status;
7195
}
7296
else
7397
{
74-
slong i, top1, top2, sz;
75-
sz = ctx->sizeof_elem;
98+
slong i, top1, top2;
7699

77-
status = gr_mul(res, poly1, poly2, ctx);
100+
if (nlo == 0)
101+
status = gr_mul(res, poly1, poly2, ctx);
78102

79-
for (i = 1; i < n; i++)
103+
for (i = FLINT_MAX(1, nlo); i < FLINT_MIN(nhi, len1 + len2 - 2); i++)
80104
{
81105
top1 = FLINT_MIN(len1 - 1, i);
82106
top2 = FLINT_MIN(len2 - 1, i);
83107

84108
status |= _gr_vec_dot_rev(GR_ENTRY(res, i, sz), NULL, 0, GR_ENTRY(poly1, i - top2, sz), GR_ENTRY(poly2, i - top1, sz), top1 + top2 - i + 1, ctx);
85109
}
86110

111+
if (nhi >= len1 + len2 - 1)
112+
status |= gr_mul(GR_ENTRY(res, len1 + len2 - 2, sz), GR_ENTRY(poly1, len1 - 1, sz), GR_ENTRY(poly2, len2 - 1, sz), ctx);
113+
87114
return status;
88115
}
89116
}
90117

118+
int
119+
gr_poly_mulmid_classical(gr_poly_t res, const gr_poly_t poly1,
120+
const gr_poly_t poly2,
121+
slong nlo, slong nhi, gr_ctx_t ctx)
122+
{
123+
slong len1 = poly1->length;
124+
slong len2 = poly2->length;
125+
int status;
126+
slong len;
127+
128+
FLINT_ASSERT(nlo >= 0);
129+
FLINT_ASSERT(nhi >= 0);
130+
131+
if (len1 == 0 || len2 == 0 || nlo >= FLINT_MIN(nhi, len1 + len2 - 1))
132+
return gr_poly_zero(res, ctx);
133+
134+
nhi = FLINT_MIN(nhi, len1 + len2 - 1);
135+
len = nhi - nlo;
136+
137+
if (res == poly1 || res == poly2)
138+
{
139+
gr_poly_t t;
140+
gr_poly_init2(t, len, ctx);
141+
status = _gr_poly_mulmid_classical(t->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
142+
gr_poly_swap(res, t, ctx);
143+
gr_poly_clear(t, ctx);
144+
}
145+
else
146+
{
147+
gr_poly_fit_length(res, len, ctx);
148+
status = _gr_poly_mulmid_classical(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
149+
}
150+
151+
_gr_poly_set_length(res, len, ctx);
152+
_gr_poly_normalise(res, ctx);
153+
return status;
154+
}
155+
156+
91157
int
92158
gr_poly_mullow_classical(gr_poly_t res, const gr_poly_t poly1,
93159
const gr_poly_t poly2,

‎src/gr_poly/mulmid.c‎

Lines changed: 18 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -12,30 +12,37 @@
1212
#include "gr_vec.h"
1313
#include "gr_poly.h"
1414

15-
/* Todo: a reasonable generic implementation. Currently many rings
16-
implement a good mullow, so fall back on that. */
1715
int
1816
_gr_poly_mulmid_generic(gr_ptr res,
1917
gr_srcptr poly1, slong len1,
2018
gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
2119
{
22-
gr_ptr tmp;
23-
slong sz = ctx->sizeof_elem;
24-
int status;
25-
2620
FLINT_ASSERT(len1 != 0);
2721
FLINT_ASSERT(len2 != 0);
2822
FLINT_ASSERT(nhi != 0);
2923
FLINT_ASSERT(nlo < nhi);
3024
FLINT_ASSERT(nlo >= 0);
3125
FLINT_ASSERT(nhi <= len1 + len2 - 1);
3226

33-
GR_TMP_INIT_VEC(tmp, nhi, ctx);
34-
status = _gr_poly_mullow(tmp, poly1, len1, poly2, len2, nhi, ctx);
35-
_gr_vec_swap(res, GR_ENTRY(tmp, nlo, sz), nhi - nlo, ctx);
36-
GR_TMP_CLEAR_VEC(tmp, nhi, ctx);
27+
/* Todo: a reasonable generic implementation. Currently many rings
28+
implement a good mullow, so fall back on that. */
29+
if (ctx->methods[GR_METHOD_POLY_MULLOW] != (gr_funcptr) _gr_poly_mullow_generic)
30+
{
31+
gr_ptr tmp;
32+
slong sz = ctx->sizeof_elem;
33+
int status;
3734

38-
return status;
35+
GR_TMP_INIT_VEC(tmp, nhi, ctx);
36+
status = _gr_poly_mullow(tmp, poly1, len1, poly2, len2, nhi, ctx);
37+
_gr_vec_swap(res, GR_ENTRY(tmp, nlo, sz), nhi - nlo, ctx);
38+
GR_TMP_CLEAR_VEC(tmp, nhi, ctx);
39+
40+
return status;
41+
}
42+
else
43+
{
44+
return _gr_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, ctx);
45+
}
3946
}
4047

4148
int

‎src/gr_poly/test/main.c‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,7 @@
6060
#include "t-mullow_complex_reorder.c"
6161
#include "t-mullow_toom_serial.c"
6262
#include "t-mulmid.c"
63+
#include "t-mulmid_classical.c"
6364
#include "t-mulmod.c"
6465
#include "t-mulmod_preinv.c"
6566
#include "t-newton_basis.c"
@@ -145,6 +146,7 @@ test_struct tests[] =
145146
TEST_FUNCTION(gr_poly_mullow_complex_reorder),
146147
TEST_FUNCTION(gr_poly_mullow_toom_serial),
147148
TEST_FUNCTION(gr_poly_mulmid),
149+
TEST_FUNCTION(gr_poly_mulmid_classical),
148150
TEST_FUNCTION(gr_poly_mulmod),
149151
TEST_FUNCTION(gr_poly_mulmod_preinv),
150152
TEST_FUNCTION(gr_poly_newton_basis),

‎src/gr_poly/test/t-mulmid_classical.c‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
/*
2+
Copyright (C) 2026 Fredrik Johansson
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 3 of the License, or
9+
(at your option) any later version. See <https://www.gnu.org/licenses/>.
10+
*/
11+
12+
#include "test_helpers.h"
13+
#include "ulong_extras.h"
14+
#include "gr_poly.h"
15+
16+
TEST_FUNCTION_START(gr_poly_mulmid_classical, state)
17+
{
18+
slong iter;
19+
20+
for (iter = 0; iter < 1000; iter++)
21+
{
22+
int status;
23+
gr_ctx_t ctx;
24+
gr_poly_t a, b, res1, res2;
25+
slong nlo, nhi, N;
26+
int aliasing;
27+
28+
gr_ctx_init_random(ctx, state);
29+
if (gr_ctx_is_finite(ctx) == T_TRUE || gr_ctx_has_real_prec(ctx) == T_TRUE)
30+
N = 10;
31+
else
32+
N = 4;
33+
34+
gr_poly_init(a, ctx);
35+
gr_poly_init(b, ctx);
36+
gr_poly_init(res1, ctx);
37+
gr_poly_init(res2, ctx);
38+
39+
status = GR_SUCCESS;
40+
41+
nlo = n_randint(state, N);
42+
nhi = n_randint(state, N);
43+
44+
status |= gr_poly_randtest(a, state, 1 + n_randint(state, N), ctx);
45+
status |= gr_poly_randtest(b, state, 1 + n_randint(state, N), ctx);
46+
status |= gr_poly_randtest(res1, state, 1 + n_randint(state, N), ctx);
47+
status |= gr_poly_randtest(res2, state, 1 + n_randint(state, N), ctx);
48+
aliasing = n_randint(state, 5);
49+
50+
if (aliasing == 0)
51+
{
52+
status |= gr_poly_mulmid_classical(res1, a, b, nlo, nhi, ctx);
53+
}
54+
else if (aliasing == 1)
55+
{
56+
status |= gr_poly_set(res1, a, ctx);
57+
status |= gr_poly_mulmid_classical(res1, res1, b, nlo, nhi, ctx);
58+
}
59+
else if (aliasing == 2)
60+
{
61+
status |= gr_poly_set(res1, b, ctx);
62+
status |= gr_poly_mulmid_classical(res1, a, res1, nlo, nhi, ctx);
63+
}
64+
else if (aliasing == 3)
65+
{
66+
status |= gr_poly_mulmid_classical(res1, a, a, nlo, nhi, ctx);
67+
}
68+
else if (aliasing == 4)
69+
{
70+
status |= gr_poly_set(res1, a, ctx);
71+
status |= gr_poly_mulmid_classical(res1, res1, res1, nlo, nhi, ctx);
72+
}
73+
74+
if (status == GR_SUCCESS)
75+
{
76+
if (aliasing == 3 || aliasing == 4)
77+
status |= gr_poly_mullow(res2, a, a, nhi, ctx);
78+
else
79+
status |= gr_poly_mullow(res2, a, b, nhi, ctx);
80+
81+
status |= gr_poly_shift_right(res2, res2, nlo, ctx);
82+
83+
if (status == GR_SUCCESS && gr_poly_equal(res1, res2, ctx) == T_FALSE)
84+
{
85+
flint_printf("FAIL: gr_poly_mulmid\n\n");
86+
gr_ctx_println(ctx);
87+
flint_printf("nlo = %wd, nhi = %wd, aliasing = %d\n", nlo, nhi, aliasing);
88+
flint_printf("a = "); gr_poly_print(a, ctx); flint_printf("\n");
89+
flint_printf("b = "); gr_poly_print(b, ctx); flint_printf("\n");
90+
flint_printf("res1 = "); gr_poly_print(res1, ctx); flint_printf("\n");
91+
flint_printf("res2 = "); gr_poly_print(res2, ctx); flint_printf("\n");
92+
flint_abort();
93+
}
94+
}
95+
96+
gr_poly_clear(a, ctx);
97+
gr_poly_clear(b, ctx);
98+
gr_poly_clear(res1, ctx);
99+
gr_poly_clear(res2, ctx);
100+
101+
gr_ctx_clear(ctx);
102+
}
103+
104+
TEST_FUNCTION_END(state);
105+
}

0 commit comments

Comments
 (0)