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