source: git/kernel/ncSAMult.cc @ 1495df4

spielwiese
Last change on this file since 1495df4 was 1495df4, checked in by Motsak Oleksandr <motsak@…>, 16 years ago
*motsak: special multiplication git-svn-id: file:///usr/local/Singular/svn/trunk@10868 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.1 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.2 2008-07-15 16:27:58 motsak Exp $
10 *******************************************************************/
11
12
13#define MYTEST 1
14#define OUTPUT 1
15
16#if MYTEST
17#define OM_CHECK 4
18#define OM_TRACK 5
19#endif
20
21#include "mod2.h"
22
23#include <ncSAMult.h> // for CMultiplier etc classes
24#include <sca.h> // for SCA
25
26
27
28
29// poly functions defined in p_Procs: ;
30static poly gnc_pp_Mult_mm(const poly p, const poly m, const ring r, poly& last)
31{
32#if OUTPUT 
33  Print("gnc_pp_Mult_mm");
34  PrintLn();
35#endif
36  assume( r->GetNC()->GetGlobalMultiplier() != NULL );
37 
38  return r->GetNC()->GetGlobalMultiplier()->MultiplyPE(p, m);
39}
40
41static poly gnc_p_Mult_mm(poly p, const poly m, const ring r)
42{
43#if OUTPUT 
44  Print("gnc_p_Mult_mm");
45  PrintLn();
46#endif
47  assume( r->GetNC()->GetGlobalMultiplier() != NULL );
48
49  return r->GetNC()->GetGlobalMultiplier()->MultiplyPEDestroy(p, m);
50}
51
52static poly gnc_mm_Mult_p(const poly m, poly p, const ring r)
53{
54#if OUTPUT 
55  Print("gnc_mm_Mult_p");
56  PrintLn();
57#endif
58  assume( r->GetNC()->GetGlobalMultiplier() != NULL );
59
60  return r->GetNC()->GetGlobalMultiplier()->MultiplyEPDestroy(m, p);
61}
62
63static poly gnc_mm_Mult_pp(const poly m, const poly p, const ring r)
64{
65#if OUTPUT 
66  Print("gnc_mm_Mult_pp");
67  PrintLn();
68#endif
69 
70  assume( r->GetNC()->GetGlobalMultiplier() != NULL );
71
72  return r->GetNC()->GetGlobalMultiplier()->MultiplyEP(m, p);
73}
74
75
76
77
78static void gnc_p_ProcsSet(ring rGR, p_Procs_s* p_Procs = NULL)
79{
80#if OUTPUT 
81  Print("\ngnc_p_ProcsSet()!!!");
82  PrintLn();
83#endif
84
85  if( p_Procs == NULL )
86    p_Procs = rGR->p_Procs;
87 
88  // "commutative"
89  p_Procs->p_Mult_mm  = rGR->p_Procs->p_Mult_mm  = gnc_p_Mult_mm;
90  p_Procs->pp_Mult_mm = rGR->p_Procs->pp_Mult_mm = gnc_pp_Mult_mm;
91
92  p_Procs->p_Minus_mm_Mult_qq = rGR->p_Procs->p_Minus_mm_Mult_qq = NULL;
93
94  // non-commutaitve multiplication by monomial from the left
95  rGR->GetNC()->p_Procs.mm_Mult_p   = gnc_mm_Mult_p;
96  rGR->GetNC()->p_Procs.mm_Mult_pp  = gnc_mm_Mult_pp;
97
98}
99
100
101
102
103
104
105
106bool ncInitSpecialPairMultiplication(ring r)
107{
108#if OUTPUT 
109  Print("ncInitSpecialPairMultiplication(ring), ring: \n");
110  rWrite(r);
111  PrintLn();
112#endif
113  assume(rIsPluralRing(r));
114  assume(!rIsSCA(r));
115
116  r->GetNC()->GetGlobalMultiplier() = new CGlobalMultiplier(r);
117
118  gnc_p_ProcsSet(r);
119  return true;
120}
121
122
123CGlobalMultiplier::CGlobalMultiplier(ring r): CMultiplier<poly>(r)
124{
125#if OUTPUT 
126  Print("CGlobalMultiplier::CGlobalMultiplier(ring)!");
127  PrintLn();
128#endif
129
130  m_cache = new CGlobalCacheHash(r);
131  m_powers = new CPowerMultiplier(r);
132
133}
134
135
136CGlobalMultiplier::~CGlobalMultiplier()
137{
138#if OUTPUT 
139  Print("CGlobalMultiplier::~CGlobalMultiplier()!");
140  PrintLn();
141#endif
142
143  delete m_cache;
144  delete m_powers;
145}
146
147
148
149// Exponent * Exponent
150// TODO: handle components!!!
151poly CGlobalMultiplier::MultiplyEE(const CExponent expLeft, const CExponent expRight)
152{
153#if OUTPUT 
154  Print("CGlobalMultiplier::MultiplyEE(expLeft, expRight)!");
155  PrintLn();
156#endif
157
158  CCacheHash<poly>::CCacheItem* pLookup;
159  int b = m_cache->LookupEE(expLeft, expRight, pLookup);
160
161  // up to now:
162  assume( b == -1 );
163
164  // TODO: use PowerMultiplier!!!!
165 
166 
167  // up to now:
168  return NULL; 
169}
170
171    // Monom * Exponent
172poly CGlobalMultiplier::MultiplyME(const poly pMonom, const CExponent expRight)
173{
174  return MultiplyEE(pMonom, expRight);
175}
176
177    // Exponent * Monom
178poly CGlobalMultiplier::MultiplyEM(const CExponent expLeft, const poly pMonom)
179{
180  return MultiplyEE(expLeft, pMonom);
181}
182
183
184
185
186
187// Exponent * Exponent
188poly CCommutativeSpecialPairMultiplier::MultiplyEE(const int expLeft, const int expRight)
189{
190  return NULL;
191}
192
193// Monom * Exponent
194poly CCommutativeSpecialPairMultiplier::MultiplyME(const poly pMonom, const int expRight)
195{
196  return NULL;
197}
198
199// Exponent * Monom
200poly CCommutativeSpecialPairMultiplier::MultiplyEM(const int expLeft, const poly pMonom)
201{
202  return NULL;
203}
204
205
206
207// factory method!
208CSpecialPairMultiplier* AnalyzePair(const ring r, int i, int j)
209{
210#if OUTPUT 
211  Print("AnalyzePair(ring, i: %d, j: %d)!", i, j);
212  PrintLn();
213#endif
214
215  const int N = r->N;
216
217  assume(i < j);
218  assume(i > 0);
219  assume(j <= N);
220
221
222  const poly c = GetC(r, i, j);
223  const poly d = GetD(r, i, j);
224
225#if OUTPUT 
226  Print("C_{%d, %d} = ", i, j); p_Write(c, r);
227  Print("D_{%d, %d} = ", i, j); p_Write(d, r);
228#endif
229
230  if( (d == NULL) && n_IsOne(p_GetCoeff(c, r), r) )
231    return new CCommutativeSpecialPairMultiplier(r, i, j);
232
233  return NULL;
234}
235
236
237
238
239
240
241CPowerMultiplier::CPowerMultiplier(ring r): CMultiplier<CPower>(r)
242{
243#if OUTPUT 
244  Print("CPowerMultiplier::CPowerMultiplier(ring)!");
245  PrintLn();
246#endif
247
248  m_specialpairs = (CSpecialPairMultiplier**)omAlloc0( ((NVars() * (NVars()-1)) / 2) * sizeof(CSpecialPairMultiplier*) );
249
250  for( int i = 1; i < NVars(); i++ )
251    for( int j = i + 1; j <= NVars(); j++ )
252      GetPair(i, j) = AnalyzePair(GetBasering(), i, j); // factory method!
253}
254
255
256CPowerMultiplier::~CPowerMultiplier()
257{
258#if OUTPUT 
259  Print("CPowerMultiplier::~CPowerMultiplier()!");
260  PrintLn();
261#endif
262
263  omFreeSize((ADDRESS)m_specialpairs, ((NVars() * (NVars()-1)) / 2) * sizeof(CSpecialPairMultiplier*) );
264}
265
266
267// Monom * Exponent
268// pMonom may NOT be of the form: var(j)^{n}!
269poly CPowerMultiplier::MultiplyME(const poly pMonom, const CExponent expRight)
270{
271  const int j = expRight.Var;
272  const int n = expRight.Power;
273 
274#if OUTPUT 
275  Print("CPowerMultiplier::MultiplyME(monom * var(%d)^{%d}!", j, n);
276  PrintLn();
277#endif
278
279  assume( (j > 0) && (j <= NVars()));
280
281  if( n == 0 )
282    return p_Copy(pMonom, GetBasering()); // Copy?!?
283 
284
285  int v = NVars();
286  int e = p_GetExp(pMonom, v, GetBasering());
287
288  while((v > j) && (e == 0))
289    e = p_GetExp(pMonom, --v, GetBasering());
290
291  // TODO: review this!
292  if( (e == 0) )
293    return expRight.GetPoly(GetBasering());
294
295  if( v == j ) 
296  {
297    poly p = p_Copy(pMonom, GetBasering());   
298    p_SetExp(p, j, e + n, GetBasering());
299    return p;
300  }
301 
302  // And now the General Case: v > j!
303
304  poly p = MultiplyEE( CPower(v, e), expRight ); // Easy way!
305
306  --v;
307 
308  while(v > 0)
309  {
310    e = p_GetExp(pMonom, v, GetBasering());
311   
312    if( e > 0 )
313      p = MultiplyEPDestroy(CPower(v, e), p);
314
315    --v;
316  }
317
318  return p;
319}
320
321// Exponent * Monom
322// pMonom may NOT be of the form: var(i)^{m}!
323poly CPowerMultiplier::MultiplyEM(const CExponent expLeft, const poly pMonom)
324{
325  // TODO: as above! (difference due to Left/Right semmantics!)
326  const int j = expLeft.Var;
327  const int n = expLeft.Power;
328
329#if OUTPUT 
330  Print("CPowerMultiplier::MultiplyEM(var(%d)^{%d} * monom)!", j, n);
331  PrintLn();
332#endif
333
334  assume( (j > 0) && (j <= NVars()));
335
336  if( n == 0 )
337    return p_Copy(pMonom, GetBasering()); // Copy?!?
338
339
340  int v = 1; // NVars();
341  int e = p_GetExp(pMonom, v, GetBasering());
342
343  while((v < j) && (e == 0))
344    e = p_GetExp(pMonom, ++v, GetBasering());
345
346  // TODO: review this!
347  if( (e == 0) )
348    return expLeft.GetPoly(GetBasering());
349
350  if( v == j ) 
351  {
352    poly p = p_Copy(pMonom, GetBasering());   
353    p_SetExp(p, j, e + n, GetBasering());
354    return p;
355  }
356
357  // And now the General Case: v > j!
358
359  poly p = MultiplyEE( expLeft, CPower(v, e) ); // Easy way!
360
361  ++v;
362  while(v <= NVars())
363  {
364    e = p_GetExp(pMonom, v, GetBasering());
365   
366    if( e > 0 )
367      p = MultiplyPEDestroy(p, CPower(v, e));
368         
369    ++v;
370  }
371 
372}
373
374
375// Exponent * Exponent
376// Computes: var(j)^{expLeft} * var(i)^{expRight}
377poly CPowerMultiplier::MultiplyEE(const CExponent expLeft, const CExponent expRight)
378{
379#if OUTPUT 
380  Print("CPowerMultiplier::MultiplyEE)!");
381  PrintLn();
382#endif
383
384  const int i = expRight.Var, j = expLeft.Var;
385  const int ei = expRight.Power, ej = expLeft.Power;
386
387#if OUTPUT 
388  Print("Input: var(%d)^{%d} * var(%d)^{%d}", j, ej, i, ei);
389  PrintLn();
390#endif
391
392  poly product;
393
394  if( i >= j )
395  {
396    // easy!
397    product = NULL;
398  } else
399  {
400    assume(1 <= i);
401    assume(i <  j);
402    assume(j <= NVars());
403    assume(ei > 0);
404    assume(ej > 0);
405
406    // No Cache Lookup!? :(
407
408    CSpecialPairMultiplier* pSpecialMultiplier = GetPair(i, j);
409
410    poly product = NULL;
411
412    // Special case?
413    if( pSpecialMultiplier != NULL )
414    {
415      assume( pSpecialMultiplier->GetI() == i );
416      assume( pSpecialMultiplier->GetJ() == j );
417      assume( pSpecialMultiplier->GetBasering() == GetBasering() );
418
419      product = pSpecialMultiplier->MultiplyEE(ej, ei);
420    } else
421    {
422      // Perform general NC Multiplication:
423      // TODO
424
425      product = NULL;
426    }
427  }
428 
429  return product; 
430}
431
432
433
434
435
436
437CSpecialPairMultiplier::CSpecialPairMultiplier(ring r, int i, int j):
438    CMultiplier<int>(r), m_i(i), m_j(j)
439{
440  assume(i < j);
441  assume(i > 0);
442  assume(j <= NVars());
443}
444
445
446// Monom * Exponent
447poly CSpecialPairMultiplier::MultiplyME(const poly pMonom, const CExponent expRight)
448{
449  return MultiplyEE(p_GetExp(pMonom, GetJ(), GetBasering()), expRight);
450}
451
452    // Exponent * Monom
453poly CSpecialPairMultiplier::MultiplyEM(const CExponent expLeft, const poly pMonom)
454{
455  return MultiplyEE(expLeft, p_GetExp(pMonom, GetI(), GetBasering()));
456}
457
458
459
460
Note: See TracBrowser for help on using the repository browser.