|
| 1 | +/* |
| 2 | + Copyright (C) 2025 Marc Mezzarobba |
| 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 | +/* todo: make it possible to compute a superset of the result, implement the |
| 13 | + * algorithm of Gerhard, Giesbrecht, Storjohann, Zima (2003) */ |
| 14 | + |
| 15 | +#include "fmpz.h" |
| 16 | +#include "fmpz_vec.h" |
| 17 | +#include "gr.h" |
| 18 | +#include "gr_poly.h" |
| 19 | +#include "gr_vec.h" |
| 20 | + |
| 21 | +int |
| 22 | +gr_poly_dispersion_resultant(fmpz_t disp, gr_vec_t disp_set, |
| 23 | + const gr_poly_t f, const gr_poly_t g, |
| 24 | + gr_ctx_t ctx) |
| 25 | +{ |
| 26 | + gr_ctx_t Pol_a, ZZ; |
| 27 | + gr_poly_t a, fa, ga, res; |
| 28 | + gr_vec_t roots, mult; |
| 29 | + |
| 30 | + truth_t is_zero = truth_and(gr_poly_is_zero(f, ctx), |
| 31 | + gr_poly_is_zero(g, ctx)); |
| 32 | + if (is_zero == T_TRUE) |
| 33 | + return GR_DOMAIN; |
| 34 | + else if (is_zero == T_UNKNOWN) |
| 35 | + return GR_UNABLE; |
| 36 | + |
| 37 | + gr_ctx_init_gr_poly(Pol_a, ctx); |
| 38 | + gr_ctx_init_fmpz(ZZ); |
| 39 | + |
| 40 | + gr_poly_init(a, ctx); |
| 41 | + gr_poly_init(fa, Pol_a); |
| 42 | + gr_poly_init(ga, Pol_a); |
| 43 | + gr_poly_init(res, ZZ); |
| 44 | + |
| 45 | + gr_vec_init(roots, 0, ZZ); |
| 46 | + gr_vec_init(mult, 0, ZZ); |
| 47 | + |
| 48 | + int status = GR_SUCCESS; |
| 49 | + |
| 50 | + status |= gr_poly_set_gr_poly_other(fa, f, ctx, Pol_a); |
| 51 | + status |= gr_poly_set_gr_poly_other(ga, g, ctx, Pol_a); |
| 52 | + |
| 53 | + status |= gr_poly_gen(a, ctx); |
| 54 | + status |= gr_poly_taylor_shift(fa, fa, a, Pol_a); |
| 55 | + |
| 56 | + /* todo: composed sum */ |
| 57 | + status |= gr_poly_resultant(res, fa, ga, Pol_a); |
| 58 | + |
| 59 | + status |= gr_poly_roots_other(roots, mult, res, ctx, 0, ZZ); |
| 60 | + |
| 61 | + if (disp_set != NULL) |
| 62 | + gr_vec_set_length(disp_set, 0, ZZ); |
| 63 | + |
| 64 | + for (slong i = 0; i < roots->length; i++) |
| 65 | + { |
| 66 | + gr_srcptr rt = gr_vec_entry_srcptr(roots, i, ZZ); |
| 67 | + if (fmpz_sgn(rt) >= 0) { |
| 68 | + if (disp_set != NULL) |
| 69 | + status |= gr_vec_append(disp_set, rt, ZZ); |
| 70 | + if (disp != NULL && fmpz_cmp(rt, disp) >= 0) |
| 71 | + fmpz_swap(disp, (fmpz *) rt); |
| 72 | + } |
| 73 | + } |
| 74 | + |
| 75 | + if (disp_set != NULL) |
| 76 | + _fmpz_vec_sort(disp_set->entries, disp_set->length); |
| 77 | + |
| 78 | + gr_vec_clear(mult, ZZ); |
| 79 | + gr_vec_clear(roots, ZZ); |
| 80 | + |
| 81 | + gr_poly_clear(res, ctx); |
| 82 | + gr_poly_clear(ga, Pol_a); |
| 83 | + gr_poly_clear(fa, Pol_a); |
| 84 | + gr_poly_clear(a, ctx); |
| 85 | + |
| 86 | + gr_ctx_clear(ZZ); |
| 87 | + gr_ctx_clear(Pol_a); |
| 88 | + |
| 89 | + return status; |
| 90 | +} |
| 91 | + |
| 92 | +int |
| 93 | +gr_poly_dispersion_from_factors(fmpz_t disp, gr_vec_t disp_set, |
| 94 | + const gr_vec_t ffac, const gr_vec_t gfac, |
| 95 | + gr_ctx_t ctx) |
| 96 | +{ |
| 97 | + gr_ctx_t ZZ, pctx; |
| 98 | + gr_ptr a, b, _alpha; |
| 99 | + fmpz_t alpha; |
| 100 | + gr_poly_t pshift; |
| 101 | + |
| 102 | + if (gr_ctx_is_unique_factorization_domain(ctx) != T_TRUE) |
| 103 | + return GR_UNABLE; |
| 104 | + |
| 105 | + /* todo: generalize to positive characteristic? (cf. Gerhard et al.) */ |
| 106 | + if (gr_ctx_is_finite_characteristic(ctx) != T_FALSE) |
| 107 | + return GR_UNABLE; |
| 108 | + |
| 109 | + gr_ctx_init_fmpz(ZZ); |
| 110 | + |
| 111 | + if (disp_set != NULL) |
| 112 | + gr_vec_set_length(disp_set, 0, ZZ); |
| 113 | + |
| 114 | + /* We assume that f and g are nonzero, so no factors means constant */ |
| 115 | + if (ffac->length == 0 || gfac->length == 0) |
| 116 | + return GR_SUCCESS; |
| 117 | + |
| 118 | + gr_ctx_init_gr_poly(pctx, ctx); |
| 119 | + GR_TMP_INIT3(a, b, _alpha, ctx); |
| 120 | + fmpz_init(alpha); |
| 121 | + gr_poly_init(pshift, ctx); |
| 122 | + |
| 123 | + int status = GR_SUCCESS; |
| 124 | + |
| 125 | + if (ffac == gfac) |
| 126 | + { |
| 127 | + if (disp_set != NULL) |
| 128 | + status |= gr_vec_append(disp_set, alpha, ZZ); |
| 129 | + if (disp != NULL) |
| 130 | + fmpz_zero(disp); |
| 131 | + } |
| 132 | + |
| 133 | + for (slong i = 0; i < ffac->length; i++) |
| 134 | + { |
| 135 | + const gr_poly_struct *p = gr_vec_entry_srcptr(ffac, i, pctx); |
| 136 | + |
| 137 | + slong j0 = (ffac == gfac) ? i + 1 : 0; |
| 138 | + for (slong j = j0; j < gfac->length; j++) |
| 139 | + { |
| 140 | + const gr_poly_struct *q = gr_vec_entry_srcptr(gfac, j, pctx); |
| 141 | + |
| 142 | + int status1 = GR_SUCCESS; |
| 143 | + |
| 144 | + status1 |= gr_poly_leading_taylor_shift(_alpha, p, q, ctx); |
| 145 | + if (status1 == GR_DOMAIN) |
| 146 | + continue; |
| 147 | + status1 |= gr_get_fmpz(alpha, _alpha, ctx); |
| 148 | + if (status1 == GR_DOMAIN) |
| 149 | + continue; |
| 150 | + status |= status1; |
| 151 | + |
| 152 | + if (status != GR_SUCCESS) |
| 153 | + goto cleanup; |
| 154 | + |
| 155 | + if (ffac == gfac) |
| 156 | + fmpz_abs(alpha, alpha); |
| 157 | + else if (fmpz_sgn(alpha) < 0) |
| 158 | + continue; |
| 159 | + |
| 160 | + if (disp_set != NULL |
| 161 | + ? gr_vec_contains(disp_set, alpha, ZZ) == T_TRUE |
| 162 | + : disp != NULL && fmpz_cmp(alpha, disp) <= 0) |
| 163 | + continue; |
| 164 | + |
| 165 | + if (p->length > 2) /* GR_SUCCESS => degree well-defined */ |
| 166 | + { |
| 167 | + status |= gr_poly_taylor_shift(pshift, p, _alpha, ctx); |
| 168 | + int eq = gr_poly_equal(pshift, q, ctx); |
| 169 | + if (eq == T_FALSE) |
| 170 | + continue; |
| 171 | + else if (eq == T_UNKNOWN) |
| 172 | + status = GR_UNABLE; |
| 173 | + } |
| 174 | + |
| 175 | + if (status != GR_SUCCESS) |
| 176 | + goto cleanup; |
| 177 | + |
| 178 | + if (disp_set != NULL) |
| 179 | + status |= gr_vec_append(disp_set, alpha, ZZ); |
| 180 | + if (disp != NULL && fmpz_cmp(alpha, disp) > 0) |
| 181 | + fmpz_swap(disp, alpha); |
| 182 | + } |
| 183 | + } |
| 184 | + |
| 185 | + if (disp_set != NULL) |
| 186 | + _fmpz_vec_sort(disp_set->entries, disp_set->length); |
| 187 | + |
| 188 | +cleanup: |
| 189 | + |
| 190 | + gr_poly_clear(pshift, ctx); |
| 191 | + fmpz_clear(alpha); |
| 192 | + GR_TMP_CLEAR3(a, b, _alpha, ctx); |
| 193 | + gr_ctx_clear(pctx); |
| 194 | + gr_ctx_clear(ZZ); |
| 195 | + |
| 196 | + return status; |
| 197 | +} |
| 198 | + |
| 199 | +int |
| 200 | +gr_poly_dispersion_factor(fmpz_t disp, gr_vec_t disp_set, |
| 201 | + const gr_poly_t f, const gr_poly_t g, |
| 202 | + gr_ctx_t ctx) |
| 203 | +{ |
| 204 | + gr_ctx_t pctx, ZZ; |
| 205 | + gr_poly_t fsqf, gsqf; |
| 206 | + gr_vec_t ffac, gfac, ignored; |
| 207 | + gr_ptr c; |
| 208 | + |
| 209 | + switch (truth_and(gr_poly_is_zero(f, ctx), gr_poly_is_zero(g, ctx))) |
| 210 | + { |
| 211 | + case T_TRUE: |
| 212 | + return GR_DOMAIN; |
| 213 | + case T_UNKNOWN: |
| 214 | + return GR_UNABLE; |
| 215 | + case T_FALSE: |
| 216 | + ; |
| 217 | + } |
| 218 | + |
| 219 | + int status = GR_SUCCESS; |
| 220 | + |
| 221 | + gr_ctx_init_fmpz(ZZ); |
| 222 | + gr_ctx_init_gr_poly(pctx, ctx); |
| 223 | + gr_poly_init(fsqf, ctx); |
| 224 | + gr_poly_init(gsqf, ctx); |
| 225 | + gr_vec_init(ffac, 0, pctx); |
| 226 | + gr_vec_init(gfac, 0, pctx); |
| 227 | + gr_vec_init(ignored, 0, ZZ); |
| 228 | + GR_TMP_INIT(c, pctx); |
| 229 | + |
| 230 | + if (gr_poly_squarefree_part(fsqf, f, ctx) != GR_SUCCESS) |
| 231 | + status |= gr_poly_set(fsqf, f, ctx); |
| 232 | + status |= gr_factor(c, ffac, ignored, f, 0, pctx); |
| 233 | + |
| 234 | + if (f == g) |
| 235 | + { |
| 236 | + status |= gr_poly_dispersion_from_factors(disp, disp_set, ffac, ffac, ctx); |
| 237 | + } |
| 238 | + else |
| 239 | + { |
| 240 | + if (gr_poly_squarefree_part(gsqf, g, ctx) != GR_SUCCESS) |
| 241 | + status |= gr_poly_set(gsqf, g, ctx); |
| 242 | + status |= gr_factor(c, gfac, ignored, g, 0, pctx); |
| 243 | + status |= gr_poly_dispersion_from_factors(disp, disp_set, ffac, gfac, ctx); |
| 244 | + } |
| 245 | + |
| 246 | + GR_TMP_CLEAR(c, pctx); |
| 247 | + gr_vec_clear(gfac, pctx); |
| 248 | + gr_vec_clear(ffac, pctx); |
| 249 | + gr_poly_clear(fsqf, ctx); |
| 250 | + gr_poly_clear(gsqf, ctx); |
| 251 | + gr_vec_clear(ignored, ZZ); |
| 252 | + gr_ctx_clear(pctx); |
| 253 | + gr_ctx_clear(ZZ); |
| 254 | + |
| 255 | + return status; |
| 256 | +} |
| 257 | + |
| 258 | +int |
| 259 | +gr_poly_dispersion(fmpz_t disp, gr_vec_t disp_set, |
| 260 | + const gr_poly_t f, const gr_poly_t g, |
| 261 | + gr_ctx_t ctx) |
| 262 | +{ |
| 263 | + return gr_poly_dispersion_factor(disp, disp_set, f, g, ctx); |
| 264 | +} |
0 commit comments