source: git/kernel/ncSAMult.cc @ 2f2bb21

spielwiese
Last change on this file since 2f2bb21 was 63774ec, checked in by Motsak Oleksandr <motsak@…>, 16 years ago
*motsak: zero inout :( git-svn-id: file:///usr/local/Singular/svn/trunk@10935 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 24.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    ncSAMult.cc
6 *  Purpose: implementation of multiplication in simple NC subalgebras
7 *  Author:  motsak
8 *  Created:
9 *  Version: $Id: ncSAMult.cc,v 1.10 2008-07-29 13:28:21 motsak Exp $
10 *******************************************************************/
11
12#define MYTEST 0
13
14#if MYTEST
15#define OM_CHECK 4
16#define OM_TRACK 5
17// these settings must be before "mod2.h" in order to work!!!
18#endif
19
20
21#include "mod2.h"
22
23#ifndef NDEBUG
24#define OUTPUT 1
25#else
26#define OUTPUT 0
27#endif
28
29
30#include <ncSAMult.h> // for CMultiplier etc classes
31#include <sca.h> // for SCA
32#include <polys.h> // for p_One
33
34
35
36// poly functions defined in p_Procs: ;
37static poly gnc_pp_Mult_mm(const poly p, const poly m, const ring r, poly& last)
38{
39  if( (p == NULL) || (m == NULL) )
40    return NULL;
41
42  assume( (p != NULL) && (m != NULL) && (r != NULL) );
43
44#if OUTPUT 
45  Print("VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV gnc_pp_Mult_mm(p, m) VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV ");
46  PrintLn();
47  PrintS("p: "); p_Write(p, r);   
48  PrintS("m: "); p_Write(m, r);     
49#endif
50  poly pResult;
51 
52  if (p_IsConstant(m, r))
53    pResult = pp_Mult_nn(p, p_GetCoeff(m,r),r);
54  else
55  { 
56    CGlobalMultiplier* const pMultiplier = r->GetNC()->GetGlobalMultiplier();
57    assume( pMultiplier != NULL );
58
59    poly pMonom = pMultiplier->LM(m, r);
60    pResult = pMultiplier->MultiplyPE(p, pMonom);
61    p_Delete(&pMonom, r);
62    p_Test(pResult, r);
63    pResult = p_Mult_nn(pResult, p_GetCoeff(m, r), r);
64  }
65
66#if OUTPUT 
67  p_Test(pResult, r);
68
69  Print("gnc_pp_Mult_mm(p, m) => "); p_Write(pResult, r);
70  PrintS("p: "); p_Write(p, r);   
71  PrintS("m: "); p_Write(m, r);     
72  Print("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ");
73  PrintLn();
74#endif
75
76  return pResult;
77
78}
79
80static poly gnc_p_Mult_mm(poly p, const poly m, const ring r)
81{
82  if( (p == NULL) || (m == NULL) )
83  {
84    p_Delete(&p, r);
85    return NULL;
86  }
87
88  assume( (p != NULL) && (m != NULL) && (r != NULL) );
89
90#if OUTPUT 
91  Print("VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV gnc_p_Mult_mm(p, m) VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV ");
92  PrintLn();
93  PrintS("p: ");
94  p_Write(p, r);   
95  PrintS("m: ");
96  p_Write(m, r); 
97#endif
98
99  poly pResult;
100
101  if (p_IsConstant(m, r))
102    pResult = p_Mult_nn(p, p_GetCoeff(m,r),r);
103  else
104  { 
105    CGlobalMultiplier* const pMultiplier = r->GetNC()->GetGlobalMultiplier();
106    assume( pMultiplier != NULL );
107
108    poly pMonom = pMultiplier->LM(m, r);
109    pResult = pMultiplier->MultiplyPEDestroy(p, pMonom);
110    p_Delete(&pMonom, r);
111    p_Test(pResult, r);
112    pResult = p_Mult_nn(pResult, p_GetCoeff(m, r), r);
113  }
114
115#if OUTPUT 
116  p_Test(pResult, r);
117
118  Print("gnc_p_Mult_mm(p, m) => "); p_Write(pResult, r);     
119//  PrintS("p: "); p_Write(p, r);   
120  PrintS("m: "); p_Write(m, r);     
121  Print("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ");
122  PrintLn();
123#endif
124 
125  return pResult;
126
127}
128
129static poly gnc_mm_Mult_p(const poly m, poly p, const ring r)
130{
131  if( (p == NULL) || (m == NULL) )
132  {
133    p_Delete(&p, r);
134    return NULL;
135  }
136
137  assume( (p != NULL) && (m != NULL) && (r != NULL) );
138
139  p_Test(m, r);
140  p_Test(p, r);
141
142#if OUTPUT 
143  Print("VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV gnc_mm_Mult_p(m, p) VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV ");
144  PrintLn();
145  PrintS("m: "); p_Write(m, r);     
146  PrintS("p: "); p_Write(p, r);   
147#endif
148
149  poly pResult;
150
151  if (p_IsConstant(m, r))
152    pResult = p_Mult_nn(p, p_GetCoeff(m,r),r);
153  else
154  { 
155    CGlobalMultiplier* const pMultiplier = r->GetNC()->GetGlobalMultiplier();
156    assume( pMultiplier != NULL );
157
158    poly pMonom = pMultiplier->LM(m, r);
159    pResult = pMultiplier->MultiplyEPDestroy(pMonom, p);
160    p_Delete(&pMonom, r);
161    p_Test(pResult, r);
162    pResult = p_Mult_nn(pResult, p_GetCoeff(m, r), r);
163  }
164 
165#if OUTPUT 
166  p_Test(pResult, r);
167
168  Print("gnc_mm_Mult_p(m, p) => "); p_Write(pResult, r);     
169//  PrintS("p: "); p_Write(p, r);   
170  PrintS("m: "); p_Write(m, r);     
171  Print("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ");
172  PrintLn();
173#endif
174 
175  return pResult;
176}
177
178static poly gnc_mm_Mult_pp(const poly m, const poly p, const ring r)
179{
180  if( (p == NULL) || (m == NULL) )
181  {
182    return NULL;
183  }
184
185  assume( (p != NULL) && (m != NULL) && (r != NULL) );
186
187  p_Test(m, r);
188  p_Test(p, r);
189 
190#if OUTPUT 
191  Print("VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV gnc_mm_Mult_pp(m, p) VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV ");
192  PrintLn();
193  PrintS("m: "); p_Write(m, r);     
194  PrintS("p: "); p_Write(p, r);   
195#endif
196 
197  poly pResult;
198
199  if (p_IsConstant(m, r))
200    pResult = pp_Mult_nn(p, p_GetCoeff(m,r),r);
201  else
202  { 
203    CGlobalMultiplier* const pMultiplier = r->GetNC()->GetGlobalMultiplier();
204    assume( pMultiplier != NULL );
205
206    poly pMonom = pMultiplier->LM(m, r);
207    pResult = pMultiplier->MultiplyEP(pMonom, p);
208    p_Delete(&pMonom, r);
209    p_Test(pResult, r);
210    pResult = p_Mult_nn(pResult, p_GetCoeff(m, r), r);
211  }
212
213#if OUTPUT 
214  p_Test(pResult, r);
215
216  Print("gnc_mm_Mult_pp(m, p) => "); p_Write(pResult, r);     
217  PrintS("p: "); p_Write(p, r);   
218  PrintS("m: "); p_Write(m, r);     
219  Print("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ");
220  PrintLn();
221#endif
222 
223  return pResult;
224}
225
226static void gnc_p_ProcsSet(ring rGR, p_Procs_s* p_Procs = NULL)
227{
228#if OUTPUT 
229  Print("|gnc_p_ProcsSet()");
230  PrintLn();
231#endif
232
233  if( p_Procs == NULL )
234    p_Procs = rGR->p_Procs;
235 
236  // "commutative"
237  p_Procs->p_Mult_mm  = rGR->p_Procs->p_Mult_mm  = gnc_p_Mult_mm;
238  p_Procs->pp_Mult_mm = rGR->p_Procs->pp_Mult_mm = gnc_pp_Mult_mm;
239
240  p_Procs->p_Minus_mm_Mult_qq = rGR->p_Procs->p_Minus_mm_Mult_qq = NULL;
241
242  // non-commutaitve multiplication by monomial from the left
243  rGR->GetNC()->p_Procs.mm_Mult_p   = gnc_mm_Mult_p;
244  rGR->GetNC()->p_Procs.mm_Mult_pp  = gnc_mm_Mult_pp;
245
246}
247
248bool ncInitSpecialPairMultiplication(ring r)
249{
250#if OUTPUT 
251  Print("ncInitSpecialPairMultiplication(ring), ring: \n");
252  rWrite(r);
253  PrintLn();
254#endif
255 
256  assume(rIsPluralRing(r));
257  assume(!rIsSCA(r));
258
259  if( r->GetNC()->GetGlobalMultiplier() != NULL )
260  {
261    WarnS("Already defined!");
262    return false;
263  }
264
265  r->GetNC()->GetGlobalMultiplier() = new CGlobalMultiplier(r);
266
267  gnc_p_ProcsSet(r);
268  return true;
269}
270
271
272CGlobalMultiplier::CGlobalMultiplier(ring r):
273    CMultiplier<poly>(r), m_RingFormulaMultiplier(GetFormulaPowerMultiplier(r))
274{
275#if OUTPUT 
276  Print("CGlobalMultiplier::CGlobalMultiplier(ring)!");
277  PrintLn();
278#endif
279
280  m_cache = new CGlobalCacheHash(r);
281  m_powers = new CPowerMultiplier(r);
282}
283
284
285CGlobalMultiplier::~CGlobalMultiplier()
286{
287#if OUTPUT 
288  Print("CGlobalMultiplier::~CGlobalMultiplier()!");
289  PrintLn();
290#endif
291
292  delete m_cache;
293  delete m_powers;
294
295  // we cannot delete m_RingFormulaMultiplier as it belongs to the ring!
296}
297
298
299
300// Exponent * Exponent
301// TODO: handle components!!!
302poly CGlobalMultiplier::MultiplyEE(const CGlobalMultiplier::CExponent expLeft, const CGlobalMultiplier::CExponent expRight)
303{
304
305  const ring r = GetBasering();
306
307#if OUTPUT 
308  Print("CGlobalMultiplier::MultiplyEE(expLeft, expRight)!");
309  PrintLn();
310  PrintS("expL: "); p_Write(expLeft, GetBasering());   
311  PrintS("expR: "); p_Write(expRight, GetBasering());   
312#endif
313
314  CCacheHash<poly>::CCacheItem* pLookup;
315 
316  int b = m_cache->LookupEE(expLeft, expRight, pLookup);
317  // TODO!!!
318
319  // up to now:
320  assume( b == -1 );
321
322  // TODO: use PowerMultiplier!!!!
323
324  poly product = NULL;
325
326  const int N = NVars();
327  int j = N;
328  int i = 1;
329
330  int ej = p_GetExp(expLeft, j, r);
331  int ei = p_GetExp(expRight, i, r);
332
333  while( (i < j) && !((ej != 0) && (ei != 0)) )
334  {
335    if( ei == 0 )
336      ei = p_GetExp(expRight, ++i, r);
337   
338    if( ej == 0 )
339      ej = p_GetExp(expLeft, --j, r);
340  }
341
342 
343#if OUTPUT 
344  Print("<CGlobalMultiplier::MultiplyEE>");
345  PrintLn();
346  Print("i: %d, j: %d", i, j); 
347  PrintLn();
348  Print("ei: %d, ej: %d", ei, ej); 
349  PrintLn();
350#endif
351
352
353  // |  expLeft   | * |  expRight  |
354  // |<<<< ej 0..0| , |0..0 ei >>>>|
355  // |<<<<  j <<<N| , |1>>>  i >>>>|
356
357  if( i >= j ) // BUG here!!!???
358  {
359    // either i == j or i = j + 1 => commutative multiple!
360    // TODO: it can be done more efficiently! ()
361    product = p_Head(expRight, r);
362
363  // |  expLeft     | * |  expRight   |
364  // |<<<< ej 0....0| , |0..00 ei >>>>|
365  // |<<<<  j i <<<N| , |1>>>j  i >>>>|
366
367    if(i > j)
368    {
369      --i;
370      ei = 0;
371    }
372     
373    if( i == j )
374    {
375      if( ej != 0 )
376        p_SetExp(product, i, ei + ej, r);
377    }
378
379    --i;
380
381    for(; i > 0; --i)
382    {
383      const int e = p_GetExp(expLeft, i, r);
384
385      if( e > 0 )
386        p_SetExp(product, i, e, r);
387    }
388
389    p_Setm(product, r);   
390
391  } else
392  { // i < j, ei != 0, ej != 0
393
394    Enum_ncSAType PairType = _ncSA_notImplemented;
395
396    if( m_RingFormulaMultiplier != NULL )
397      PairType = m_RingFormulaMultiplier->GetPair(i, j);
398
399
400    if( PairType == _ncSA_notImplemented )
401      product = m_powers->MultiplyEE( CPower(j, ej), CPower(i, ei) );
402//    return gnc_uu_Mult_ww_vert(i, a, j, b, r);
403    else
404 //    return m_RingFormulaMultiplier->Multiply(j, i, b, a);
405      product = CFormulaPowerMultiplier::Multiply( PairType, i, j, ei, ej, GetBasering());
406   
407
408#if OUTPUT 
409    Print("<CGlobalMultiplier::MultiplyEE> ==> ");
410    PrintLn();
411    Print("i: %d, j: %d", i, j); 
412    PrintLn();
413    Print("ei: %d, ej: %d", ei, ej); 
414    PrintLn();
415    Print("<product>: "); p_Write(product, GetBasering()); 
416#endif
417   
418
419    // TODO: Choose some multiplication strategy!!!
420   
421    while( (product != NULL) && !((i == NVars()) && (j == 1)) )
422    {
423
424      // make some choice here!:
425
426      if( i < NVars() )
427      {
428        ei = p_GetExp(expRight, ++i, r);
429       
430        while( (ei == 0) && (i < NVars()) )
431          ei = p_GetExp(expRight, ++i, r);
432
433        if( ei != 0 )
434          product = m_powers->MultiplyPEDestroy(product, CPower(i, ei));
435      } 
436
437      if( j > 1 )
438      {
439        ej = p_GetExp(expLeft, --j, r);
440
441        while( (ej == 0) && (1 < j) )
442          ej = p_GetExp(expLeft, --j, r);
443
444        if( ej != 0 )
445          product = m_powers->MultiplyEPDestroy(CPower(j, ej), product);
446      }
447
448
449#if OUTPUT 
450      Print("<CGlobalMultiplier::MultiplyEE> ==> ");
451      PrintLn();
452      Print("i: %d, j: %d", i, j); 
453      PrintLn();
454      Print("ei: %d, ej: %d", ei, ej); 
455      PrintLn();
456      Print("<product>: "); p_Write(product, GetBasering()); 
457#endif
458     
459    }
460
461  }
462
463  // TODO!     
464 
465
466  m_cache->StoreEE( expLeft, expRight, product);
467  // up to now:
468  return product; 
469}
470
471    // Monom * Exponent
472poly CGlobalMultiplier::MultiplyME(const poly pMonom, const CGlobalMultiplier::CExponent expRight)
473{
474#if OUTPUT 
475  Print("CGlobalMultiplier::MultiplyME(monom, expR)!"); 
476  PrintLn();
477  PrintS("Monom: "); p_Write(pMonom, GetBasering());   
478  PrintS("expR: "); p_Write(expRight, GetBasering());   
479#endif
480
481  return MultiplyEE(pMonom, expRight);
482}
483
484    // Exponent * Monom
485poly CGlobalMultiplier::MultiplyEM(const CGlobalMultiplier::CExponent expLeft, const poly pMonom)
486{
487#if OUTPUT 
488  Print("CGlobalMultiplier::MultiplyEM(expL, monom)!"); 
489  PrintLn();
490  PrintS("expL: "); p_Write(expLeft, GetBasering());   
491  PrintS("Monom: "); p_Write(pMonom, GetBasering());   
492#endif
493
494  return MultiplyEE(expLeft, pMonom);
495}
496
497
498
499
500
501///////////////////////////////////////////////////////////////////////////////////////////
502CCommutativeSpecialPairMultiplier::CCommutativeSpecialPairMultiplier(ring r, int i, int j):
503    CSpecialPairMultiplier(r, i, j)
504{
505#if OUTPUT 
506  Print("CCommutativeSpecialPairMultiplier::CCommutativeSpecialPairMultiplier(ring, i: %d, j: %d)!", i, j);
507  PrintLn();
508#endif
509};
510
511
512CCommutativeSpecialPairMultiplier::~CCommutativeSpecialPairMultiplier()
513{
514#if OUTPUT 
515  Print("CCommutativeSpecialPairMultiplier::~CCommutativeSpecialPairMultiplier()");
516  PrintLn();
517#endif
518}
519
520// Exponent * Exponent
521poly CCommutativeSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
522{
523#if OUTPUT 
524  Print("CCommutativeSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
525  PrintLn();
526#endif
527
528  const ring r = GetBasering();
529
530  return CFormulaPowerMultiplier::ncSA_1xy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
531}
532
533///////////////////////////////////////////////////////////////////////////////////////////
534CAntiCommutativeSpecialPairMultiplier::CAntiCommutativeSpecialPairMultiplier(ring r, int i, int j):
535                CSpecialPairMultiplier(r, i, j)
536{
537#if OUTPUT 
538        Print("CAntiCommutativeSpecialPairMultiplier::CAntiCommutativeSpecialPairMultiplier(ring, i: %d, j: %d)!", i, j);
539        PrintLn();
540#endif
541};
542
543
544CAntiCommutativeSpecialPairMultiplier::~CAntiCommutativeSpecialPairMultiplier()
545{
546#if OUTPUT 
547        Print("CAntiCommutativeSpecialPairMultiplier::~CAntiCommutativeSpecialPairMultiplier()");
548        PrintLn();
549#endif
550}
551
552// Exponent * Exponent
553poly CAntiCommutativeSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
554{
555#if OUTPUT 
556        Print("CAntiCommutativeSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
557        PrintLn();
558#endif
559
560        const ring r = GetBasering();
561
562        return CFormulaPowerMultiplier::ncSA_Mxy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
563}
564
565///////////////////////////////////////////////////////////////////////////////////////////
566CQuasiCommutativeSpecialPairMultiplier::CQuasiCommutativeSpecialPairMultiplier(ring r, int i, int j, number q):
567                CSpecialPairMultiplier(r, i, j), m_q(q)
568{
569#if OUTPUT 
570        Print("CQuasiCommutativeSpecialPairMultiplier::CQuasiCommutativeSpecialPairMultiplier(ring, i: %d, j: %d, q)!", i, j);
571        PrintLn();
572        PrintS("Parameter q: ");
573        n_Write(q, r);
574#endif
575};
576
577
578CQuasiCommutativeSpecialPairMultiplier::~CQuasiCommutativeSpecialPairMultiplier()
579{
580#if OUTPUT 
581        Print("CQuasiCommutativeSpecialPairMultiplier::~CQuasiCommutativeSpecialPairMultiplier()");
582        PrintLn();
583#endif
584}
585
586// Exponent * Exponent
587poly CQuasiCommutativeSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
588{
589#if OUTPUT 
590        Print("CQuasiCommutativeSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
591        PrintLn();
592#endif
593
594        const ring r = GetBasering();
595
596        return CFormulaPowerMultiplier::ncSA_Qxy0x0y0(GetI(), GetJ(), expRight, expLeft, m_q, r);
597}
598
599
600///////////////////////////////////////////////////////////////////////////////////////////
601CWeylSpecialPairMultiplier::CWeylSpecialPairMultiplier(ring r, int i, int j, number g):
602    CSpecialPairMultiplier(r, i, j), m_g(g)
603{
604#if OUTPUT 
605  Print("CWeylSpecialPairMultiplier::CWeylSpecialPairMultiplier(ring, i: %d, j: %d, g)!", i, j);
606  PrintLn();
607  PrintS("Parameter g: ");
608  n_Write(g, r);
609#endif
610};
611
612
613CWeylSpecialPairMultiplier::~CWeylSpecialPairMultiplier()
614{
615#if OUTPUT 
616  Print("CWeylSpecialPairMultiplier::~CWeylSpecialPairMultiplier()");
617  PrintLn();
618#endif
619}
620
621// Exponent * Exponent
622poly CWeylSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
623{
624#if OUTPUT 
625  Print("CWeylSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
626  PrintLn();
627#endif
628  // Char == 0, otherwise - problem!
629
630 
631  const ring r = GetBasering();
632
633  assume( expLeft*expRight > 0 );
634
635  return CFormulaPowerMultiplier::ncSA_1xy0x0yG(GetI(), GetJ(), expRight, expLeft, m_g, r);
636}
637
638
639///////////////////////////////////////////////////////////////////////////////////////////
640CShiftSpecialPairMultiplier::CShiftSpecialPairMultiplier(ring r, int i, int j, int s, number c):
641    CSpecialPairMultiplier(r, i, j), m_shiftCoef(c), m_shiftVar(s)
642{
643#if OUTPUT 
644  Print("CShiftSpecialPairMultiplier::CShiftSpecialPairMultiplier(ring, i: %d, j: %d, s: %d, c)!", i, j, s);
645  PrintLn();
646  PrintS("Parameter c: "); n_Write(c, r);
647#endif
648};
649
650
651CShiftSpecialPairMultiplier::~CShiftSpecialPairMultiplier()
652{
653#if OUTPUT 
654  Print("CShiftSpecialPairMultiplier::~CShiftSpecialPairMultiplier()");
655  PrintLn();
656#endif
657}
658
659// Exponent * Exponent
660poly CShiftSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
661{
662#if OUTPUT 
663  Print("CShiftSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
664  PrintLn();
665#endif
666  // Char == 0, otherwise - problem!
667
668  assume( expLeft*expRight > 0 );
669
670  const ring r = GetBasering();
671
672  if( m_shiftVar != GetI() ) // YX = XY + b*Y?
673    return CFormulaPowerMultiplier::ncSA_1xy0xBy0(GetI(), GetJ(), expRight, expLeft, m_shiftCoef, r); // case: (1, 0, beta, 0, 0)
674  else
675    return CFormulaPowerMultiplier::ncSA_1xyAx0y0(GetI(), GetJ(), expRight, expLeft, m_shiftCoef, r); // case: (1, alpha, 0, 0)
676
677}
678
679
680
681///////////////////////////////////////////////////////////////////////////////////////////
682CExternalSpecialPairMultiplier::CExternalSpecialPairMultiplier(ring r, int i, int j, Enum_ncSAType type):
683    CSpecialPairMultiplier(r, i, j), m_ncSAtype(type)
684{
685#if OUTPUT 
686  Print("CExternalSpecialPairMultiplier::CExternalSpecialPairMultiplier(ring, i: %d, j: %d, type: %d, c)!", i, j, (int)type);
687  PrintLn();
688#endif
689};
690
691
692CExternalSpecialPairMultiplier::~CExternalSpecialPairMultiplier()
693{
694#if OUTPUT 
695  Print("CExternalSpecialPairMultiplier::~CExternalSpecialPairMultiplier()");
696  PrintLn();
697#endif
698}
699
700// Exponent * Exponent
701poly CExternalSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
702{
703#if OUTPUT 
704  Print("CExternalSpecialPairMultiplier::MultiplyEE(var(%d)^{%d}, var(%d)^{%d})!", GetJ(), expLeft, GetI(), expRight); 
705  PrintLn();
706#endif
707  // Char == 0, otherwise - problem!
708
709  assume( expLeft*expRight > 0 );
710
711  const ring r = GetBasering();
712
713  return CFormulaPowerMultiplier::Multiply(m_ncSAtype, GetI(), GetJ(), expRight, expLeft, r); 
714
715}
716
717
718
719///////////////////////////////////////////////////////////////////////////////////////////
720
721// factory method!
722CSpecialPairMultiplier* AnalyzePair(const ring r, int i, int j)
723{
724#if OUTPUT 
725  Print("AnalyzePair(ring, i: %d, j: %d)!", i, j);
726  PrintLn();
727#endif
728
729  Enum_ncSAType type = CFormulaPowerMultiplier::AnalyzePair(r, i, j);
730
731  if( type == _ncSA_notImplemented ) return NULL;
732
733
734  // last possibility:
735  return new CExternalSpecialPairMultiplier(r, i, j, type); // For tests!
736
737 
738  if( type == _ncSA_1xy0x0y0 )
739    return new CCommutativeSpecialPairMultiplier(r, i, j);
740
741  if( type == _ncSA_Mxy0x0y0 )
742    return new CAntiCommutativeSpecialPairMultiplier(r, i, j);
743
744  if( type == _ncSA_Qxy0x0y0 )
745  {
746    const number q = p_GetCoeff(GetC(r, i, j), r);
747    return new CQuasiCommutativeSpecialPairMultiplier(r, i, j, q);
748  }
749 
750  const poly d = GetD(r, i, j);
751 
752  assume( d != NULL );
753  assume( pNext(d) == NULL );
754
755  const number g = p_GetCoeff(d, r);
756
757  if( type == _ncSA_1xy0x0yG ) // Weyl
758    return new CWeylSpecialPairMultiplier(r, i, j, g);         
759
760  if( type == _ncSA_1xyAx0y0 ) // Shift 1
761    return new CShiftSpecialPairMultiplier(r, i, j, i, g);
762
763  if( type == _ncSA_1xy0xBy0 ) // Shift 2
764    return new CShiftSpecialPairMultiplier(r, i, j, j, g);
765
766}
767
768
769
770
771
772
773CPowerMultiplier::CPowerMultiplier(ring r): CMultiplier<CPower>(r)
774{
775#if OUTPUT 
776  Print("CPowerMultiplier::CPowerMultiplier(ring)!");
777  PrintLn();
778#endif
779
780  m_specialpairs = (CSpecialPairMultiplier**)omAlloc0( ((NVars() * (NVars()-1)) / 2) * sizeof(CSpecialPairMultiplier*) );
781
782  for( int i = 1; i < NVars(); i++ )
783    for( int j = i + 1; j <= NVars(); j++ )
784      GetPair(i, j) = AnalyzePair(GetBasering(), i, j); // factory method!
785}
786
787
788CPowerMultiplier::~CPowerMultiplier()
789{
790#if OUTPUT 
791  Print("CPowerMultiplier::~CPowerMultiplier()!");
792  PrintLn();
793#endif
794
795  omFreeSize((ADDRESS)m_specialpairs, ((NVars() * (NVars()-1)) / 2) * sizeof(CSpecialPairMultiplier*) );
796}
797
798
799// Monom * Exponent
800// pMonom may NOT be of the form: var(j)^{n}!
801poly CPowerMultiplier::MultiplyME(const poly pMonom, const CExponent expRight)
802{
803  const int j = expRight.Var;
804  const int n = expRight.Power;
805
806  const ring r = GetBasering();
807 
808#if OUTPUT 
809  Print("CPowerMultiplier::MultiplyME(monom * var(%d)^{%d})!", j, n);
810  PrintLn();
811  PrintS("Monom: "); p_Write(pMonom, r); 
812#endif
813
814  assume( (j > 0) && (j <= NVars()));
815
816  if( n == 0 )
817    return p_Head(pMonom, r); // Copy?!?
818 
819
820  int v = NVars();
821  int e = p_GetExp(pMonom, v, r);
822
823  while((v > j) && (e == 0))
824    e = p_GetExp(pMonom, --v, r);
825
826  // TODO: review this!
827  if( (v == j) )
828  {
829    poly p = p_Head(pMonom, r);   
830    p_SetExp(p, v, e + n, r);
831    p_Setm(p, r);   
832
833    return p;
834  }
835
836  assume( v > j );
837  assume( e > 0 );
838 
839  // And now the General Case: v > j!
840
841  poly p = MultiplyEE( CPower(v, e), expRight ); // Easy way!
842
843  --v;
844 
845  while(v > 0)
846  {
847    e = p_GetExp(pMonom, v, GetBasering());
848   
849    if( e > 0 )
850      p = MultiplyEPDestroy(CPower(v, e), p);
851
852    --v;
853  }
854
855#if OUTPUT 
856  Print("CPowerMultiplier::MultiplyME() ===> ");
857  p_Write(p, GetBasering()); 
858#endif
859 
860  return p;
861}
862
863// Exponent * Monom
864// pMonom may NOT be of the form: var(i)^{m}!
865poly CPowerMultiplier::MultiplyEM(const CExponent expLeft, const poly pMonom)
866{
867  const ring r = GetBasering();
868
869  // TODO: as above! (difference due to Left/Right semmantics!)
870  const int j = expLeft.Var;
871  const int n = expLeft.Power;
872
873#if OUTPUT 
874  Print("CPowerMultiplier::MultiplyEM(var(%d)^{%d} * monom)!", j, n);
875  PrintLn();
876  PrintS("Monom: "); p_Write(pMonom, r); 
877#endif
878
879  assume( (j > 0) && (j <= NVars()));
880
881  if( n == 0 )
882    return p_Head(pMonom, r); // Copy?!?
883
884
885  int v = 1; // NVars();
886  int e = p_GetExp(pMonom, v, r);
887
888  while((v < j) && (e == 0))
889    e = p_GetExp(pMonom, ++v, r);
890
891  if( v == j ) 
892  {
893    poly p = p_Head(pMonom, r);   
894    p_SetExp(p, j, e + n, r);
895    p_Setm(p, r);   
896
897    return p;
898  }
899
900  assume( v < j );
901  assume( e > 0 );
902
903 
904  // And now the General Case: v > j!
905
906  poly p = MultiplyEE( expLeft, CPower(v, e) ); // Easy way!
907
908  ++v;
909
910  while(v <= NVars())
911  {
912    e = p_GetExp(pMonom, v, r);
913   
914    if( e > 0 )
915      p = MultiplyPEDestroy(p, CPower(v, e));
916         
917    ++v;
918  }
919
920#if OUTPUT 
921  Print("CPowerMultiplier::MultiplyEM() ===> ");
922  p_Write(p, r); 
923#endif
924
925  return p;
926 
927}
928
929
930// Exponent * Exponent
931// Computes: var(j)^{expLeft} * var(i)^{expRight}
932poly CPowerMultiplier::MultiplyEE(const CExponent expLeft, const CExponent expRight)
933{
934#if OUTPUT 
935  Print("CPowerMultiplier::MultiplyEE)!");
936  PrintLn();
937#endif
938
939  const int i = expRight.Var, j = expLeft.Var;
940  const int ei = expRight.Power, ej = expLeft.Power;
941
942#if OUTPUT 
943  Print("Input: var(%d)^{%d} * var(%d)^{%d}", j, ej, i, ei);
944  PrintLn();
945#endif
946
947  assume(1 <= i);
948  assume(j <= NVars());
949  assume(1 <= j);
950  assume(i <= NVars());
951  assume(ei > 0);
952  assume(ej > 0);
953 
954  if( i >= j )
955  {
956    const ring r = GetBasering();
957
958    poly product = p_One(r);
959    p_SetExp(product, j, ej, r);
960    p_SetExp(product, i, ei, r);
961    p_Setm(product, r);
962
963    return product;
964
965  } else
966  {
967    assume(i <  j);
968
969    // No Cache Lookup!? :(
970
971    CSpecialPairMultiplier* pSpecialMultiplier = GetPair(i, j);
972
973    poly product = NULL;
974
975    // Special case?
976    if( pSpecialMultiplier != NULL )
977    {
978      assume( pSpecialMultiplier->GetI() == i );
979      assume( pSpecialMultiplier->GetJ() == j );
980      assume( pSpecialMultiplier->GetBasering() == GetBasering() );
981
982      return pSpecialMultiplier->MultiplyEE(ej, ei);
983    } else
984    {
985      // Perform general NC Multiplication:
986      // TODO
987
988//      product = NULL;
989    }
990  }
991 
992  return NULL; 
993}
994
995
996
997
998
999
1000CSpecialPairMultiplier::CSpecialPairMultiplier(ring r, int i, int j):
1001    CMultiplier<int>(r), m_i(i), m_j(j)
1002{
1003#if OUTPUT 
1004  Print("CSpecialPairMultiplier::CSpecialPairMultiplier(ring, i: %d, j: %d)!", i, j);
1005  PrintLn();
1006#endif
1007 
1008  assume(i < j);
1009  assume(i > 0);
1010  assume(j <= NVars());
1011}
1012
1013
1014CSpecialPairMultiplier::~CSpecialPairMultiplier()
1015{
1016#if OUTPUT 
1017  Print("CSpecialPairMultiplier::~CSpecialPairMultiplier()!");
1018  PrintLn();
1019#endif
1020}
1021
1022
1023
1024// Monom * Exponent
1025poly CSpecialPairMultiplier::MultiplyME(const poly pMonom, const CExponent expRight)
1026{
1027#if OUTPUT 
1028  Print("CSpecialPairMultiplier::MultiplyME(monom, var(%d)^{%d})!", GetI(), expRight); 
1029  PrintLn();
1030  PrintS("Monom: "); p_Write(pMonom, GetBasering());
1031#endif
1032 
1033  return MultiplyEE(p_GetExp(pMonom, GetJ(), GetBasering()), expRight);
1034}
1035
1036    // Exponent * Monom
1037poly CSpecialPairMultiplier::MultiplyEM(const CExponent expLeft, const poly pMonom)
1038{
1039#if OUTPUT 
1040  Print("CSpecialPairMultiplier::MultiplyEM(var(%d)^{%d}, monom)!", GetJ(), expLeft); 
1041  PrintLn();
1042  PrintS("Monom: "); p_Write(pMonom, GetBasering());
1043#endif
1044
1045  return MultiplyEE(expLeft, p_GetExp(pMonom, GetI(), GetBasering()));
1046}
1047
1048
1049
1050
Note: See TracBrowser for help on using the repository browser.