[7bf145] | 1 | /*****************************************************************************\ |
---|
[806c18] | 2 | * Computer Algebra System SINGULAR |
---|
[7bf145] | 3 | \*****************************************************************************/ |
---|
| 4 | /** @file facFqBivarUtil.cc |
---|
[806c18] | 5 | * |
---|
| 6 | * This file provides utility functions for bivariate factorization |
---|
[7bf145] | 7 | * |
---|
| 8 | * @author Martin Lee |
---|
| 9 | * |
---|
| 10 | **/ |
---|
| 11 | /*****************************************************************************/ |
---|
| 12 | |
---|
[9f7665] | 13 | |
---|
[e4fe2b] | 14 | #include "config.h" |
---|
[9f7665] | 15 | |
---|
[7bf145] | 16 | |
---|
[5e4636] | 17 | #include "timing.h" |
---|
| 18 | |
---|
[7bf145] | 19 | #include "cf_map.h" |
---|
| 20 | #include "cf_map_ext.h" |
---|
[6db552] | 21 | #include "templates/ftmpl_functions.h" |
---|
[7bf145] | 22 | #include "ExtensionInfo.h" |
---|
[fecc08] | 23 | #include "cf_algorithm.h" |
---|
| 24 | #include "cf_factory.h" |
---|
[07718c3] | 25 | #include "cf_util.h" |
---|
[fecc08] | 26 | #include "imm.h" |
---|
| 27 | #include "cf_iter.h" |
---|
[7bf145] | 28 | #include "facFqBivarUtil.h" |
---|
[f876a66] | 29 | #include "cfNewtonPolygon.h" |
---|
| 30 | #include "facHensel.h" |
---|
[0e2e23] | 31 | #include "facMul.h" |
---|
[7bf145] | 32 | |
---|
[cb7b21] | 33 | #ifdef HAVE_FLINT |
---|
| 34 | #include "FLINTconvert.h" |
---|
| 35 | #endif |
---|
| 36 | |
---|
[5e4636] | 37 | TIMING_DEFINE_PRINT(fac_log_deriv_div) |
---|
| 38 | TIMING_DEFINE_PRINT(fac_log_deriv_mul) |
---|
| 39 | TIMING_DEFINE_PRINT(fac_log_deriv_pre) |
---|
[7bf145] | 40 | |
---|
[806c18] | 41 | void append (CFList& factors1, const CFList& factors2) |
---|
[7bf145] | 42 | { |
---|
[806c18] | 43 | for (CFListIterator i= factors2; i.hasItem(); i++) |
---|
[7bf145] | 44 | { |
---|
| 45 | if (!i.getItem().inCoeffDomain()) |
---|
| 46 | factors1.append (i.getItem()); |
---|
| 47 | } |
---|
| 48 | return; |
---|
| 49 | } |
---|
| 50 | |
---|
[806c18] | 51 | void decompress (CFList& factors, const CFMap& N) |
---|
[7bf145] | 52 | { |
---|
[806c18] | 53 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
[7bf145] | 54 | i.getItem()= N (i.getItem()); |
---|
[f876a66] | 55 | } |
---|
| 56 | |
---|
| 57 | void decompress (CFFList& factors, const CFMap& N) |
---|
| 58 | { |
---|
| 59 | for (CFFListIterator i= factors; i.hasItem(); i++) |
---|
| 60 | i.getItem()= CFFactor (N (i.getItem().factor()), i.getItem().exp()); |
---|
[7bf145] | 61 | } |
---|
| 62 | |
---|
[b728bf] | 63 | void decompress (CFAFList& factors, const CFMap& N) |
---|
| 64 | { |
---|
| 65 | for (CFAFListIterator i= factors; i.hasItem(); i++) |
---|
| 66 | i.getItem()= CFAFactor (N (i.getItem().factor()), i.getItem().minpoly(), |
---|
| 67 | i.getItem().exp()); |
---|
| 68 | } |
---|
| 69 | |
---|
[806c18] | 70 | void appendSwapDecompress (CFList& factors1, const CFList& factors2, |
---|
| 71 | const CFList& factors3, const bool swap1, |
---|
| 72 | const bool swap2, const CFMap& N) |
---|
[7bf145] | 73 | { |
---|
| 74 | Variable x= Variable (1); |
---|
| 75 | Variable y= Variable (2); |
---|
[806c18] | 76 | for (CFListIterator i= factors1; i.hasItem(); i++) |
---|
[7bf145] | 77 | { |
---|
[806c18] | 78 | if (swap1) |
---|
[7bf145] | 79 | { |
---|
[806c18] | 80 | if (!swap2) |
---|
[7bf145] | 81 | i.getItem()= swapvar (i.getItem(), x, y); |
---|
| 82 | } |
---|
[806c18] | 83 | else |
---|
[7bf145] | 84 | { |
---|
[806c18] | 85 | if (swap2) |
---|
[7bf145] | 86 | i.getItem()= swapvar (i.getItem(), y, x); |
---|
| 87 | } |
---|
| 88 | i.getItem()= N (i.getItem()); |
---|
| 89 | } |
---|
[806c18] | 90 | for (CFListIterator i= factors2; i.hasItem(); i++) |
---|
[7bf145] | 91 | factors1.append (N (i.getItem())); |
---|
[806c18] | 92 | for (CFListIterator i= factors3; i.hasItem(); i++) |
---|
[7bf145] | 93 | factors1.append (N (i.getItem())); |
---|
| 94 | return; |
---|
| 95 | } |
---|
| 96 | |
---|
[806c18] | 97 | void swapDecompress (CFList& factors, const bool swap, const CFMap& N) |
---|
[7bf145] | 98 | { |
---|
| 99 | Variable x= Variable (1); |
---|
| 100 | Variable y= Variable (2); |
---|
[806c18] | 101 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
[7bf145] | 102 | { |
---|
[806c18] | 103 | if (swap) |
---|
[7bf145] | 104 | i.getItem()= swapvar (i.getItem(), x, y); |
---|
| 105 | i.getItem()= N (i.getItem()); |
---|
| 106 | } |
---|
| 107 | return; |
---|
| 108 | } |
---|
| 109 | |
---|
[806c18] | 110 | static inline |
---|
| 111 | bool GFInExtensionHelper (const CanonicalForm& F, const int number) |
---|
[7bf145] | 112 | { |
---|
| 113 | if (F.isOne()) return false; |
---|
| 114 | InternalCF* buf; |
---|
| 115 | int exp; |
---|
| 116 | bool result= false; |
---|
[806c18] | 117 | if (F.inBaseDomain()) |
---|
[7bf145] | 118 | { |
---|
| 119 | buf= F.getval(); |
---|
| 120 | exp= imm2int (buf); |
---|
[806c18] | 121 | if (exp%number != 0) |
---|
[7bf145] | 122 | return true; |
---|
[806c18] | 123 | else |
---|
[7bf145] | 124 | return result; |
---|
| 125 | } |
---|
[806c18] | 126 | else |
---|
[7bf145] | 127 | { |
---|
[806c18] | 128 | for (CFIterator i= F; i.hasTerms(); i++) |
---|
[7bf145] | 129 | { |
---|
| 130 | result= GFInExtensionHelper (i.coeff(), number); |
---|
| 131 | if (result == true) |
---|
| 132 | return result; |
---|
| 133 | } |
---|
| 134 | } |
---|
| 135 | return result; |
---|
| 136 | } |
---|
| 137 | |
---|
[806c18] | 138 | static inline |
---|
[f876a66] | 139 | bool FqInExtensionHelper (const CanonicalForm& F, const CanonicalForm& gamma, |
---|
| 140 | const CanonicalForm& delta, CFList& source, |
---|
| 141 | CFList& dest) |
---|
[7bf145] | 142 | { |
---|
| 143 | bool result= false; |
---|
[806c18] | 144 | if (F.inBaseDomain()) |
---|
[7bf145] | 145 | return result; |
---|
[806c18] | 146 | else if (F.inCoeffDomain()) |
---|
[7bf145] | 147 | { |
---|
[806c18] | 148 | if (!fdivides (gamma, F)) |
---|
[7bf145] | 149 | return true; |
---|
[806c18] | 150 | else |
---|
[f876a66] | 151 | { |
---|
| 152 | int pos= findItem (source, F); |
---|
| 153 | if (pos > 0) |
---|
| 154 | return false; |
---|
| 155 | Variable a; |
---|
| 156 | hasFirstAlgVar (F, a); |
---|
| 157 | int bound= ipower (getCharacteristic(), degree (getMipo (a))); |
---|
| 158 | CanonicalForm buf= 1; |
---|
| 159 | for (int i= 1; i < bound; i++) |
---|
| 160 | { |
---|
| 161 | buf *= gamma; |
---|
| 162 | if (buf == F) |
---|
| 163 | { |
---|
| 164 | source.append (buf); |
---|
| 165 | dest.append (power (delta, i)); |
---|
| 166 | return false; |
---|
| 167 | } |
---|
| 168 | } |
---|
| 169 | return true; |
---|
| 170 | } |
---|
[7bf145] | 171 | } |
---|
[806c18] | 172 | else |
---|
[7bf145] | 173 | { |
---|
[806c18] | 174 | for (CFIterator i= F; i.hasTerms(); i++) |
---|
[7bf145] | 175 | { |
---|
[f876a66] | 176 | result= FqInExtensionHelper (i.coeff(), gamma, delta, source, dest); |
---|
[7bf145] | 177 | if (result == true) |
---|
[806c18] | 178 | return result; |
---|
[7bf145] | 179 | } |
---|
| 180 | } |
---|
| 181 | return result; |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | bool isInExtension (const CanonicalForm& F, const CanonicalForm& gamma, |
---|
[f876a66] | 185 | const int k, const CanonicalForm& delta, |
---|
| 186 | CFList& source, CFList& dest) |
---|
[7bf145] | 187 | { |
---|
| 188 | bool result; |
---|
[806c18] | 189 | if (CFFactory::gettype() == GaloisFieldDomain) |
---|
[7bf145] | 190 | { |
---|
| 191 | int p= getCharacteristic(); |
---|
| 192 | int orderFieldExtension= ipower (p, getGFDegree()) - 1; |
---|
| 193 | int order= ipower (p, k) - 1; |
---|
| 194 | int number= orderFieldExtension/order; |
---|
| 195 | result= GFInExtensionHelper (F, number); |
---|
| 196 | return result; |
---|
| 197 | } |
---|
[806c18] | 198 | else |
---|
[7bf145] | 199 | { |
---|
[f876a66] | 200 | result= FqInExtensionHelper (F, gamma, delta, source, dest); |
---|
[7bf145] | 201 | return result; |
---|
| 202 | } |
---|
| 203 | } |
---|
| 204 | |
---|
[806c18] | 205 | CanonicalForm |
---|
| 206 | mapDown (const CanonicalForm& F, const ExtensionInfo& info, CFList& source, |
---|
| 207 | CFList& dest) |
---|
[7bf145] | 208 | { |
---|
| 209 | int k= info.getGFDegree(); |
---|
| 210 | Variable beta= info.getAlpha(); |
---|
| 211 | CanonicalForm primElem= info.getGamma(); |
---|
| 212 | CanonicalForm imPrimElem= info.getDelta(); |
---|
| 213 | if (k > 1) |
---|
| 214 | return GFMapDown (F, k); |
---|
| 215 | else if (k == 1) |
---|
| 216 | return F; |
---|
[c1b9927] | 217 | if (/*k==0 &&*/ beta == Variable (1)) |
---|
[7bf145] | 218 | return F; |
---|
[c1b9927] | 219 | else /*if (k==0 && beta != Variable (1))*/ |
---|
[618da5] | 220 | return mapDown (F, imPrimElem, primElem, beta, source, dest); |
---|
[7bf145] | 221 | } |
---|
| 222 | |
---|
[806c18] | 223 | void appendTestMapDown (CFList& factors, const CanonicalForm& f, |
---|
| 224 | const ExtensionInfo& info, CFList& source, CFList& dest) |
---|
[7bf145] | 225 | { |
---|
| 226 | int k= info.getGFDegree(); |
---|
| 227 | Variable beta= info.getBeta(); |
---|
| 228 | Variable alpha= info.getAlpha(); |
---|
| 229 | CanonicalForm delta= info.getDelta(); |
---|
| 230 | CanonicalForm gamma= info.getGamma(); |
---|
| 231 | CanonicalForm g= f; |
---|
[88faa1] | 232 | int degMipoBeta=0; |
---|
[1372ae] | 233 | if (!k && beta.level() == 1) |
---|
[7bf145] | 234 | degMipoBeta= 1; |
---|
[1372ae] | 235 | else if (!k && beta.level() != 1) |
---|
[7bf145] | 236 | degMipoBeta= degree (getMipo (beta)); |
---|
[806c18] | 237 | if (k > 1) |
---|
[7bf145] | 238 | { |
---|
[f876a66] | 239 | if (!isInExtension (g, gamma, k, delta, source, dest)) |
---|
[7bf145] | 240 | { |
---|
| 241 | g= GFMapDown (g, k); |
---|
| 242 | factors.append (g); |
---|
| 243 | } |
---|
| 244 | } |
---|
[806c18] | 245 | else if (k == 1) |
---|
| 246 | { |
---|
[1709e1] | 247 | if (!isInExtension (g, gamma, k, delta, source, dest)) |
---|
[7bf145] | 248 | factors.append (g); |
---|
| 249 | } |
---|
[806c18] | 250 | else if (!k && beta == Variable (1)) |
---|
[7bf145] | 251 | { |
---|
| 252 | if (degree (g, alpha) < degMipoBeta) |
---|
| 253 | factors.append (g); |
---|
| 254 | } |
---|
[806c18] | 255 | else if (!k && beta != Variable (1)) |
---|
[7bf145] | 256 | { |
---|
[f876a66] | 257 | if (!isInExtension (g, gamma, k, delta, source, dest)) |
---|
[7bf145] | 258 | { |
---|
[618da5] | 259 | g= mapDown (g, delta, gamma, alpha, source, dest); |
---|
[7bf145] | 260 | factors.append (g); |
---|
[806c18] | 261 | } |
---|
[7bf145] | 262 | } |
---|
| 263 | return; |
---|
| 264 | } |
---|
| 265 | |
---|
[806c18] | 266 | void |
---|
| 267 | appendMapDown (CFList& factors, const CanonicalForm& g, |
---|
| 268 | const ExtensionInfo& info, CFList& source, CFList& dest) |
---|
[7bf145] | 269 | { |
---|
| 270 | int k= info.getGFDegree(); |
---|
| 271 | Variable beta= info.getBeta(); |
---|
| 272 | Variable alpha= info.getAlpha(); |
---|
| 273 | CanonicalForm gamma= info.getGamma(); |
---|
| 274 | CanonicalForm delta= info.getDelta(); |
---|
[806c18] | 275 | if (k > 1) |
---|
[7bf145] | 276 | factors.append (GFMapDown (g, k)); |
---|
[806c18] | 277 | else if (k == 1) |
---|
[7bf145] | 278 | factors.append (g); |
---|
[806c18] | 279 | else if (!k && beta == Variable (1)) |
---|
[7bf145] | 280 | factors.append (g); |
---|
[806c18] | 281 | else if (!k && beta != Variable (1)) |
---|
[618da5] | 282 | factors.append (mapDown (g, delta, gamma, alpha, source, dest)); |
---|
[7bf145] | 283 | return; |
---|
| 284 | } |
---|
| 285 | |
---|
[806c18] | 286 | void normalize (CFList& factors) |
---|
[7bf145] | 287 | { |
---|
[1682691] | 288 | CanonicalForm lcinv; |
---|
[806c18] | 289 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
[1682691] | 290 | { |
---|
| 291 | lcinv= 1/Lc (i.getItem()); |
---|
| 292 | i.getItem() *= lcinv; |
---|
| 293 | } |
---|
[7bf145] | 294 | return; |
---|
| 295 | } |
---|
| 296 | |
---|
[f876a66] | 297 | void normalize (CFFList& factors) |
---|
| 298 | { |
---|
[1682691] | 299 | CanonicalForm lcinv; |
---|
[c1b9927] | 300 | for (CFFListIterator i= factors; i.hasItem(); i++) |
---|
[1682691] | 301 | { |
---|
| 302 | lcinv= 1/ Lc (i.getItem().factor()); |
---|
| 303 | i.getItem()= CFFactor (i.getItem().factor()*lcinv, |
---|
[f876a66] | 304 | i.getItem().exp()); |
---|
[1682691] | 305 | } |
---|
[f876a66] | 306 | return; |
---|
| 307 | } |
---|
| 308 | |
---|
[806c18] | 309 | CFList subset (int index [], const int& s, const CFArray& elements, |
---|
| 310 | bool& noSubset) |
---|
[7bf145] | 311 | { |
---|
| 312 | int r= elements.size(); |
---|
| 313 | int i= 0; |
---|
| 314 | CFList result; |
---|
| 315 | noSubset= false; |
---|
[806c18] | 316 | if (index[s - 1] == 0) |
---|
[7bf145] | 317 | { |
---|
[806c18] | 318 | while (i < s) |
---|
[7bf145] | 319 | { |
---|
| 320 | index[i]= i + 1; |
---|
| 321 | result.append (elements[i]); |
---|
| 322 | i++; |
---|
[806c18] | 323 | } |
---|
[7bf145] | 324 | return result; |
---|
| 325 | } |
---|
| 326 | int buf; |
---|
| 327 | int k; |
---|
| 328 | bool found= false; |
---|
[806c18] | 329 | if (index[s - 1] == r) |
---|
[7bf145] | 330 | { |
---|
[806c18] | 331 | if (index[0] == r - s + 1) |
---|
[7bf145] | 332 | { |
---|
| 333 | noSubset= true; |
---|
| 334 | return result; |
---|
| 335 | } |
---|
| 336 | else { |
---|
[806c18] | 337 | while (found == false) |
---|
[7bf145] | 338 | { |
---|
| 339 | if (index[s - 2 - i] < r - i - 1) |
---|
| 340 | found= true; |
---|
[806c18] | 341 | i++; |
---|
[7bf145] | 342 | } |
---|
| 343 | buf= index[s - i - 1]; |
---|
| 344 | k= 0; |
---|
[806c18] | 345 | while (s - i - 1 + k < s) |
---|
[7bf145] | 346 | { |
---|
| 347 | index[s - i - 1 + k]= buf + k + 1; |
---|
| 348 | k++; |
---|
[806c18] | 349 | } |
---|
[7bf145] | 350 | } |
---|
[806c18] | 351 | for (int j= 0; j < s; j++) |
---|
[7bf145] | 352 | result.append (elements[index[j] - 1]); |
---|
| 353 | return result; |
---|
| 354 | } |
---|
[806c18] | 355 | else |
---|
[7bf145] | 356 | { |
---|
| 357 | index[s - 1] += 1; |
---|
| 358 | for (int j= 0; j < s; j++) |
---|
| 359 | result.append (elements[index[j] - 1]); |
---|
| 360 | return result; |
---|
| 361 | } |
---|
[806c18] | 362 | } |
---|
[7bf145] | 363 | |
---|
[806c18] | 364 | CFArray copy (const CFList& list) |
---|
[7bf145] | 365 | { |
---|
| 366 | CFArray array= CFArray (list.length()); |
---|
| 367 | int j= 0; |
---|
[806c18] | 368 | for (CFListIterator i= list; i.hasItem(); i++, j++) |
---|
[7bf145] | 369 | array[j]= i.getItem(); |
---|
| 370 | return array; |
---|
| 371 | } |
---|
| 372 | |
---|
[806c18] | 373 | void indexUpdate (int index [], const int& subsetSize, const int& setSize, |
---|
| 374 | bool& noSubset) |
---|
[7bf145] | 375 | { |
---|
| 376 | noSubset= false; |
---|
[806c18] | 377 | if (subsetSize > setSize) |
---|
[7bf145] | 378 | { |
---|
| 379 | noSubset= true; |
---|
| 380 | return; |
---|
| 381 | } |
---|
[1673386] | 382 | int * v= new int [setSize]; |
---|
[806c18] | 383 | for (int i= 0; i < setSize; i++) |
---|
[7bf145] | 384 | v[i]= index[i]; |
---|
[806c18] | 385 | if (subsetSize == 1) |
---|
[7bf145] | 386 | { |
---|
| 387 | v[0]= v[0] - 1; |
---|
[806c18] | 388 | if (v[0] >= setSize) |
---|
[7bf145] | 389 | { |
---|
| 390 | noSubset= true; |
---|
[1673386] | 391 | delete [] v; |
---|
[7bf145] | 392 | return; |
---|
| 393 | } |
---|
| 394 | } |
---|
[806c18] | 395 | else |
---|
[7bf145] | 396 | { |
---|
[806c18] | 397 | if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1) |
---|
[7bf145] | 398 | { |
---|
[806c18] | 399 | if (v[0] + subsetSize - 1 > setSize) |
---|
| 400 | { |
---|
[7bf145] | 401 | noSubset= true; |
---|
[1673386] | 402 | delete [] v; |
---|
[7bf145] | 403 | return; |
---|
| 404 | } |
---|
| 405 | v[0]= v[0] - 1; |
---|
[806c18] | 406 | for (int i= 1; i < subsetSize - 1; i++) |
---|
[7bf145] | 407 | v[i]= v[i - 1] + 1; |
---|
| 408 | v[subsetSize - 1]= v[subsetSize - 2]; |
---|
| 409 | } |
---|
[806c18] | 410 | else |
---|
[7bf145] | 411 | { |
---|
[806c18] | 412 | if (v[0] + subsetSize - 1 > setSize) |
---|
| 413 | { |
---|
[7bf145] | 414 | noSubset= true; |
---|
[1673386] | 415 | delete [] v; |
---|
[7bf145] | 416 | return; |
---|
| 417 | } |
---|
| 418 | for (int i= 1; i < subsetSize - 1; i++) |
---|
| 419 | v[i]= v[i - 1] + 1; |
---|
| 420 | v[subsetSize - 1]= v[subsetSize - 2]; |
---|
| 421 | } |
---|
| 422 | } |
---|
| 423 | for (int i= 0; i < setSize; i++) |
---|
| 424 | index[i]= v[i]; |
---|
[1673386] | 425 | delete [] v; |
---|
[7bf145] | 426 | } |
---|
| 427 | |
---|
[806c18] | 428 | int subsetDegree (const CFList& S) |
---|
[7bf145] | 429 | { |
---|
| 430 | int result= 0; |
---|
| 431 | for (CFListIterator i= S; i.hasItem(); i++) |
---|
| 432 | result += degree (i.getItem(), Variable (1)); |
---|
| 433 | return result; |
---|
| 434 | } |
---|
| 435 | |
---|
| 436 | CFFList multiplicity (CanonicalForm& F, const CFList& factors) |
---|
| 437 | { |
---|
| 438 | if (F.inCoeffDomain()) |
---|
| 439 | return CFFList (CFFactor (F, 1)); |
---|
| 440 | CFFList result; |
---|
| 441 | int multi= 0; |
---|
[21b8f4c] | 442 | CanonicalForm quot; |
---|
[7bf145] | 443 | for (CFListIterator i= factors; i.hasItem(); i++) |
---|
| 444 | { |
---|
[21b8f4c] | 445 | while (fdivides (i.getItem(), F, quot)) |
---|
[7bf145] | 446 | { |
---|
| 447 | multi++; |
---|
[21b8f4c] | 448 | F= quot; |
---|
[7bf145] | 449 | } |
---|
[575e1d] | 450 | if (multi > 0) |
---|
| 451 | result.append (CFFactor (i.getItem(), multi)); |
---|
[7bf145] | 452 | multi= 0; |
---|
| 453 | } |
---|
| 454 | return result; |
---|
| 455 | } |
---|
| 456 | |
---|
[615ca8] | 457 | #ifdef HAVE_NTL |
---|
[f876a66] | 458 | CFArray |
---|
| 459 | logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l, |
---|
| 460 | CanonicalForm& Q |
---|
| 461 | ) |
---|
| 462 | { |
---|
| 463 | Variable x= Variable (2); |
---|
| 464 | Variable y= Variable (1); |
---|
| 465 | CanonicalForm xToL= power (x, l); |
---|
| 466 | CanonicalForm q,r; |
---|
| 467 | CanonicalForm logDeriv; |
---|
| 468 | |
---|
[5e4636] | 469 | TIMING_START (fac_log_deriv_div); |
---|
[f876a66] | 470 | q= newtonDiv (F, G, xToL); |
---|
[5e4636] | 471 | TIMING_END_AND_PRINT (fac_log_deriv_div, "time for division in logderiv1: "); |
---|
[f876a66] | 472 | |
---|
[5e4636] | 473 | TIMING_START (fac_log_deriv_mul); |
---|
[f876a66] | 474 | logDeriv= mulMod2 (q, deriv (G, y), xToL); |
---|
[5e4636] | 475 | TIMING_END_AND_PRINT (fac_log_deriv_mul, "time to multiply in logderiv1: "); |
---|
[f876a66] | 476 | |
---|
[a3b213] | 477 | if (degree (logDeriv, x) == 0) |
---|
| 478 | { |
---|
| 479 | Q= q; |
---|
| 480 | return CFArray(); |
---|
| 481 | } |
---|
| 482 | |
---|
[9614bb] | 483 | int j= degree (logDeriv, y) + 1; |
---|
[f876a66] | 484 | CFArray result= CFArray (j); |
---|
[9614bb] | 485 | CFIterator ii; |
---|
[f876a66] | 486 | for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++) |
---|
| 487 | { |
---|
[0b447e] | 488 | if (i.coeff().inCoeffDomain()) |
---|
| 489 | result[0] += i.coeff()*power (x,i.exp()); |
---|
| 490 | else |
---|
| 491 | { |
---|
| 492 | for (ii= i.coeff(); ii.hasTerms(); ii++) |
---|
| 493 | result[ii.exp()] += ii.coeff()*power (x,i.exp()); |
---|
| 494 | } |
---|
[f876a66] | 495 | } |
---|
| 496 | Q= q; |
---|
| 497 | return result; |
---|
| 498 | } |
---|
| 499 | |
---|
| 500 | CFArray |
---|
| 501 | logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l, |
---|
| 502 | int oldL, const CanonicalForm& oldQ, CanonicalForm& Q |
---|
| 503 | ) |
---|
| 504 | { |
---|
| 505 | Variable x= Variable (2); |
---|
| 506 | Variable y= Variable (1); |
---|
| 507 | CanonicalForm xToL= power (x, l); |
---|
| 508 | CanonicalForm xToOldL= power (x, oldL); |
---|
| 509 | CanonicalForm xToLOldL= power (x, l-oldL); |
---|
| 510 | CanonicalForm q,r; |
---|
| 511 | CanonicalForm logDeriv; |
---|
| 512 | |
---|
[ac3fcca] | 513 | CanonicalForm bufF; |
---|
[5e4636] | 514 | TIMING_START (fac_log_deriv_pre); |
---|
[ac3fcca] | 515 | if ((oldL > 100 && l - oldL < 50) || (oldL < 100 && l - oldL < 30)) |
---|
| 516 | { |
---|
| 517 | bufF= F; |
---|
| 518 | CanonicalForm oldF= mulMod2 (G, oldQ, xToL); |
---|
| 519 | bufF -= oldF; |
---|
| 520 | bufF= div (bufF, xToOldL); |
---|
| 521 | } |
---|
| 522 | else |
---|
| 523 | { |
---|
| 524 | //middle product style computation of [G*oldQ]^{l}_{oldL} |
---|
| 525 | CanonicalForm G3= div (G, xToOldL); |
---|
| 526 | CanonicalForm Up= mulMod2 (G3, oldQ, xToLOldL); |
---|
[fd80670] | 527 | CanonicalForm xToOldL2= power (x, (oldL+1)/2); |
---|
[ac3fcca] | 528 | CanonicalForm G2= mod (G, xToOldL); |
---|
| 529 | CanonicalForm G1= div (G2, xToOldL2); |
---|
| 530 | CanonicalForm G0= mod (G2, xToOldL2); |
---|
| 531 | CanonicalForm oldQ1= div (oldQ, xToOldL2); |
---|
| 532 | CanonicalForm oldQ0= mod (oldQ, xToOldL2); |
---|
[fd80670] | 533 | CanonicalForm Mid; |
---|
| 534 | if (oldL % 2 == 1) |
---|
| 535 | Mid= mulMod2 (G1, oldQ1*x, xToLOldL); |
---|
| 536 | else |
---|
| 537 | Mid= mulMod2 (G1, oldQ1, xToLOldL); |
---|
[ac3fcca] | 538 | //computation of Low might be faster using a real middle product? |
---|
| 539 | CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL); |
---|
[fd80670] | 540 | Low= div (Low, power (x, oldL/2)); |
---|
| 541 | Low= mod (Low, xToLOldL); |
---|
[ac3fcca] | 542 | Up += Mid + Low; |
---|
| 543 | bufF= div (F, xToOldL); |
---|
| 544 | bufF -= Up; |
---|
| 545 | } |
---|
[5e4636] | 546 | TIMING_END_AND_PRINT (fac_log_deriv_pre, "time to preprocess: "); |
---|
[f876a66] | 547 | |
---|
[5e4636] | 548 | TIMING_START (fac_log_deriv_div); |
---|
[72f1e4b] | 549 | if (l-oldL > 0) |
---|
| 550 | q= newtonDiv (bufF, G, xToLOldL); |
---|
| 551 | else |
---|
| 552 | q= 0; |
---|
[f876a66] | 553 | q *= xToOldL; |
---|
| 554 | q += oldQ; |
---|
[5e4636] | 555 | TIMING_END_AND_PRINT (fac_log_deriv_div, "time for div in logderiv2: "); |
---|
[f876a66] | 556 | |
---|
[5e4636] | 557 | TIMING_START (fac_log_deriv_mul); |
---|
[f876a66] | 558 | logDeriv= mulMod2 (q, deriv (G, y), xToL); |
---|
[5e4636] | 559 | TIMING_END_AND_PRINT (fac_log_deriv_mul, "time for mul in logderiv2: "); |
---|
[f876a66] | 560 | |
---|
[a3b213] | 561 | if (degree (logDeriv, x) == 0) |
---|
| 562 | { |
---|
| 563 | Q= q; |
---|
| 564 | return CFArray(); |
---|
| 565 | } |
---|
| 566 | |
---|
[9614bb] | 567 | int j= degree (logDeriv,y) + 1; |
---|
[f876a66] | 568 | CFArray result= CFArray (j); |
---|
[9614bb] | 569 | CFIterator ii; |
---|
[f876a66] | 570 | for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++) |
---|
| 571 | { |
---|
[0b447e] | 572 | if (i.coeff().inCoeffDomain()) |
---|
| 573 | result[0] += i.coeff()*power (x,i.exp()); |
---|
| 574 | else |
---|
| 575 | { |
---|
| 576 | for (ii= i.coeff(); ii.hasTerms(); ii++) |
---|
| 577 | result[ii.exp()] += ii.coeff()*power (x,i.exp()); |
---|
| 578 | } |
---|
[f876a66] | 579 | } |
---|
| 580 | Q= q; |
---|
| 581 | return result; |
---|
| 582 | } |
---|
[615ca8] | 583 | #endif |
---|
[f876a66] | 584 | |
---|
| 585 | void |
---|
| 586 | writeInMatrix (CFMatrix& M, const CFArray& A, const int column, |
---|
| 587 | const int startIndex |
---|
| 588 | ) |
---|
| 589 | { |
---|
[050d1b] | 590 | ASSERT (A.size () - startIndex >= 0, "wrong starting index"); |
---|
| 591 | ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index"); |
---|
[f876a66] | 592 | ASSERT (column > 0 && column <= M.columns(), "wrong column"); |
---|
| 593 | if (A.size() - startIndex <= 0) return; |
---|
| 594 | int j= 1; |
---|
| 595 | for (int i= startIndex; i < A.size(); i++, j++) |
---|
| 596 | M (j, column)= A [i]; |
---|
| 597 | } |
---|
| 598 | |
---|
| 599 | CFArray getCoeffs (const CanonicalForm& F, const int k) |
---|
| 600 | { |
---|
[050d1b] | 601 | ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected"); |
---|
[f876a66] | 602 | if (degree (F, 2) < k) |
---|
| 603 | return CFArray(); |
---|
| 604 | |
---|
| 605 | CFArray result= CFArray (degree (F) - k + 1); |
---|
| 606 | CFIterator j= F; |
---|
| 607 | for (int i= degree (F); i >= k; i--) |
---|
| 608 | { |
---|
| 609 | if (j.exp() == i) |
---|
| 610 | { |
---|
| 611 | result [i - k]= j.coeff(); |
---|
[c1b9927] | 612 | j++; |
---|
[f876a66] | 613 | if (!j.hasTerms()) |
---|
| 614 | return result; |
---|
| 615 | } |
---|
| 616 | else |
---|
| 617 | result[i - k]= 0; |
---|
| 618 | } |
---|
| 619 | return result; |
---|
| 620 | } |
---|
| 621 | |
---|
| 622 | CFArray getCoeffs (const CanonicalForm& F, const int k, const Variable& alpha) |
---|
| 623 | { |
---|
[050d1b] | 624 | ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected"); |
---|
[f876a66] | 625 | if (degree (F, 2) < k) |
---|
| 626 | return CFArray (); |
---|
| 627 | |
---|
| 628 | int d= degree (getMipo (alpha)); |
---|
| 629 | CFArray result= CFArray ((degree (F) - k + 1)*d); |
---|
| 630 | CFIterator j= F; |
---|
| 631 | CanonicalForm buf; |
---|
| 632 | CFIterator iter; |
---|
| 633 | for (int i= degree (F); i >= k; i--) |
---|
| 634 | { |
---|
| 635 | if (j.exp() == i) |
---|
| 636 | { |
---|
| 637 | iter= j.coeff(); |
---|
| 638 | for (int l= degree (j.coeff(), alpha); l >= 0; l--) |
---|
| 639 | { |
---|
| 640 | if (iter.exp() == l) |
---|
| 641 | { |
---|
| 642 | result [(i - k)*d + l]= iter.coeff(); |
---|
| 643 | iter++; |
---|
| 644 | if (!iter.hasTerms()) |
---|
| 645 | break; |
---|
| 646 | } |
---|
| 647 | } |
---|
| 648 | j++; |
---|
| 649 | if (!j.hasTerms()) |
---|
| 650 | return result; |
---|
| 651 | } |
---|
| 652 | else |
---|
| 653 | { |
---|
| 654 | for (int l= 0; l < d; l++) |
---|
| 655 | result[(i - k)*d + l]= 0; |
---|
| 656 | } |
---|
| 657 | } |
---|
| 658 | return result; |
---|
| 659 | } |
---|
| 660 | |
---|
| 661 | #ifdef HAVE_NTL |
---|
| 662 | CFArray |
---|
| 663 | getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo, |
---|
| 664 | const Variable& alpha, const CanonicalForm& evaluation, |
---|
| 665 | const mat_zz_p& M) |
---|
| 666 | { |
---|
[050d1b] | 667 | ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected"); |
---|
[f876a66] | 668 | CanonicalForm F= G (G.mvar() - evaluation, G.mvar()); |
---|
| 669 | if (F.isZero()) |
---|
| 670 | return CFArray (); |
---|
| 671 | |
---|
| 672 | Variable y= Variable (2); |
---|
| 673 | F= F (power (y, degMipo), y); |
---|
| 674 | F= F (y, alpha); |
---|
| 675 | zz_pX NTLF= convertFacCF2NTLzzpX (F); |
---|
| 676 | NTLF.rep.SetLength (l*degMipo); |
---|
| 677 | NTLF.rep= M*NTLF.rep; |
---|
| 678 | NTLF.normalize(); |
---|
| 679 | F= convertNTLzzpX2CF (NTLF, y); |
---|
| 680 | |
---|
| 681 | if (degree (F, 2) < k) |
---|
| 682 | return CFArray(); |
---|
[cb7b21] | 683 | |
---|
| 684 | CFArray result= CFArray (degree (F) - k + 1); |
---|
| 685 | |
---|
| 686 | CFIterator j= F; |
---|
| 687 | for (int i= degree (F); i >= k; i--) |
---|
| 688 | { |
---|
| 689 | if (j.exp() == i) |
---|
| 690 | { |
---|
| 691 | result [i - k]= j.coeff(); |
---|
| 692 | j++; |
---|
| 693 | if (!j.hasTerms()) |
---|
| 694 | return result; |
---|
| 695 | } |
---|
| 696 | else |
---|
| 697 | result[i - k]= 0; |
---|
| 698 | } |
---|
| 699 | return result; |
---|
| 700 | } |
---|
| 701 | #endif |
---|
| 702 | |
---|
| 703 | #ifdef HAVE_FLINT |
---|
| 704 | CFArray |
---|
| 705 | getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo, |
---|
| 706 | const Variable& alpha, const CanonicalForm& evaluation, |
---|
| 707 | const nmod_mat_t M) |
---|
| 708 | { |
---|
| 709 | ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected"); |
---|
| 710 | CanonicalForm F= G (G.mvar() - evaluation, G.mvar()); |
---|
| 711 | if (F.isZero()) |
---|
| 712 | return CFArray (); |
---|
| 713 | |
---|
| 714 | Variable y= Variable (2); |
---|
| 715 | F= F (power (y, degMipo), y); |
---|
| 716 | F= F (y, alpha); |
---|
| 717 | |
---|
| 718 | nmod_poly_t FLINTF; |
---|
| 719 | nmod_mat_t MFLINTF, mulResult; |
---|
| 720 | nmod_mat_init (MFLINTF, l*degMipo, 1, getCharacteristic()); |
---|
| 721 | nmod_mat_init (mulResult, l*degMipo, 1, getCharacteristic()); |
---|
| 722 | |
---|
| 723 | convertFacCF2nmod_poly_t (FLINTF, F); |
---|
| 724 | |
---|
[6f1268] | 725 | #ifndef slong |
---|
| 726 | #define slong long |
---|
| 727 | #endif |
---|
[cb7b21] | 728 | slong i; |
---|
| 729 | |
---|
| 730 | for (i= 0; i < FLINTF->length; i++) |
---|
| 731 | nmod_mat_entry (MFLINTF, i, 0)= FLINTF->coeffs[i]; |
---|
| 732 | |
---|
| 733 | for (; i < MFLINTF->r; i++) |
---|
[6f1268] | 734 | nmod_mat_entry (MFLINTF, i, 0)= 0; |
---|
[cb7b21] | 735 | |
---|
| 736 | nmod_mat_mul (mulResult, M, MFLINTF); |
---|
| 737 | |
---|
| 738 | F= 0; |
---|
| 739 | for (i= 0; i < mulResult->r; i++) |
---|
| 740 | F += CanonicalForm ((long) nmod_mat_entry (mulResult, i, 0))*power (y, i); |
---|
| 741 | |
---|
| 742 | nmod_mat_clear (MFLINTF); |
---|
| 743 | nmod_mat_clear (mulResult); |
---|
[8a0e72] | 744 | nmod_poly_clear(FLINTF); |
---|
[cb7b21] | 745 | |
---|
| 746 | if (degree (F, 2) < k) |
---|
| 747 | return CFArray(); |
---|
[f876a66] | 748 | |
---|
| 749 | CFArray result= CFArray (degree (F) - k + 1); |
---|
| 750 | |
---|
| 751 | CFIterator j= F; |
---|
| 752 | for (int i= degree (F); i >= k; i--) |
---|
| 753 | { |
---|
| 754 | if (j.exp() == i) |
---|
| 755 | { |
---|
| 756 | result [i - k]= j.coeff(); |
---|
[c1b9927] | 757 | j++; |
---|
[f876a66] | 758 | if (!j.hasTerms()) |
---|
| 759 | return result; |
---|
| 760 | } |
---|
| 761 | else |
---|
| 762 | result[i - k]= 0; |
---|
| 763 | } |
---|
| 764 | return result; |
---|
| 765 | } |
---|
| 766 | #endif |
---|
| 767 | |
---|
[542864] | 768 | int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible) |
---|
[f876a66] | 769 | { |
---|
| 770 | n= degree (F, 1); |
---|
| 771 | int* result= new int [n]; |
---|
| 772 | int sizeOfNewtonPolygon; |
---|
| 773 | int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon); |
---|
| 774 | |
---|
[542864] | 775 | isIrreducible= false; |
---|
| 776 | if (sizeOfNewtonPolygon == 3) |
---|
| 777 | { |
---|
| 778 | bool check1= |
---|
| 779 | (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0); |
---|
| 780 | if (check1) |
---|
| 781 | { |
---|
| 782 | bool check2= |
---|
| 783 | (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0); |
---|
| 784 | if (check2) |
---|
| 785 | { |
---|
| 786 | int p=getCharacteristic(); |
---|
| 787 | int d=1; |
---|
| 788 | char bufGFName='Z'; |
---|
| 789 | bool GF= (CFFactory::gettype()==GaloisFieldDomain); |
---|
| 790 | if (GF) |
---|
| 791 | { |
---|
| 792 | d= getGFDegree(); |
---|
| 793 | bufGFName=gf_name; |
---|
| 794 | } |
---|
| 795 | setCharacteristic(0); |
---|
[36ef97a] | 796 | CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]); // maybe it's better to use plain intgcd |
---|
[542864] | 797 | tmp= gcd (tmp, newtonPolyg[1][0]); |
---|
| 798 | tmp= gcd (tmp, newtonPolyg[1][1]); |
---|
| 799 | tmp= gcd (tmp, newtonPolyg[2][0]); |
---|
| 800 | tmp= gcd (tmp, newtonPolyg[2][1]); |
---|
| 801 | isIrreducible= (tmp==1); |
---|
| 802 | if (GF) |
---|
| 803 | setCharacteristic (p, d, bufGFName); |
---|
| 804 | else |
---|
| 805 | setCharacteristic(p); |
---|
| 806 | } |
---|
| 807 | } |
---|
| 808 | } |
---|
| 809 | |
---|
[f876a66] | 810 | int minX, minY, maxX, maxY; |
---|
| 811 | minX= newtonPolyg [0] [0]; |
---|
| 812 | minY= newtonPolyg [0] [1]; |
---|
| 813 | maxX= minX; |
---|
| 814 | maxY= minY; |
---|
[f0befc] | 815 | int indZero= 0; |
---|
[f876a66] | 816 | for (int i= 1; i < sizeOfNewtonPolygon; i++) |
---|
| 817 | { |
---|
[f0befc] | 818 | if (newtonPolyg[i][1] == 0) |
---|
| 819 | { |
---|
| 820 | if (newtonPolyg[indZero][1] == 0) |
---|
| 821 | { |
---|
| 822 | if (newtonPolyg[indZero][0] < newtonPolyg[i][0]) |
---|
| 823 | indZero= i; |
---|
| 824 | } |
---|
| 825 | else |
---|
| 826 | indZero= i; |
---|
| 827 | } |
---|
[f876a66] | 828 | if (minX > newtonPolyg [i] [0]) |
---|
| 829 | minX= newtonPolyg [i] [0]; |
---|
| 830 | if (maxX < newtonPolyg [i] [0]) |
---|
| 831 | maxX= newtonPolyg [i] [0]; |
---|
| 832 | if (minY > newtonPolyg [i] [1]) |
---|
| 833 | minY= newtonPolyg [i] [1]; |
---|
| 834 | if (maxY < newtonPolyg [i] [1]) |
---|
| 835 | maxY= newtonPolyg [i] [1]; |
---|
| 836 | } |
---|
| 837 | |
---|
[f0befc] | 838 | int slopeNum, slopeDen, constTerm; |
---|
| 839 | bool negativeSlope=false; |
---|
| 840 | if (indZero != sizeOfNewtonPolygon - 1) |
---|
| 841 | { |
---|
| 842 | slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0]; |
---|
| 843 | slopeDen= newtonPolyg[indZero+1][1]; |
---|
| 844 | constTerm= newtonPolyg[indZero][0]; |
---|
| 845 | } |
---|
| 846 | else |
---|
| 847 | { |
---|
| 848 | slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0]; |
---|
| 849 | slopeDen= newtonPolyg[0][1]; |
---|
| 850 | constTerm= newtonPolyg[indZero][0]; |
---|
| 851 | } |
---|
| 852 | if (slopeNum < 0) |
---|
| 853 | { |
---|
| 854 | slopeNum= -slopeNum; |
---|
| 855 | negativeSlope= true; |
---|
| 856 | } |
---|
| 857 | int k= 0; |
---|
[9fbf02] | 858 | int* point= new int [2]; |
---|
[f876a66] | 859 | for (int i= 0; i < n; i++) |
---|
| 860 | { |
---|
[f0befc] | 861 | if (((indZero+1) < sizeOfNewtonPolygon && (i+1) > newtonPolyg[indZero+1][1]) |
---|
| 862 | || ((indZero+1) >= sizeOfNewtonPolygon && (i+1) > newtonPolyg[0][1])) |
---|
| 863 | { |
---|
| 864 | if (indZero + 1 != sizeOfNewtonPolygon) |
---|
| 865 | indZero++; |
---|
| 866 | else |
---|
| 867 | indZero= 0; |
---|
| 868 | if (indZero != sizeOfNewtonPolygon - 1) |
---|
| 869 | { |
---|
| 870 | slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0]; |
---|
| 871 | slopeDen= newtonPolyg[indZero+1][1]-newtonPolyg[indZero][1]; |
---|
| 872 | constTerm= newtonPolyg[indZero][0]; |
---|
| 873 | } |
---|
| 874 | else |
---|
| 875 | { |
---|
| 876 | slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0]; |
---|
| 877 | slopeDen= newtonPolyg[0][1]-newtonPolyg[indZero][1]; |
---|
| 878 | constTerm= newtonPolyg[indZero][0]; |
---|
| 879 | } |
---|
| 880 | if (slopeNum < 0) |
---|
| 881 | { |
---|
| 882 | negativeSlope= true; |
---|
| 883 | slopeNum= - slopeNum; |
---|
| 884 | k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/ |
---|
| 885 | slopeDen) + constTerm; |
---|
| 886 | } |
---|
| 887 | else |
---|
| 888 | k= (int) (((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen) |
---|
| 889 | + constTerm; |
---|
| 890 | } |
---|
| 891 | else |
---|
| 892 | { |
---|
| 893 | if (negativeSlope) |
---|
| 894 | k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/ |
---|
| 895 | slopeDen) + constTerm; |
---|
| 896 | else |
---|
| 897 | k= (int) ((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen |
---|
| 898 | + constTerm; |
---|
| 899 | } |
---|
[f876a66] | 900 | if (i + 1 > maxY || i + 1 < minY) |
---|
| 901 | { |
---|
| 902 | result [i]= 0; |
---|
| 903 | continue; |
---|
| 904 | } |
---|
| 905 | point [0]= k; |
---|
| 906 | point [1]= i + 1; |
---|
[f0befc] | 907 | if (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0) |
---|
| 908 | k= 0; |
---|
[f876a66] | 909 | result [i]= k; |
---|
| 910 | } |
---|
| 911 | |
---|
[9fbf02] | 912 | delete [] point; |
---|
| 913 | |
---|
[f0befc] | 914 | for (int i= 0; i < sizeOfNewtonPolygon; i++) |
---|
| 915 | delete [] newtonPolyg[i]; |
---|
| 916 | delete [] newtonPolyg; |
---|
| 917 | |
---|
[f876a66] | 918 | return result; |
---|
| 919 | } |
---|
| 920 | |
---|
[d8a7da] | 921 | int * |
---|
| 922 | computeBoundsWrtDiffMainvar (const CanonicalForm& F, int& n, |
---|
| 923 | bool& isIrreducible) |
---|
| 924 | { |
---|
| 925 | n= degree (F, 2); |
---|
| 926 | int* result= new int [n]; |
---|
| 927 | int sizeOfNewtonPolygon; |
---|
| 928 | int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon); |
---|
| 929 | |
---|
| 930 | isIrreducible= false; |
---|
| 931 | if (sizeOfNewtonPolygon == 3) |
---|
| 932 | { |
---|
| 933 | bool check1= |
---|
| 934 | (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0); |
---|
| 935 | if (check1) |
---|
| 936 | { |
---|
| 937 | bool check2= |
---|
| 938 | (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0); |
---|
| 939 | if (check2) |
---|
| 940 | { |
---|
| 941 | int p=getCharacteristic(); |
---|
| 942 | int d=1; |
---|
| 943 | char bufGFName='Z'; |
---|
| 944 | bool GF= (CFFactory::gettype()==GaloisFieldDomain); |
---|
| 945 | if (GF) |
---|
| 946 | { |
---|
| 947 | d= getGFDegree(); |
---|
| 948 | bufGFName=gf_name; |
---|
| 949 | } |
---|
| 950 | setCharacteristic(0); |
---|
| 951 | CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]); |
---|
| 952 | tmp= gcd (tmp, newtonPolyg[1][0]); |
---|
| 953 | tmp= gcd (tmp, newtonPolyg[1][1]); |
---|
| 954 | tmp= gcd (tmp, newtonPolyg[2][0]); |
---|
| 955 | tmp= gcd (tmp, newtonPolyg[2][1]); |
---|
| 956 | isIrreducible= (tmp==1); |
---|
| 957 | if (GF) |
---|
| 958 | setCharacteristic (p, d, bufGFName); |
---|
| 959 | else |
---|
| 960 | setCharacteristic(p); |
---|
| 961 | } |
---|
| 962 | } |
---|
| 963 | } |
---|
| 964 | |
---|
[9fbf02] | 965 | int swap; |
---|
| 966 | for (int i= 0; i < sizeOfNewtonPolygon; i++) |
---|
| 967 | { |
---|
| 968 | swap= newtonPolyg[i][1]; |
---|
| 969 | newtonPolyg[i][1]=newtonPolyg[i][0]; |
---|
| 970 | newtonPolyg[i][0]= swap; |
---|
| 971 | } |
---|
| 972 | |
---|
| 973 | sizeOfNewtonPolygon= polygon(newtonPolyg, sizeOfNewtonPolygon); |
---|
| 974 | |
---|
[d8a7da] | 975 | int minX, minY, maxX, maxY; |
---|
| 976 | minX= newtonPolyg [0] [0]; |
---|
| 977 | minY= newtonPolyg [0] [1]; |
---|
| 978 | maxX= minX; |
---|
| 979 | maxY= minY; |
---|
[f0befc] | 980 | int indZero= 0; |
---|
[d8a7da] | 981 | for (int i= 1; i < sizeOfNewtonPolygon; i++) |
---|
| 982 | { |
---|
[9fbf02] | 983 | if (newtonPolyg[i][1] == 0) |
---|
[f0befc] | 984 | { |
---|
[9fbf02] | 985 | if (newtonPolyg[indZero][1] == 0) |
---|
[f0befc] | 986 | { |
---|
[9fbf02] | 987 | if (newtonPolyg[indZero][0] < newtonPolyg[i][0]) |
---|
[f0befc] | 988 | indZero= i; |
---|
| 989 | } |
---|
| 990 | else |
---|
| 991 | indZero= i; |
---|
| 992 | } |
---|
[d8a7da] | 993 | if (minX > newtonPolyg [i] [0]) |
---|
| 994 | minX= newtonPolyg [i] [0]; |
---|
| 995 | if (maxX < newtonPolyg [i] [0]) |
---|
| 996 | maxX= newtonPolyg [i] [0]; |
---|
| 997 | if (minY > newtonPolyg [i] [1]) |
---|
| 998 | minY= newtonPolyg [i] [1]; |
---|
| 999 | if (maxY < newtonPolyg [i] [1]) |
---|
| 1000 | maxY= newtonPolyg [i] [1]; |
---|
| 1001 | } |
---|
| 1002 | |
---|
[f0befc] | 1003 | int slopeNum, slopeDen, constTerm; |
---|
| 1004 | bool negativeSlope=false; |
---|
| 1005 | if (indZero != sizeOfNewtonPolygon - 1) |
---|
| 1006 | { |
---|
[9fbf02] | 1007 | slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0]; |
---|
| 1008 | slopeDen= newtonPolyg[indZero+1][1]; |
---|
| 1009 | constTerm= newtonPolyg[indZero][0]; |
---|
[f0befc] | 1010 | } |
---|
| 1011 | else |
---|
| 1012 | { |
---|
[9fbf02] | 1013 | slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0]; |
---|
| 1014 | slopeDen= newtonPolyg[0][1]; |
---|
| 1015 | constTerm= newtonPolyg[indZero][0]; |
---|
[f0befc] | 1016 | } |
---|
| 1017 | if (slopeNum < 0) |
---|
| 1018 | { |
---|
| 1019 | slopeNum= -slopeNum; |
---|
| 1020 | negativeSlope= true; |
---|
| 1021 | } |
---|
| 1022 | int k= 0; |
---|
[9fbf02] | 1023 | |
---|
| 1024 | int* point= new int [2]; |
---|
[d8a7da] | 1025 | for (int i= 0; i < n; i++) |
---|
| 1026 | { |
---|
[9fbf02] | 1027 | if (((indZero+1) < sizeOfNewtonPolygon && (i+1) > newtonPolyg[indZero+1][1]) |
---|
| 1028 | || ((indZero+1) >= sizeOfNewtonPolygon && (i+1) > newtonPolyg[0][1])) |
---|
[f0befc] | 1029 | { |
---|
| 1030 | if (indZero + 1 != sizeOfNewtonPolygon) |
---|
| 1031 | indZero++; |
---|
| 1032 | else |
---|
| 1033 | indZero= 0; |
---|
| 1034 | if (indZero != sizeOfNewtonPolygon - 1) |
---|
| 1035 | { |
---|
[9fbf02] | 1036 | slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0]; |
---|
| 1037 | slopeDen= newtonPolyg[indZero+1][1]-newtonPolyg[indZero][1]; |
---|
| 1038 | constTerm= newtonPolyg[indZero][0]; |
---|
[f0befc] | 1039 | } |
---|
| 1040 | else |
---|
| 1041 | { |
---|
[9fbf02] | 1042 | slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0]; |
---|
| 1043 | slopeDen= newtonPolyg[0][1]-newtonPolyg[indZero][1]; |
---|
| 1044 | constTerm= newtonPolyg[indZero][0]; |
---|
[f0befc] | 1045 | } |
---|
| 1046 | if (slopeNum < 0) |
---|
| 1047 | { |
---|
| 1048 | negativeSlope= true; |
---|
| 1049 | slopeNum= - slopeNum; |
---|
[9fbf02] | 1050 | k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/ |
---|
[f0befc] | 1051 | slopeDen) + constTerm; |
---|
| 1052 | } |
---|
| 1053 | else |
---|
[9fbf02] | 1054 | k= (int) (((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen) |
---|
[f0befc] | 1055 | + constTerm; |
---|
| 1056 | } |
---|
| 1057 | else |
---|
| 1058 | { |
---|
| 1059 | if (negativeSlope) |
---|
[9fbf02] | 1060 | k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/ |
---|
[f0befc] | 1061 | slopeDen) + constTerm; |
---|
| 1062 | else |
---|
[9fbf02] | 1063 | k= (int) ((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen |
---|
[f0befc] | 1064 | + constTerm; |
---|
| 1065 | } |
---|
[9fbf02] | 1066 | if (i + 1 > maxY || i + 1 < minY) |
---|
[d8a7da] | 1067 | { |
---|
| 1068 | result [i]= 0; |
---|
| 1069 | continue; |
---|
| 1070 | } |
---|
[9fbf02] | 1071 | |
---|
| 1072 | point [0]= k; |
---|
| 1073 | point [1]= i + 1; |
---|
[f0befc] | 1074 | if (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0) |
---|
| 1075 | k= 0; |
---|
[d8a7da] | 1076 | result [i]= k; |
---|
| 1077 | } |
---|
| 1078 | |
---|
[9fbf02] | 1079 | delete [] point; |
---|
| 1080 | |
---|
[f0befc] | 1081 | for (int i= 0; i < sizeOfNewtonPolygon; i++) |
---|
| 1082 | delete [] newtonPolyg[i]; |
---|
| 1083 | delete [] newtonPolyg; |
---|
| 1084 | |
---|
[d8a7da] | 1085 | return result; |
---|
| 1086 | } |
---|
| 1087 | |
---|
[dd8047] | 1088 | int |
---|
| 1089 | substituteCheck (const CanonicalForm& F, const Variable& x) |
---|
| 1090 | { |
---|
| 1091 | if (F.inCoeffDomain()) |
---|
| 1092 | return 0; |
---|
| 1093 | if (degree (F, x) < 0) |
---|
| 1094 | return 0; |
---|
| 1095 | CanonicalForm f= swapvar (F, F.mvar(), x); |
---|
| 1096 | int sizef= 0; |
---|
| 1097 | for (CFIterator i= f; i.hasTerms(); i++, sizef++) |
---|
| 1098 | { |
---|
| 1099 | if (i.exp() == 1) |
---|
| 1100 | return 0; |
---|
| 1101 | } |
---|
| 1102 | int * expf= new int [sizef]; |
---|
| 1103 | int j= 0; |
---|
| 1104 | for (CFIterator i= f; i.hasTerms(); i++, j++) |
---|
| 1105 | expf [j]= i.exp(); |
---|
| 1106 | |
---|
| 1107 | int indf= sizef - 1; |
---|
| 1108 | if (expf[indf] == 0) |
---|
| 1109 | indf--; |
---|
| 1110 | |
---|
| 1111 | int result= expf[indf]; |
---|
| 1112 | for (int i= indf - 1; i >= 0; i--) |
---|
| 1113 | { |
---|
| 1114 | if (expf [i]%result != 0) |
---|
| 1115 | { |
---|
| 1116 | delete [] expf; |
---|
| 1117 | return 0; |
---|
| 1118 | } |
---|
| 1119 | } |
---|
| 1120 | |
---|
| 1121 | delete [] expf; |
---|
| 1122 | return result; |
---|
| 1123 | } |
---|
| 1124 | |
---|
[686ce3] | 1125 | static int |
---|
| 1126 | substituteCheck (const CanonicalForm& F, const CanonicalForm& G) |
---|
| 1127 | { |
---|
| 1128 | if (F.inCoeffDomain() || G.inCoeffDomain()) |
---|
| 1129 | return 0; |
---|
| 1130 | Variable x= Variable (1); |
---|
| 1131 | if (degree (F, x) <= 1 || degree (G, x) <= 1) |
---|
| 1132 | return 0; |
---|
| 1133 | CanonicalForm f= swapvar (F, F.mvar(), x); |
---|
| 1134 | CanonicalForm g= swapvar (G, G.mvar(), x); |
---|
| 1135 | int sizef= 0; |
---|
| 1136 | int sizeg= 0; |
---|
| 1137 | for (CFIterator i= f; i.hasTerms(); i++, sizef++) |
---|
| 1138 | { |
---|
| 1139 | if (i.exp() == 1) |
---|
| 1140 | return 0; |
---|
| 1141 | } |
---|
| 1142 | for (CFIterator i= g; i.hasTerms(); i++, sizeg++) |
---|
| 1143 | { |
---|
| 1144 | if (i.exp() == 1) |
---|
| 1145 | return 0; |
---|
| 1146 | } |
---|
| 1147 | int * expf= new int [sizef]; |
---|
| 1148 | int * expg= new int [sizeg]; |
---|
| 1149 | int j= 0; |
---|
| 1150 | for (CFIterator i= f; i.hasTerms(); i++, j++) |
---|
| 1151 | { |
---|
| 1152 | expf [j]= i.exp(); |
---|
| 1153 | } |
---|
| 1154 | j= 0; |
---|
| 1155 | for (CFIterator i= g; i.hasTerms(); i++, j++) |
---|
| 1156 | { |
---|
| 1157 | expg [j]= i.exp(); |
---|
| 1158 | } |
---|
| 1159 | |
---|
| 1160 | int indf= sizef - 1; |
---|
| 1161 | int indg= sizeg - 1; |
---|
| 1162 | if (expf[indf] == 0) |
---|
| 1163 | indf--; |
---|
| 1164 | if (expg[indg] == 0) |
---|
| 1165 | indg--; |
---|
| 1166 | |
---|
| 1167 | if ((expg[indg]%expf [indf] != 0 && expf[indf]%expg[indg] != 0) || |
---|
| 1168 | (expg[indg] == 1 && expf[indf] == 1)) |
---|
| 1169 | { |
---|
| 1170 | delete [] expg; |
---|
| 1171 | delete [] expf; |
---|
| 1172 | return 0; |
---|
| 1173 | } |
---|
| 1174 | |
---|
| 1175 | int result; |
---|
| 1176 | if (expg [indg]%expf [indf] == 0) |
---|
| 1177 | result= expf[indf]; |
---|
| 1178 | else |
---|
| 1179 | result= expg[indg]; |
---|
| 1180 | for (int i= indf - 1; i >= 0; i--) |
---|
| 1181 | { |
---|
| 1182 | if (expf [i]%result != 0) |
---|
| 1183 | { |
---|
| 1184 | delete [] expf; |
---|
| 1185 | delete [] expg; |
---|
| 1186 | return 0; |
---|
| 1187 | } |
---|
| 1188 | } |
---|
| 1189 | |
---|
| 1190 | for (int i= indg - 1; i >= 0; i--) |
---|
| 1191 | { |
---|
| 1192 | if (expg [i]%result != 0) |
---|
| 1193 | { |
---|
| 1194 | delete [] expf; |
---|
| 1195 | delete [] expg; |
---|
| 1196 | return 0; |
---|
| 1197 | } |
---|
| 1198 | } |
---|
| 1199 | |
---|
| 1200 | delete [] expg; |
---|
| 1201 | delete [] expf; |
---|
| 1202 | return result; |
---|
| 1203 | } |
---|
| 1204 | |
---|
| 1205 | int recSubstituteCheck (const CanonicalForm& F, const int d) |
---|
| 1206 | { |
---|
| 1207 | if (F.inCoeffDomain()) |
---|
| 1208 | return 0; |
---|
| 1209 | Variable x= Variable (1); |
---|
| 1210 | if (degree (F, x) <= 1) |
---|
| 1211 | return 0; |
---|
| 1212 | CanonicalForm f= swapvar (F, F.mvar(), x); |
---|
| 1213 | int sizef= 0; |
---|
| 1214 | for (CFIterator i= f; i.hasTerms(); i++, sizef++) |
---|
| 1215 | { |
---|
| 1216 | if (i.exp() == 1) |
---|
| 1217 | return 0; |
---|
| 1218 | } |
---|
| 1219 | int * expf= new int [sizef]; |
---|
| 1220 | int j= 0; |
---|
| 1221 | for (CFIterator i= f; i.hasTerms(); i++, j++) |
---|
| 1222 | { |
---|
| 1223 | expf [j]= i.exp(); |
---|
| 1224 | } |
---|
| 1225 | |
---|
| 1226 | int indf= sizef - 1; |
---|
| 1227 | if (expf[indf] == 0) |
---|
| 1228 | indf--; |
---|
| 1229 | |
---|
| 1230 | if ((d%expf [indf] != 0 && expf[indf]%d != 0) || (expf[indf] == 1)) |
---|
| 1231 | { |
---|
| 1232 | delete [] expf; |
---|
| 1233 | return 0; |
---|
| 1234 | } |
---|
| 1235 | |
---|
| 1236 | int result; |
---|
| 1237 | if (d%expf [indf] == 0) |
---|
| 1238 | result= expf[indf]; |
---|
| 1239 | else |
---|
| 1240 | result= d; |
---|
| 1241 | for (int i= indf - 1; i >= 0; i--) |
---|
| 1242 | { |
---|
| 1243 | if (expf [i]%result != 0) |
---|
| 1244 | { |
---|
| 1245 | delete [] expf; |
---|
| 1246 | return 0; |
---|
| 1247 | } |
---|
| 1248 | } |
---|
| 1249 | |
---|
| 1250 | delete [] expf; |
---|
| 1251 | return result; |
---|
| 1252 | } |
---|
| 1253 | |
---|
| 1254 | int substituteCheck (const CFList& L) |
---|
| 1255 | { |
---|
| 1256 | ASSERT (L.length() > 1, "expected a list of at least two elements"); |
---|
| 1257 | if (L.length() < 2) |
---|
| 1258 | return 0; |
---|
| 1259 | CFListIterator i= L; |
---|
| 1260 | i++; |
---|
| 1261 | int result= substituteCheck (L.getFirst(), i.getItem()); |
---|
| 1262 | if (result <= 1) |
---|
| 1263 | return result; |
---|
| 1264 | i++; |
---|
| 1265 | for (;i.hasItem(); i++) |
---|
| 1266 | { |
---|
| 1267 | result= recSubstituteCheck (i.getItem(), result); |
---|
| 1268 | if (result <= 1) |
---|
| 1269 | return result; |
---|
| 1270 | } |
---|
| 1271 | return result; |
---|
| 1272 | } |
---|
| 1273 | |
---|
| 1274 | void |
---|
| 1275 | subst (const CanonicalForm& F, CanonicalForm& A, const int d, const Variable& x) |
---|
| 1276 | { |
---|
| 1277 | if (d <= 1) |
---|
| 1278 | { |
---|
| 1279 | A= F; |
---|
| 1280 | return; |
---|
| 1281 | } |
---|
| 1282 | if (degree (F, x) <= 0) |
---|
| 1283 | { |
---|
| 1284 | A= F; |
---|
| 1285 | return; |
---|
| 1286 | } |
---|
| 1287 | CanonicalForm C= 0; |
---|
| 1288 | CanonicalForm f= swapvar (F, x, F.mvar()); |
---|
| 1289 | for (CFIterator i= f; i.hasTerms(); i++) |
---|
| 1290 | C += i.coeff()*power (f.mvar(), i.exp()/ d); |
---|
| 1291 | A= swapvar (C, x, F.mvar()); |
---|
| 1292 | } |
---|
| 1293 | |
---|
| 1294 | CanonicalForm |
---|
| 1295 | reverseSubst (const CanonicalForm& F, const int d, const Variable& x) |
---|
| 1296 | { |
---|
| 1297 | if (d <= 1) |
---|
| 1298 | return F; |
---|
| 1299 | if (degree (F, x) <= 0) |
---|
| 1300 | return F; |
---|
| 1301 | CanonicalForm f= swapvar (F, x, F.mvar()); |
---|
| 1302 | CanonicalForm result= 0; |
---|
| 1303 | for (CFIterator i= f; i.hasTerms(); i++) |
---|
| 1304 | result += i.coeff()*power (f.mvar(), d*i.exp()); |
---|
| 1305 | return swapvar (result, x, F.mvar()); |
---|
| 1306 | } |
---|
| 1307 | |
---|
| 1308 | void |
---|
| 1309 | reverseSubst (CFList& L, const int d, const Variable& x) |
---|
| 1310 | { |
---|
| 1311 | for (CFListIterator i= L; i.hasItem(); i++) |
---|
| 1312 | i.getItem()= reverseSubst (i.getItem(), d, x); |
---|
| 1313 | } |
---|
| 1314 | |
---|