1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
/*
Copyright (C) 2020 Daniel Schultz
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "mpoly.h"
#include "long_extras.h"
/*
for f in F[y][x] of degree n - 1 in x
on input: l[i] = 1 + deg_y(coeff(f, x^i)) for 0 <= i <= n
on output: l[i] > deg_y(coeff(f*g_x/g, x^i)) for any divisor g of f
algo is a 2x scan to find the upper hull and then to fill in the output
*/
void mpoly_bivar_cld_bounds(slong * l, slong n)
{
slong * P;
slong Plen = 0;
slong i, j, x0, y0, x1, y1;
TMP_INIT;
FLINT_ASSERT(n > 0);
TMP_START;
P = (slong *) TMP_ALLOC(2*n*sizeof(slong));
n--;
P[0] = n; P[1] = l[n];
Plen = 1;
for (i = n - 1; i >= 0; i--)
{
if (l[i] < 1)
continue;
x0 = i; y0 = l[i];
while (Plen >= 2 && !z_mat22_det_is_negative(
P[2*(Plen - 1) + 0] - x0, P[2*(Plen - 1) + 1] - y0,
P[2*(Plen - 2) + 0] - x0, P[2*(Plen - 2) + 1] - y0))
Plen--;
P[2*Plen + 0] = x0;
P[2*Plen + 1] = y0;
Plen++;
}
i = Plen - 1;
x0 = P[2*i + 0]; y0 = P[2*i + 1];
for (j = 1; j <= x0; j++)
l[j - 1] = j < x0 ? 0 : y0;
FLINT_ASSERT(j == x0 + 1);
while (i > 0)
{
x1 = P[2*(i - 1) + 0]; y1 = P[2*(i - 1) + 1];
for ( ; j <= x1; j++)
{
ulong t1, t2, t3, t4;
FLINT_ASSERT(x0 < j && j <= x1);
umul_ppmm(t1, t2, j - x0, y1);
umul_ppmm(t3, t4, x1 - j, y0);
add_ssaaaa(t1, t2, t1, t2, t3, t4);
udiv_qrnnd(l[j - 1], t3, t1, t2, x1 - x0);
}
x0 = x1; y0 = y1;
i--;
}
FLINT_ASSERT(j == n + 1);
l[j - 1] = 0;
TMP_END;
}