[a7fbdd] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
| 4 | /*************************************************************** |
---|
| 5 | * File: ncSAFormula.cc |
---|
| 6 | * Purpose: implementation of multiplication by formulas in simple NC subalgebras |
---|
| 7 | * Author: motsak |
---|
| 8 | * Created: |
---|
| 9 | *******************************************************************/ |
---|
| 10 | |
---|
| 11 | #define MYTEST 0 |
---|
| 12 | |
---|
| 13 | #if MYTEST |
---|
| 14 | #define OM_CHECK 4 |
---|
| 15 | #define OM_TRACK 5 |
---|
| 16 | // these settings must be before "mod2.h" in order to work!!! |
---|
| 17 | #endif |
---|
| 18 | |
---|
| 19 | |
---|
[d6a97c3] | 20 | #include "config.h" |
---|
| 21 | #include <misc/auxiliary.h> |
---|
[a7fbdd] | 22 | |
---|
[6e05dc] | 23 | #ifdef HAVE_PLURAL |
---|
| 24 | |
---|
[d6a97c3] | 25 | #define PLURAL_INTERNAL_DECLARATIONS |
---|
| 26 | |
---|
[a7fbdd] | 27 | #ifndef NDEBUG |
---|
[1f5565d] | 28 | #define OUTPUT MYTEST |
---|
[a7fbdd] | 29 | #else |
---|
| 30 | #define OUTPUT 0 |
---|
| 31 | #endif |
---|
| 32 | |
---|
[d6a97c3] | 33 | #include <reporter/reporter.h> |
---|
| 34 | |
---|
| 35 | #include <coeffs/numbers.h> |
---|
| 36 | #include "coeffrings.h" |
---|
| 37 | |
---|
| 38 | #include "nc/ncSAFormula.h" |
---|
| 39 | // for CFormulaPowerMultiplier |
---|
[a7fbdd] | 40 | |
---|
[d6a97c3] | 41 | #include "monomials/ring.h" |
---|
| 42 | #include "monomials/p_polys.h" |
---|
[a7fbdd] | 43 | |
---|
[d6a97c3] | 44 | #include "nc/sca.h" |
---|
[a7fbdd] | 45 | |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | bool ncInitSpecialPowersMultiplication(ring r) |
---|
| 50 | { |
---|
| 51 | #if OUTPUT |
---|
| 52 | Print("ncInitSpecialPowersMultiplication(ring), ring: \n"); |
---|
[1f5565d] | 53 | rWrite(r, TRUE); |
---|
[a7fbdd] | 54 | PrintLn(); |
---|
| 55 | #endif |
---|
| 56 | |
---|
| 57 | assume(rIsPluralRing(r)); |
---|
| 58 | assume(!rIsSCA(r)); |
---|
| 59 | |
---|
[b902246] | 60 | |
---|
| 61 | if( r->GetNC()->GetFormulaPowerMultiplier() != NULL ) |
---|
| 62 | { |
---|
| 63 | WarnS("Already defined!"); |
---|
| 64 | return false; |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | |
---|
[a7fbdd] | 68 | r->GetNC()->GetFormulaPowerMultiplier() = new CFormulaPowerMultiplier(r); |
---|
| 69 | |
---|
| 70 | return true; |
---|
| 71 | |
---|
[a60e0b] | 72 | } |
---|
[a7fbdd] | 73 | |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | |
---|
[1f5565d] | 78 | |
---|
| 79 | |
---|
| 80 | // TODO: return q-coeff? |
---|
| 81 | static inline BOOLEAN AreCommutingVariables(const ring r, int i, int j/*, number *qq*/) |
---|
| 82 | { |
---|
| 83 | #if OUTPUT |
---|
| 84 | Print("AreCommutingVariables(ring, k: %d, i: %d)!", j, i); |
---|
| 85 | PrintLn(); |
---|
| 86 | #endif |
---|
| 87 | |
---|
| 88 | assume(i != j); |
---|
| 89 | |
---|
| 90 | assume(i > 0); |
---|
| 91 | assume(i <= r->N); |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | assume(j > 0); |
---|
| 95 | assume(j <= r->N); |
---|
| 96 | |
---|
| 97 | const BOOLEAN reverse = (i > j); |
---|
| 98 | |
---|
| 99 | if (reverse) { int k = j; j = i; i = k; } |
---|
| 100 | |
---|
| 101 | assume(i < j); |
---|
| 102 | |
---|
| 103 | { |
---|
| 104 | const poly d = GetD(r, i, j); |
---|
| 105 | |
---|
| 106 | #if OUTPUT |
---|
| 107 | Print("D_{%d, %d} = ", i, j); p_Write(d, r); |
---|
| 108 | #endif |
---|
| 109 | |
---|
| 110 | if( d != NULL) |
---|
| 111 | return FALSE; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | { |
---|
| 116 | const number q = p_GetCoeff(GetC(r, i, j), r); |
---|
| 117 | |
---|
| 118 | if( !n_IsOne(q, r) ) |
---|
| 119 | return FALSE; |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | return TRUE; // [VAR(I), VAR(J)] = 0!! |
---|
| 123 | |
---|
| 124 | /* |
---|
| 125 | if (reverse) |
---|
| 126 | *qq = n_Invers(q, r); |
---|
| 127 | else |
---|
| 128 | *qq = n_Copy(q, r); |
---|
| 129 | return TRUE; |
---|
| 130 | */ |
---|
| 131 | } |
---|
| 132 | |
---|
[a7fbdd] | 133 | static inline Enum_ncSAType AnalyzePairType(const ring r, int i, int j) |
---|
| 134 | { |
---|
| 135 | #if OUTPUT |
---|
[1f5565d] | 136 | Print("AnalyzePair(ring, i: %d, j: %d):", i, j); |
---|
[a7fbdd] | 137 | PrintLn(); |
---|
| 138 | #endif |
---|
| 139 | |
---|
| 140 | const int N = r->N; |
---|
| 141 | |
---|
| 142 | assume(i < j); |
---|
| 143 | assume(i > 0); |
---|
| 144 | assume(j <= N); |
---|
| 145 | |
---|
| 146 | |
---|
[1f5565d] | 147 | const poly c = GetC(r, i, j); |
---|
| 148 | const number q = p_GetCoeff(c, r); |
---|
[a7fbdd] | 149 | const poly d = GetD(r, i, j); |
---|
| 150 | |
---|
[1f5565d] | 151 | #if 0 && OUTPUT |
---|
| 152 | Print("C_{%d, %d} = ", i, j); p_Write(c, r); PrintLn(); |
---|
[a7fbdd] | 153 | Print("D_{%d, %d} = ", i, j); p_Write(d, r); |
---|
| 154 | #endif |
---|
| 155 | |
---|
| 156 | // const number q = p_GetCoeff(c, r); |
---|
| 157 | |
---|
| 158 | if( d == NULL) |
---|
| 159 | { |
---|
| 160 | |
---|
| 161 | if( n_IsOne(q, r) ) // commutative |
---|
| 162 | return _ncSA_1xy0x0y0; |
---|
| 163 | |
---|
[1f5565d] | 164 | if( n_IsMOne(q, r) ) // anti-commutative |
---|
[a7fbdd] | 165 | return _ncSA_Mxy0x0y0; |
---|
| 166 | |
---|
[1f5565d] | 167 | return _ncSA_Qxy0x0y0; // quasi-commutative |
---|
[a7fbdd] | 168 | } else |
---|
| 169 | { |
---|
| 170 | if( n_IsOne(q, r) ) // "Lie" case |
---|
| 171 | { |
---|
[1f5565d] | 172 | if( pNext(d) == NULL ) // Our Main Special Case: d is only a term! |
---|
[a7fbdd] | 173 | { |
---|
[d2f720] | 174 | // const number g = p_GetCoeff(d, r); // not used for now |
---|
[1f5565d] | 175 | if( p_LmIsConstantComp(d, r) ) // Weyl |
---|
[a7fbdd] | 176 | return _ncSA_1xy0x0yG; |
---|
| 177 | |
---|
| 178 | const int k = p_IsPurePower(d, r); // k if not pure power |
---|
| 179 | |
---|
[1f5565d] | 180 | if( k > 0 ) // d = var(k)^?? |
---|
| 181 | { |
---|
| 182 | const int exp = p_GetExp(d, k, r); |
---|
| 183 | |
---|
| 184 | if (exp == 1) |
---|
[a7fbdd] | 185 | { |
---|
[1f5565d] | 186 | if(k == i) // 2 -ubalgebra in var(i) & var(j), with linear relation...? |
---|
[a7fbdd] | 187 | return _ncSA_1xyAx0y0; |
---|
| 188 | |
---|
| 189 | if(k == j) |
---|
| 190 | return _ncSA_1xy0xBy0; |
---|
[1f5565d] | 191 | } else if ( exp == 2 && k!= i && k != j) // Homogenized Weyl algebra [x, Dx] = t^2? |
---|
| 192 | { |
---|
| 193 | // number qi, qj; |
---|
| 194 | if (AreCommutingVariables(r, k, i/*, &qi*/) && AreCommutingVariables(r, k, j/*, &qj*/) ) // [x, t] = [Dx, t] = 0? |
---|
| 195 | { |
---|
| 196 | const number g = p_GetCoeff(d, r); |
---|
| 197 | |
---|
| 198 | if (n_IsOne(g, r)) |
---|
| 199 | return _ncSA_1xy0x0yT2; // save k!?, G = LC(d) == qi == qj == 1!!! |
---|
| 200 | } |
---|
[a7fbdd] | 201 | } |
---|
[1f5565d] | 202 | } |
---|
[a7fbdd] | 203 | } |
---|
| 204 | } |
---|
[1f5565d] | 205 | // Hmm, what about a more general case of q != 1??? |
---|
[a7fbdd] | 206 | } |
---|
[1f5565d] | 207 | #if OUTPUT |
---|
| 208 | Print("C_{%d, %d} = ", i, j); p_Write(c, r); |
---|
| 209 | Print("D_{%d, %d} = ", i, j); p_Write(d, r); |
---|
| 210 | PrintS("====>>>>_ncSA_notImplemented\n"); |
---|
| 211 | #endif |
---|
[a7fbdd] | 212 | |
---|
| 213 | return _ncSA_notImplemented; |
---|
[a60e0b] | 214 | } |
---|
[a7fbdd] | 215 | |
---|
| 216 | |
---|
[cc4cc80] | 217 | CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring r): m_NVars(r->N), m_BaseRing(r) |
---|
[a7fbdd] | 218 | { |
---|
| 219 | #if OUTPUT |
---|
| 220 | Print("CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring)!"); |
---|
| 221 | PrintLn(); |
---|
| 222 | #endif |
---|
| 223 | |
---|
| 224 | m_SAPairTypes = (Enum_ncSAType*)omAlloc0( ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) ); |
---|
| 225 | |
---|
| 226 | for( int i = 1; i < NVars(); i++ ) |
---|
| 227 | for( int j = i + 1; j <= NVars(); j++ ) |
---|
| 228 | GetPair(i, j) = AnalyzePairType(GetBasering(), i, j); |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | |
---|
| 233 | |
---|
| 234 | CFormulaPowerMultiplier::~CFormulaPowerMultiplier() |
---|
| 235 | { |
---|
| 236 | #if OUTPUT |
---|
| 237 | Print("CFormulaPowerMultiplier::~CFormulaPowerMultiplier()!"); |
---|
| 238 | PrintLn(); |
---|
| 239 | #endif |
---|
| 240 | |
---|
| 241 | omFreeSize((ADDRESS)m_SAPairTypes, ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) ); |
---|
| 242 | } |
---|
| 243 | |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | |
---|
| 247 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 248 | static inline void CorrectPolyWRTOrdering(poly &pResult, const ring r) |
---|
| 249 | { |
---|
| 250 | if( pNext(pResult) != NULL ) |
---|
| 251 | { |
---|
| 252 | const int cmp = p_LmCmp(pResult, pNext(pResult), r); |
---|
| 253 | assume( cmp != 0 ); // must not be equal!!! |
---|
| 254 | if( cmp != 1 ) // Wrong order!!! |
---|
| 255 | pResult = pReverse(pResult); // Reverse!!! |
---|
| 256 | } |
---|
| 257 | p_Test(pResult, r); |
---|
| 258 | } |
---|
| 259 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 260 | static inline poly ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r) |
---|
| 261 | { |
---|
| 262 | #if OUTPUT |
---|
| 263 | Print("ncSA_1xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); |
---|
| 264 | PrintLn(); |
---|
| 265 | #endif |
---|
| 266 | |
---|
[b902246] | 267 | poly p = p_One( r); |
---|
[a7fbdd] | 268 | p_SetExp(p, j, m, r); |
---|
| 269 | p_SetExp(p, i, n, r); |
---|
| 270 | p_Setm(p, r); |
---|
| 271 | |
---|
| 272 | p_Test(p, r); |
---|
| 273 | |
---|
| 274 | return p; |
---|
| 275 | |
---|
| 276 | } // return ncSA_1xy0x0y0(GetI(), GetJ(), expRight, expLeft, r); |
---|
| 277 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 278 | static inline poly ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r) |
---|
| 279 | { |
---|
| 280 | #if OUTPUT |
---|
| 281 | Print("ncSA_{M = -1}xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); |
---|
| 282 | PrintLn(); |
---|
| 283 | #endif |
---|
| 284 | |
---|
| 285 | const int sign = 1 - ((n & (m & 1)) << 1); |
---|
| 286 | poly p = p_ISet(sign, r); |
---|
| 287 | p_SetExp(p, j, m, r); |
---|
| 288 | p_SetExp(p, i, n, r); |
---|
| 289 | p_Setm(p, r); |
---|
| 290 | |
---|
| 291 | |
---|
| 292 | p_Test(p, r); |
---|
| 293 | |
---|
| 294 | return p; |
---|
| 295 | |
---|
| 296 | } // return ncSA_Mxy0x0y0(GetI(), GetJ(), expRight, expLeft, r); |
---|
| 297 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 298 | static inline poly ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r) |
---|
| 299 | { |
---|
| 300 | #if OUTPUT |
---|
| 301 | Print("ncSA_Qxy0x0y0(var(%d)^{%d}, var(%d)^{%d}, Q, r)!", j, m, i, n); |
---|
| 302 | PrintLn(); |
---|
| 303 | #endif |
---|
| 304 | |
---|
| 305 | int min, max; |
---|
| 306 | |
---|
| 307 | if( n < m ) |
---|
| 308 | { |
---|
| 309 | min = n; |
---|
| 310 | max = m; |
---|
| 311 | } else |
---|
| 312 | { |
---|
| 313 | min = m; |
---|
| 314 | max = n; |
---|
| 315 | } |
---|
| 316 | |
---|
| 317 | number qN; |
---|
| 318 | |
---|
| 319 | if( max == 1 ) |
---|
| 320 | qN = n_Copy(m_q, r); |
---|
| 321 | else |
---|
| 322 | { |
---|
| 323 | number t; |
---|
| 324 | n_Power(m_q, max, &t, r); |
---|
| 325 | |
---|
| 326 | if( min > 1 ) |
---|
| 327 | { |
---|
| 328 | n_Power(t, min, &qN, r); |
---|
| 329 | n_Delete(&t, r); |
---|
| 330 | } |
---|
| 331 | else |
---|
| 332 | qN = t; |
---|
| 333 | } |
---|
| 334 | |
---|
| 335 | poly p = p_NSet(qN, r); |
---|
| 336 | p_SetExp(p, j, m, r); |
---|
| 337 | p_SetExp(p, i, n, r); |
---|
| 338 | p_Setm(p, r); |
---|
| 339 | |
---|
| 340 | |
---|
| 341 | p_Test(p, r); |
---|
| 342 | |
---|
| 343 | return p; |
---|
| 344 | |
---|
| 345 | } // return ncSA_Qxy0x0y0(GetI(), GetJ(), expRight, expLeft, m_q, r); |
---|
| 346 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 347 | static inline poly ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r) |
---|
| 348 | { |
---|
| 349 | #if OUTPUT |
---|
| 350 | Print("ncSA_1xy0x0yG(var(%d)^{%d}, var(%d)^{%d}, G, r)!", j, m, i, n); |
---|
| 351 | PrintLn(); |
---|
| 352 | number t = n_Copy(m_g, r); |
---|
| 353 | PrintS("Parameter G: "); n_Write(t, r); |
---|
| 354 | n_Delete(&t, r); |
---|
| 355 | #endif |
---|
| 356 | |
---|
| 357 | int kn = n; |
---|
| 358 | int km = m; |
---|
| 359 | |
---|
| 360 | number c = n_Init(1, r); |
---|
| 361 | |
---|
[b902246] | 362 | poly p = p_One( r); |
---|
[a7fbdd] | 363 | |
---|
| 364 | p_SetExp(p, j, km--, r); // y ^ (m-k) |
---|
| 365 | p_SetExp(p, i, kn--, r); // x ^ (n-k) |
---|
| 366 | |
---|
| 367 | p_Setm(p, r); // pResult = x^n * y^m |
---|
| 368 | |
---|
| 369 | |
---|
| 370 | poly pResult = p; |
---|
| 371 | poly pLast = p; |
---|
| 372 | |
---|
| 373 | int min = si_min(m, n); |
---|
| 374 | |
---|
| 375 | int k = 1; |
---|
| 376 | |
---|
| 377 | for(; k < min; k++ ) |
---|
| 378 | { |
---|
| 379 | number t = n_Init(km + 1, r); |
---|
| 380 | n_InpMult(t, m_g, r); // t = ((m - k) + 1) * gamma |
---|
| 381 | n_InpMult(c, t, r); // c = c'* ((m - k) + 1) * gamma |
---|
| 382 | n_Delete(&t, r); |
---|
| 383 | |
---|
| 384 | t = n_Init(kn + 1, r); |
---|
| 385 | n_InpMult(c, t, r); // c = (c'* ((m - k) + 1) * gamma) * ((n - k) + 1) |
---|
| 386 | n_Delete(&t, r); |
---|
| 387 | |
---|
| 388 | t = n_Init(k, r); |
---|
| 389 | c = n_Div(c, t, r); |
---|
| 390 | n_Delete(&t, r); |
---|
| 391 | |
---|
| 392 | // n_Normalize(c, r); |
---|
| 393 | |
---|
| 394 | t = n_Copy(c, r); // not the last! |
---|
| 395 | |
---|
| 396 | p = p_NSet(t, r); |
---|
| 397 | |
---|
| 398 | p_SetExp(p, j, km--, r); // y ^ (m-k) |
---|
| 399 | p_SetExp(p, i, kn--, r); // x ^ (n-k) |
---|
| 400 | |
---|
| 401 | p_Setm(p, r); // pResult = x^n * y^m |
---|
| 402 | |
---|
| 403 | pNext(pLast) = p; |
---|
| 404 | pLast = p; |
---|
| 405 | } |
---|
| 406 | |
---|
| 407 | assume(k == min); |
---|
| 408 | assume((km == 0) || (kn == 0) ); |
---|
| 409 | |
---|
| 410 | { |
---|
| 411 | n_InpMult(c, m_g, r); // c = c'* gamma |
---|
| 412 | |
---|
| 413 | if( km > 0 ) |
---|
| 414 | { |
---|
| 415 | number t = n_Init(km + 1, r); |
---|
| 416 | n_InpMult(c, t, r); // c = (c'* gamma) * (m - k + 1) |
---|
| 417 | n_Delete(&t, r); |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | if( kn > 0 ) |
---|
| 421 | { |
---|
| 422 | number t = n_Init(kn + 1, r); |
---|
| 423 | n_InpMult(c, t, r); // c = (c'* gamma) * (n - k + 1) |
---|
| 424 | n_Delete(&t, r); |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | number t = n_Init(k, r); // c = ((c'* gamma) * ((n - k + 1) * (m - k + 1))) / k; |
---|
| 428 | c = n_Div(c, t, r); |
---|
| 429 | n_Delete(&t, r); |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | p = p_NSet(c, r); |
---|
| 433 | |
---|
| 434 | p_SetExp(p, j, km, r); // y ^ (m-k) |
---|
| 435 | p_SetExp(p, i, kn, r); // x ^ (n-k) |
---|
| 436 | |
---|
| 437 | p_Setm(p, r); // |
---|
| 438 | |
---|
| 439 | pNext(pLast) = p; |
---|
| 440 | |
---|
| 441 | CorrectPolyWRTOrdering(pResult, r); |
---|
| 442 | |
---|
| 443 | return pResult; |
---|
| 444 | } |
---|
| 445 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
[1f5565d] | 446 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 447 | static inline poly ncSA_1xy0x0yT2(const int i, const int j, const int n, const int m, const int m_k, const ring r) |
---|
| 448 | { |
---|
| 449 | #if OUTPUT |
---|
| 450 | Print("ncSA_1xy0x0yT2(var(%d)^{%d}, var(%d)^{%d}, t: var(%d), r)!", j, m, i, n, m_k); |
---|
| 451 | PrintLn(); |
---|
| 452 | #endif |
---|
| 453 | |
---|
| 454 | int kn = n; |
---|
| 455 | int km = m; |
---|
| 456 | |
---|
| 457 | // k == 0! |
---|
| 458 | number c = n_Init(1, r); |
---|
| 459 | |
---|
| 460 | poly p = p_One( r ); |
---|
| 461 | |
---|
| 462 | p_SetExp(p, j, km--, r); // y ^ (m) |
---|
| 463 | p_SetExp(p, i, kn--, r); // x ^ (n) |
---|
| 464 | // p_SetExp(p, m_k, k << 1, r); // homogenization with var(m_k) ^ (2*k) |
---|
| 465 | |
---|
| 466 | p_Setm(p, r); // pResult = x^n * y^m |
---|
| 467 | |
---|
| 468 | |
---|
| 469 | poly pResult = p; |
---|
| 470 | poly pLast = p; |
---|
| 471 | |
---|
| 472 | int min = si_min(m, n); |
---|
| 473 | |
---|
| 474 | int k = 1; |
---|
| 475 | |
---|
| 476 | for(; k < min; k++ ) |
---|
| 477 | { |
---|
| 478 | number t = n_Init(km + 1, r); |
---|
| 479 | // n_InpMult(t, m_g, r); // t = ((m - k) + 1) * gamma |
---|
| 480 | n_InpMult(c, t, r); // c = c'* ((m - k) + 1) * gamma |
---|
| 481 | n_Delete(&t, r); |
---|
| 482 | |
---|
| 483 | t = n_Init(kn + 1, r); |
---|
| 484 | n_InpMult(c, t, r); // c = (c'* ((m - k) + 1) * gamma) * ((n - k) + 1) |
---|
| 485 | n_Delete(&t, r); |
---|
| 486 | |
---|
| 487 | t = n_Init(k, r); |
---|
| 488 | c = n_Div(c, t, r); |
---|
| 489 | n_Delete(&t, r); |
---|
| 490 | |
---|
| 491 | // // n_Normalize(c, r); |
---|
| 492 | |
---|
| 493 | t = n_Copy(c, r); // not the last! |
---|
| 494 | |
---|
| 495 | p = p_NSet(t, r); |
---|
| 496 | |
---|
| 497 | p_SetExp(p, j, km--, r); // y ^ (m-k) |
---|
| 498 | p_SetExp(p, i, kn--, r); // x ^ (n-k) |
---|
| 499 | |
---|
| 500 | p_SetExp(p, m_k, k << 1, r); // homogenization with var(m_k) ^ (2*k) |
---|
| 501 | |
---|
| 502 | p_Setm(p, r); // pResult = x^(n-k) * y^(m-k) |
---|
| 503 | |
---|
| 504 | pNext(pLast) = p; |
---|
| 505 | pLast = p; |
---|
| 506 | } |
---|
| 507 | |
---|
| 508 | assume(k == min); |
---|
| 509 | assume((km == 0) || (kn == 0) ); |
---|
| 510 | |
---|
| 511 | { |
---|
| 512 | // n_InpMult(c, m_g, r); // c = c'* gamma |
---|
| 513 | |
---|
| 514 | if( km > 0 ) |
---|
| 515 | { |
---|
| 516 | number t = n_Init(km + 1, r); |
---|
| 517 | n_InpMult(c, t, r); // c = (c'* gamma) * (m - k + 1) |
---|
| 518 | n_Delete(&t, r); |
---|
| 519 | } |
---|
| 520 | |
---|
| 521 | if( kn > 0 ) |
---|
| 522 | { |
---|
| 523 | number t = n_Init(kn + 1, r); |
---|
| 524 | n_InpMult(c, t, r); // c = (c'* gamma) * (n - k + 1) |
---|
| 525 | n_Delete(&t, r); |
---|
| 526 | } |
---|
| 527 | |
---|
| 528 | number t = n_Init(k, r); // c = ((c'* gamma) * ((n - k + 1) * (m - k + 1))) / k; |
---|
| 529 | c = n_Div(c, t, r); |
---|
| 530 | n_Delete(&t, r); |
---|
| 531 | } |
---|
| 532 | |
---|
| 533 | p = p_NSet(c, r); |
---|
| 534 | |
---|
| 535 | p_SetExp(p, j, km, r); // y ^ (m-k) |
---|
| 536 | p_SetExp(p, i, kn, r); // x ^ (n-k) |
---|
| 537 | |
---|
| 538 | p_SetExp(p, m_k, k << 1, r); // homogenization with var(m_k) ^ (2*k) |
---|
| 539 | |
---|
| 540 | p_Setm(p, r); // |
---|
| 541 | |
---|
| 542 | pNext(pLast) = p; |
---|
| 543 | |
---|
| 544 | CorrectPolyWRTOrdering(pResult, r); |
---|
| 545 | |
---|
| 546 | return pResult; |
---|
| 547 | } |
---|
| 548 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 549 | |
---|
| 550 | |
---|
| 551 | |
---|
[a7fbdd] | 552 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 553 | static inline poly ncSA_ShiftAx(int i, int j, int n, int m, const number m_shiftCoef, const ring r) |
---|
| 554 | { |
---|
| 555 | // Char == 0, otherwise - problem! |
---|
| 556 | |
---|
| 557 | int k = m; // to 0 |
---|
| 558 | |
---|
| 559 | number c = n_Init(1, r); // k = m, C_k = 1 |
---|
[b902246] | 560 | poly p = p_One( r); |
---|
[a7fbdd] | 561 | |
---|
| 562 | p_SetExp(p, j, k, r); // Y^{k} |
---|
| 563 | p_SetExp(p, i, n, r); |
---|
| 564 | |
---|
| 565 | p_Setm(p, r); // pResult = C_k * x^n * y^k, k == m |
---|
| 566 | |
---|
| 567 | |
---|
| 568 | poly pResult = p; |
---|
| 569 | poly pLast = p; |
---|
| 570 | |
---|
| 571 | number nn = n_Init(n, r); // number(n)! |
---|
| 572 | n_InpMult(nn, m_shiftCoef, r); // nn = (alpha*n) |
---|
| 573 | |
---|
| 574 | --k; |
---|
| 575 | |
---|
| 576 | int mk = 1; // mk = (m - k) |
---|
| 577 | |
---|
| 578 | for(; k > 0; k-- ) |
---|
| 579 | { |
---|
| 580 | number t = n_Init(k + 1, r); // t = k+1 |
---|
| 581 | n_InpMult(c, t, r); // c = c' * (k+1) |
---|
| 582 | n_InpMult(c, nn, r); // c = (c' * (k+1)) * (alpha * n) |
---|
| 583 | |
---|
| 584 | n_Delete(&t, r); |
---|
| 585 | t = n_Init(mk++, r); |
---|
| 586 | c = n_Div(c, t, r); // c = ((c' * (k+1)) * (alpha * n)) / (m-k); |
---|
| 587 | n_Delete(&t, r); |
---|
| 588 | |
---|
| 589 | // n_Normalize(c, r); |
---|
| 590 | |
---|
| 591 | t = n_Copy(c, r); // not the last! |
---|
| 592 | |
---|
| 593 | p = p_NSet(t, r); |
---|
| 594 | |
---|
| 595 | p_SetExp(p, j, k, r); // y^k |
---|
| 596 | p_SetExp(p, i, n, r); // x^n |
---|
| 597 | |
---|
| 598 | p_Setm(p, r); // pResult = x^n * y^m |
---|
| 599 | |
---|
| 600 | pNext(pLast) = p; |
---|
| 601 | pLast = p; |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | assume(k == 0); |
---|
| 605 | |
---|
| 606 | { |
---|
| 607 | n_InpMult(c, nn, r); // c = (c' * (0+1)) * (alpha * n) |
---|
| 608 | |
---|
| 609 | number t = n_Init(m, r); |
---|
| 610 | c = n_Div(c, t, r); // c = ((c' * (0+1)) * (alpha * n)) / (m-0); |
---|
| 611 | n_Delete(&t, r); |
---|
| 612 | } |
---|
| 613 | |
---|
| 614 | n_Delete(&nn, r); |
---|
| 615 | |
---|
| 616 | p = p_NSet(c, r); |
---|
| 617 | |
---|
| 618 | p_SetExp(p, j, k, r); // y^k |
---|
| 619 | p_SetExp(p, i, n, r); // x^n |
---|
| 620 | |
---|
| 621 | p_Setm(p, r); // |
---|
| 622 | |
---|
| 623 | pNext(pLast) = p; |
---|
| 624 | |
---|
| 625 | CorrectPolyWRTOrdering(pResult, r); |
---|
| 626 | |
---|
| 627 | return pResult; |
---|
| 628 | } |
---|
| 629 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 630 | static inline poly ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r) |
---|
| 631 | { |
---|
| 632 | #if OUTPUT |
---|
| 633 | Print("ncSA_1xyAx0y0(var(%d)^{%d}, var(%d)^{%d}, A, r)!", j, m, i, n); |
---|
| 634 | PrintLn(); |
---|
| 635 | number t = n_Copy(m_shiftCoef, r); |
---|
| 636 | PrintS("Parameter A: "); n_Write(t, r); |
---|
| 637 | n_Delete(&t, r); |
---|
| 638 | #endif |
---|
| 639 | |
---|
| 640 | return ncSA_ShiftAx(i, j, n, m, m_shiftCoef, r); |
---|
| 641 | } |
---|
| 642 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 643 | static inline poly ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r) |
---|
| 644 | { |
---|
| 645 | #if OUTPUT |
---|
| 646 | Print("ncSA_1xy0xBy0(var(%d)^{%d}, var(%d)^{%d}, B, r)!", j, m, i, n); |
---|
| 647 | PrintLn(); |
---|
| 648 | number t = n_Copy(m_shiftCoef, r); |
---|
| 649 | PrintS("Parameter B: "); n_Write(t, r); |
---|
| 650 | n_Delete(&t, r); |
---|
| 651 | #endif |
---|
| 652 | |
---|
| 653 | return ncSA_ShiftAx(j, i, m, n, m_shiftCoef, r); |
---|
| 654 | } |
---|
| 655 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 656 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 657 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 658 | /////////////////////////////////////////////////////////////////////////////////////////// |
---|
| 659 | |
---|
| 660 | |
---|
| 661 | static inline poly ncSA_Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r) |
---|
| 662 | { |
---|
| 663 | #if OUTPUT |
---|
| 664 | Print("ncSA_Multiply(type: %d, ring, (var(%d)^{%d} * var(%d)^{%d}, r)!", (int)type, j, m, i, n); |
---|
| 665 | PrintLn(); |
---|
| 666 | #endif |
---|
| 667 | |
---|
| 668 | assume( type != _ncSA_notImplemented ); |
---|
| 669 | assume( (n > 0) && (m > 0) ); |
---|
| 670 | |
---|
| 671 | if( type == _ncSA_1xy0x0y0 ) |
---|
| 672 | return ::ncSA_1xy0x0y0(i, j, n, m, r); |
---|
| 673 | |
---|
| 674 | if( type == _ncSA_Mxy0x0y0 ) |
---|
| 675 | return ::ncSA_Mxy0x0y0(i, j, n, m, r); |
---|
| 676 | |
---|
| 677 | if( type == _ncSA_Qxy0x0y0 ) |
---|
| 678 | { |
---|
| 679 | const number q = p_GetCoeff(GetC(r, i, j), r); |
---|
| 680 | return ::ncSA_Qxy0x0y0(i, j, n, m, q, r); |
---|
| 681 | } |
---|
| 682 | |
---|
[1f5565d] | 683 | const poly d = GetD(r, i, j); |
---|
| 684 | const number g = p_GetCoeff(d, r); |
---|
[a7fbdd] | 685 | |
---|
| 686 | if( type == _ncSA_1xy0x0yG ) // Weyl |
---|
| 687 | return ::ncSA_1xy0x0yG(i, j, n, m, g, r); |
---|
| 688 | |
---|
[1f5565d] | 689 | if( type == _ncSA_1xy0x0yT2 ) // Homogenous Weyl... |
---|
| 690 | return ::ncSA_1xy0x0yT2(i, j, n, m, p_IsPurePower(d, r), r); |
---|
| 691 | |
---|
[a7fbdd] | 692 | if( type == _ncSA_1xyAx0y0 ) // Shift 1 |
---|
| 693 | return ::ncSA_1xyAx0y0(i, j, n, m, g, r); |
---|
| 694 | |
---|
| 695 | if( type == _ncSA_1xy0xBy0 ) // Shift 2 |
---|
| 696 | return ::ncSA_1xy0xBy0(i, j, n, m, g, r); |
---|
| 697 | |
---|
| 698 | assume( type == _ncSA_notImplemented ); |
---|
| 699 | |
---|
| 700 | return NULL; |
---|
| 701 | } |
---|
| 702 | |
---|
| 703 | |
---|
| 704 | poly CFormulaPowerMultiplier::Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r) |
---|
| 705 | { |
---|
| 706 | return ncSA_Multiply( type, i, j, n, m, r); |
---|
| 707 | } |
---|
| 708 | |
---|
| 709 | |
---|
| 710 | Enum_ncSAType CFormulaPowerMultiplier::AnalyzePair(const ring r, int i, int j) |
---|
| 711 | { |
---|
| 712 | return ::AnalyzePairType( r, i, j); |
---|
| 713 | } |
---|
| 714 | |
---|
| 715 | poly CFormulaPowerMultiplier::Multiply( int i, int j, const int n, const int m) |
---|
| 716 | { |
---|
| 717 | return ncSA_Multiply( GetPair(i, j), i, j, n, m, GetBasering()); |
---|
| 718 | } |
---|
| 719 | |
---|
| 720 | |
---|
| 721 | |
---|
| 722 | |
---|
| 723 | poly CFormulaPowerMultiplier::ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r) |
---|
| 724 | { |
---|
| 725 | return ::ncSA_1xy0x0y0(i, j, n, m, r); |
---|
| 726 | } |
---|
| 727 | |
---|
| 728 | poly CFormulaPowerMultiplier::ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r) |
---|
| 729 | { |
---|
| 730 | return ::ncSA_Mxy0x0y0(i, j, n, m, r); |
---|
| 731 | } |
---|
| 732 | |
---|
| 733 | poly CFormulaPowerMultiplier::ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r) |
---|
| 734 | { |
---|
| 735 | return ::ncSA_Qxy0x0y0(i, j, n, m, m_q, r); |
---|
[a60e0b] | 736 | } |
---|
[a7fbdd] | 737 | |
---|
| 738 | poly CFormulaPowerMultiplier::ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r) |
---|
| 739 | { |
---|
| 740 | return ::ncSA_1xy0x0yG(i, j, n, m, m_g, r); |
---|
| 741 | } |
---|
| 742 | |
---|
[1f5565d] | 743 | poly CFormulaPowerMultiplier::ncSA_1xy0x0yT2(const int i, const int j, const int n, const int m, const int k, const ring r) |
---|
| 744 | { |
---|
| 745 | return ::ncSA_1xy0x0yT2(i, j, n, m, k, r); |
---|
| 746 | } |
---|
| 747 | |
---|
[a7fbdd] | 748 | poly CFormulaPowerMultiplier::ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r) |
---|
| 749 | { |
---|
| 750 | return ::ncSA_1xyAx0y0(i, j, n, m, m_shiftCoef, r); |
---|
| 751 | } |
---|
| 752 | |
---|
| 753 | poly CFormulaPowerMultiplier::ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r) |
---|
| 754 | { |
---|
| 755 | return ::ncSA_1xy0xBy0(i, j, n, m, m_shiftCoef, r); |
---|
| 756 | } |
---|
[6e05dc] | 757 | #endif |
---|