This codes assumes that there are at least two variables in the current base ring. No assumption is made regarding the monomial ordering. */
| 1217 | /* This codes assumes that there are at least two variables in the current |
| 1218 | base ring. No assumption is made regarding the monomial ordering. */ |
| 1219 | void henselFactors(const int xIndex, const int yIndex, const poly h, |
| 1220 | const poly f0, const poly g0, const int d, poly &f, poly &g) |
| 1221 | { |
| 1222 | int n = (int)p_Deg(f0,currRing); |
| 1223 | int m = (int)p_Deg(g0,currRing); |
| 1224 | matrix aMat = mpNew(n + m, n + m); /* matrix A for linear system */ |
| 1225 | matrix pMat; matrix lMat; matrix uMat; /* for the decomposition of A */ |
| 1226 | f = pCopy(f0); g = pCopy(g0); /* initially: h = f*g mod <x^1> */ |
| 1227 | |
| 1228 | /* initial step: read off coefficients of f0, and g0 */ |
| 1229 | poly p = f0; poly matEntry; number c; |
| 1230 | while (p != NULL) |
| 1231 | { |
| 1232 | c = nCopy(pGetCoeff(p)); |
| 1233 | matEntry = pOne(); pSetCoeff(matEntry, c); |
| 1234 | MATELEM(aMat, pGetExp(p, yIndex) + 1, 1) = matEntry; |
| 1235 | p = pNext(p); |
| 1236 | } |
| 1237 | p = g0; |
| 1238 | while (p != NULL) |
| 1239 | { |
| 1240 | c = nCopy(pGetCoeff(p)); |
| 1241 | matEntry = pOne(); pSetCoeff(matEntry, c); |
| 1242 | MATELEM(aMat, pGetExp(p, yIndex) + 1, m + 1) = matEntry; |
| 1243 | p = pNext(p); |
| 1244 | } |
| 1245 | /* fill the rest of A */ |
| 1246 | for (int row = 2; row <= n + 1; row++) |
| 1247 | for (int col = 2; col <= m; col++) |
| 1248 | { |
| 1249 | if (col > row) break; |
| 1250 | MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1)); |
| 1251 | } |
| 1252 | for (int row = n + 2; row <= n + m; row++) |
| 1253 | for (int col = row - n; col <= m; col++) |
| 1254 | MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1)); |
| 1255 | for (int row = 2; row <= m + 1; row++) |
| 1256 | for (int col = m + 2; col <= m + n; col++) |
| 1257 | { |
| 1258 | if (col - m > row) break; |
| 1259 | MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1)); |
| 1260 | } |
| 1261 | for (int row = m + 2; row <= n + m; row++) |
| 1262 | for (int col = row; col <= m + n; col++) |
| 1263 | MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1)); |
| 1264 | |
| 1265 | /* constructing the LU-decomposition of A */ |
| 1266 | luDecomp(aMat, pMat, lMat, uMat); |
| 1267 | |
| 1268 | /* Before the xExp-th loop, we know that h = f*g mod <x^xExp>. |
| 1269 | Afterwards the algorithm ensures h = f*g mod <x^(xExp + 1)>. |
| 1270 | Hence in the end we obtain f and g as required, i.e., |
| 1271 | h = f*g mod <x^(d+1)>. |
| 1272 | The algorithm works by solving a (m+n)x(m+n) linear equation system |
| 1273 | A*x = b with constant matrix A (as decomposed above). By theory, the |
| 1274 | system is guaranteed to have a unique solution. */ |
| 1275 | poly fg = ppMult_qq(f, g); /* for storing the product of f and g */ |
| 1276 | for (int xExp = 1; xExp <= d; xExp++) |
no test coverage detected