source: git/libpolys/polys/nc/ncSAFormula.cc @ 76fead

spielwiese
Last change on this file since 76fead was a60e0b, checked in by Mohamed Barakat <mohamed.barakat@…>, 13 years ago
fixed the syntax for compilers in pedantic mood :)
  • Property mode set to 100644
File size: 13.0 KB
Line 
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 *  Version: $Id$
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 "config.h"
22#include <misc/auxiliary.h>
23
24#ifdef HAVE_PLURAL
25
26#define PLURAL_INTERNAL_DECLARATIONS
27
28#ifndef NDEBUG
29#define OUTPUT 1
30#else
31#define OUTPUT 0
32#endif
33
34#include <reporter/reporter.h>
35
36#include <coeffs/numbers.h>
37#include "coeffrings.h"
38
39#include "nc/ncSAFormula.h"
40// for CFormulaPowerMultiplier
41
42#include "monomials/ring.h"
43#include "monomials/p_polys.h"
44
45#include "nc/sca.h"
46
47
48
49
50bool ncInitSpecialPowersMultiplication(ring r)
51{
52#if OUTPUT 
53  Print("ncInitSpecialPowersMultiplication(ring), ring: \n");
54  rWrite(r);
55  PrintLn();
56#endif
57
58  assume(rIsPluralRing(r));
59  assume(!rIsSCA(r));
60
61
62  if( r->GetNC()->GetFormulaPowerMultiplier() != NULL )
63  {
64    WarnS("Already defined!");
65    return false;
66  }
67
68 
69  r->GetNC()->GetFormulaPowerMultiplier() = new CFormulaPowerMultiplier(r);
70
71  return true;
72 
73}
74
75
76
77
78
79static inline Enum_ncSAType AnalyzePairType(const ring r, int i, int j)
80{
81#if OUTPUT 
82  Print("AnalyzePair(ring, i: %d, j: %d)!", i, j);
83  PrintLn();
84#endif
85
86  const int N = r->N;
87
88  assume(i < j);
89  assume(i > 0);
90  assume(j <= N);
91
92
93  const number q = p_GetCoeff(GetC(r, i, j), r);
94  const poly d = GetD(r, i, j);
95
96#if OUTPUT 
97  Print("C_{%d, %d} = ", i, j);  { number t = n_Copy(q, r);  n_Write(t, r);  n_Delete(&t, r); };
98  Print("D_{%d, %d} = ", i, j); p_Write(d, r);
99#endif
100
101//  const number q = p_GetCoeff(c, r);
102
103  if( d == NULL)
104  {
105
106    if( n_IsOne(q, r) ) // commutative
107      return _ncSA_1xy0x0y0;
108
109    if( n_IsMOne(q, r) )
110      return _ncSA_Mxy0x0y0;
111
112    return _ncSA_Qxy0x0y0;
113  } else
114  {
115    if( n_IsOne(q, r) ) // "Lie" case
116    {
117      if( pNext(d) == NULL ) // Our Main Special Case!
118      {
119        const number g = p_GetCoeff(d, r);
120
121        if( p_LmIsConstantComp(d, r) ) // Weyl
122          return _ncSA_1xy0x0yG;
123
124        const int k = p_IsPurePower(d, r); // k if not pure power
125
126        if( k > 0 )
127          if( p_GetExp(d, k, r) == 1 )
128          {
129            if(k == i)
130              return _ncSA_1xyAx0y0;
131
132            if(k == j)
133              return _ncSA_1xy0xBy0;             
134          }
135      }
136    }
137  }
138
139  return _ncSA_notImplemented;
140}
141
142
143
144CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring r): m_BaseRing(r), m_NVars(r->N)
145{
146#if OUTPUT 
147  Print("CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring)!");
148  PrintLn();
149#endif
150
151  m_SAPairTypes = (Enum_ncSAType*)omAlloc0( ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) );
152
153  for( int i = 1; i < NVars(); i++ )
154    for( int j = i + 1; j <= NVars(); j++ )
155      GetPair(i, j) = AnalyzePairType(GetBasering(), i, j);
156}
157
158
159
160
161CFormulaPowerMultiplier::~CFormulaPowerMultiplier()
162{
163#if OUTPUT 
164  Print("CFormulaPowerMultiplier::~CFormulaPowerMultiplier()!");
165  PrintLn();
166#endif
167
168  omFreeSize((ADDRESS)m_SAPairTypes, ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) );
169}
170
171
172
173
174///////////////////////////////////////////////////////////////////////////////////////////
175static inline void CorrectPolyWRTOrdering(poly &pResult, const ring r)
176{
177  if( pNext(pResult) != NULL )
178  {
179    const int cmp = p_LmCmp(pResult, pNext(pResult), r);
180    assume( cmp != 0 ); // must not be equal!!!
181    if( cmp != 1 ) // Wrong order!!!
182      pResult = pReverse(pResult); // Reverse!!!
183  }
184  p_Test(pResult, r);
185}
186///////////////////////////////////////////////////////////////////////////////////////////
187static inline poly ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r)
188{
189#if OUTPUT 
190  Print("ncSA_1xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); 
191  PrintLn();
192#endif
193
194  poly p = p_One( r);
195  p_SetExp(p, j, m, r);
196  p_SetExp(p, i, n, r);
197  p_Setm(p, r);
198
199  p_Test(p, r);
200
201  return p;
202
203} //    return ncSA_1xy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
204///////////////////////////////////////////////////////////////////////////////////////////
205static inline poly ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r)
206{
207#if OUTPUT 
208  Print("ncSA_{M = -1}xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); 
209  PrintLn();
210#endif
211
212  const int  sign = 1 - ((n & (m & 1)) << 1);
213  poly p = p_ISet(sign, r);
214  p_SetExp(p, j, m, r);
215  p_SetExp(p, i, n, r);
216  p_Setm(p, r);
217
218
219  p_Test(p, r);
220
221  return p;
222
223} //    return ncSA_Mxy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
224///////////////////////////////////////////////////////////////////////////////////////////
225static inline poly ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r)
226{
227#if OUTPUT 
228  Print("ncSA_Qxy0x0y0(var(%d)^{%d}, var(%d)^{%d}, Q, r)!", j, m, i, n); 
229  PrintLn();
230#endif
231
232  int min, max;
233
234  if( n < m )
235  {
236    min = n;
237    max = m;
238  } else
239  {
240    min = m;
241    max = n;
242  }
243
244  number qN;
245
246  if( max == 1 )
247    qN = n_Copy(m_q, r);
248  else
249  {
250    number t;
251    n_Power(m_q, max, &t, r);
252
253    if( min > 1 )
254    {
255      n_Power(t, min, &qN, r);
256      n_Delete(&t, r);     
257    }
258    else
259      qN = t;
260  }
261
262  poly p = p_NSet(qN, r);
263  p_SetExp(p, j, m, r);
264  p_SetExp(p, i, n, r);
265  p_Setm(p, r);
266
267
268  p_Test(p, r);
269
270  return p;
271
272} //    return ncSA_Qxy0x0y0(GetI(), GetJ(), expRight, expLeft, m_q, r);
273///////////////////////////////////////////////////////////////////////////////////////////
274static inline poly ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r)
275{
276#if OUTPUT 
277  Print("ncSA_1xy0x0yG(var(%d)^{%d}, var(%d)^{%d}, G, r)!", j, m, i, n); 
278  PrintLn();
279  number t = n_Copy(m_g, r);
280  PrintS("Parameter G: "); n_Write(t, r);
281  n_Delete(&t, r);
282#endif
283
284  int kn = n;
285  int km = m;
286
287  number c = n_Init(1, r);
288
289  poly p = p_One( r);
290
291  p_SetExp(p, j, km--, r); // y ^ (m-k)
292  p_SetExp(p, i, kn--, r); // x ^ (n-k)
293
294  p_Setm(p, r); // pResult = x^n * y^m
295
296
297  poly pResult = p;
298  poly pLast = p;
299
300  int min = si_min(m, n);
301
302  int k = 1;
303
304  for(; k < min; k++ )
305  {
306    number t = n_Init(km + 1, r);
307    n_InpMult(t, m_g, r); // t = ((m - k) + 1) * gamma
308    n_InpMult(c, t, r);   // c = c'* ((m - k) + 1) * gamma
309    n_Delete(&t, r);
310
311    t = n_Init(kn + 1, r);
312    n_InpMult(c, t, r);   // c = (c'* ((m - k) + 1) * gamma) * ((n - k) + 1)
313    n_Delete(&t, r);
314
315    t = n_Init(k, r);
316    c = n_Div(c, t, r);
317    n_Delete(&t, r);
318
319//    n_Normalize(c, r);
320
321    t = n_Copy(c, r); // not the last!
322
323    p = p_NSet(t, r);
324
325    p_SetExp(p, j, km--, r); // y ^ (m-k)
326    p_SetExp(p, i, kn--, r); // x ^ (n-k)
327
328    p_Setm(p, r); // pResult = x^n * y^m
329
330    pNext(pLast) = p;
331    pLast = p;
332  }
333
334  assume(k == min);
335  assume((km == 0) || (kn == 0) );
336
337  {
338    n_InpMult(c, m_g, r);   // c = c'* gamma
339
340    if( km > 0 )
341    {
342      number t = n_Init(km + 1, r);
343      n_InpMult(c, t, r);   // c = (c'* gamma) * (m - k + 1)
344      n_Delete(&t, r);
345    }
346
347    if( kn > 0 )
348    {
349      number t = n_Init(kn + 1, r);
350      n_InpMult(c, t, r);   // c = (c'* gamma) * (n - k + 1)
351      n_Delete(&t, r);
352    }
353
354    number t = n_Init(k, r); // c = ((c'* gamma) * ((n - k + 1) * (m - k + 1))) / k;
355    c = n_Div(c, t, r);
356    n_Delete(&t, r);
357  }
358
359  p = p_NSet(c, r);
360
361  p_SetExp(p, j, km, r); // y ^ (m-k)
362  p_SetExp(p, i, kn, r); // x ^ (n-k)
363
364  p_Setm(p, r); //
365
366  pNext(pLast) = p;
367
368  CorrectPolyWRTOrdering(pResult, r);
369
370  return pResult;
371}
372///////////////////////////////////////////////////////////////////////////////////////////
373///////////////////////////////////////////////////////////////////////////////////////////
374static inline poly ncSA_ShiftAx(int i, int j, int n, int m, const number m_shiftCoef, const ring r)
375{
376  // Char == 0, otherwise - problem!
377
378  int k = m; // to 0
379
380  number c = n_Init(1, r); // k = m, C_k = 1
381  poly p = p_One( r);
382
383  p_SetExp(p, j, k, r); // Y^{k}
384  p_SetExp(p, i, n, r); 
385
386  p_Setm(p, r); // pResult = C_k * x^n * y^k, k == m
387
388
389  poly pResult = p;
390  poly pLast = p;
391
392  number nn =  n_Init(n, r); // number(n)!
393  n_InpMult(nn, m_shiftCoef, r); // nn = (alpha*n)
394
395  --k;
396
397  int mk = 1; // mk = (m - k)
398
399  for(; k > 0; k-- )
400  {
401    number t = n_Init(k + 1, r);  // t = k+1
402    n_InpMult(c, t, r);           // c = c' * (k+1)
403    n_InpMult(c, nn, r);          // c = (c' * (k+1)) * (alpha * n)
404
405    n_Delete(&t, r);
406    t = n_Init(mk++, r);         
407    c = n_Div(c, t, r);           // c = ((c' * (k+1))  * (alpha * n)) / (m-k);
408    n_Delete(&t, r);
409
410//    n_Normalize(c, r);
411
412    t = n_Copy(c, r); // not the last!
413
414    p = p_NSet(t, r);
415
416    p_SetExp(p, j, k, r); // y^k
417    p_SetExp(p, i, n, r); // x^n
418
419    p_Setm(p, r); // pResult = x^n * y^m
420
421    pNext(pLast) = p;
422    pLast = p;
423  }
424
425  assume(k == 0);
426
427  {
428    n_InpMult(c, nn, r);          // c = (c' * (0+1)) * (alpha * n)
429
430    number t = n_Init(m, r);         
431    c = n_Div(c, t, r);           // c = ((c' * (0+1))  * (alpha * n)) / (m-0);
432    n_Delete(&t, r);
433  }
434
435  n_Delete(&nn, r);
436
437  p = p_NSet(c, r);
438
439  p_SetExp(p, j, k, r); // y^k
440  p_SetExp(p, i, n, r); // x^n
441
442  p_Setm(p, r); //
443
444  pNext(pLast) = p;
445
446  CorrectPolyWRTOrdering(pResult, r);
447
448  return pResult;
449}
450///////////////////////////////////////////////////////////////////////////////////////////
451static inline poly ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
452{
453#if OUTPUT 
454  Print("ncSA_1xyAx0y0(var(%d)^{%d}, var(%d)^{%d}, A, r)!", j, m, i, n); 
455  PrintLn();
456  number t = n_Copy(m_shiftCoef, r);
457  PrintS("Parameter A: "); n_Write(t, r);
458  n_Delete(&t, r); 
459#endif
460
461  return ncSA_ShiftAx(i, j, n, m, m_shiftCoef, r);
462}
463///////////////////////////////////////////////////////////////////////////////////////////
464static inline poly ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
465{
466#if OUTPUT 
467  Print("ncSA_1xy0xBy0(var(%d)^{%d}, var(%d)^{%d}, B, r)!", j, m, i, n); 
468  PrintLn();
469  number t = n_Copy(m_shiftCoef, r);
470  PrintS("Parameter B: "); n_Write(t, r);
471  n_Delete(&t, r); 
472#endif
473
474  return ncSA_ShiftAx(j, i, m, n, m_shiftCoef, r);
475}
476///////////////////////////////////////////////////////////////////////////////////////////
477///////////////////////////////////////////////////////////////////////////////////////////
478///////////////////////////////////////////////////////////////////////////////////////////
479///////////////////////////////////////////////////////////////////////////////////////////
480
481
482static inline poly ncSA_Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r)
483{
484#if OUTPUT 
485  Print("ncSA_Multiply(type: %d, ring, (var(%d)^{%d} * var(%d)^{%d}, r)!", (int)type, j, m, i, n); 
486  PrintLn();
487#endif
488
489  assume( type != _ncSA_notImplemented );
490  assume( (n > 0) && (m > 0) );
491
492  if( type == _ncSA_1xy0x0y0 )
493    return ::ncSA_1xy0x0y0(i, j, n, m, r);
494
495  if( type == _ncSA_Mxy0x0y0 )
496    return ::ncSA_Mxy0x0y0(i, j, n, m, r);
497
498  if( type == _ncSA_Qxy0x0y0 )
499  {
500    const number q = p_GetCoeff(GetC(r, i, j), r);
501    return ::ncSA_Qxy0x0y0(i, j, n, m, q, r);
502  }
503
504  const number g = p_GetCoeff(GetD(r, i, j), r);
505
506  if( type == _ncSA_1xy0x0yG ) // Weyl
507    return ::ncSA_1xy0x0yG(i, j, n, m, g, r);
508
509  if( type == _ncSA_1xyAx0y0 ) // Shift 1
510    return ::ncSA_1xyAx0y0(i, j, n, m, g, r);
511
512  if( type == _ncSA_1xy0xBy0 ) // Shift 2
513    return ::ncSA_1xy0xBy0(i, j, n, m, g, r);
514
515  assume( type == _ncSA_notImplemented );
516
517  return NULL;
518}
519
520
521poly CFormulaPowerMultiplier::Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r)
522{
523  return ncSA_Multiply( type, i, j, n, m, r);
524}
525
526
527Enum_ncSAType CFormulaPowerMultiplier::AnalyzePair(const ring r, int i, int j)
528{
529  return ::AnalyzePairType( r, i, j);
530}
531
532poly CFormulaPowerMultiplier::Multiply( int i, int j, const int n, const int m)
533{
534  return ncSA_Multiply( GetPair(i, j), i, j, n, m, GetBasering());
535}
536
537
538
539
540poly CFormulaPowerMultiplier::ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r)
541{
542  return ::ncSA_1xy0x0y0(i, j, n, m, r);
543}
544
545poly CFormulaPowerMultiplier::ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r)
546{
547  return ::ncSA_Mxy0x0y0(i, j, n, m, r);
548}
549
550poly CFormulaPowerMultiplier::ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r)
551{
552  return ::ncSA_Qxy0x0y0(i, j, n, m, m_q, r);
553}
554
555poly CFormulaPowerMultiplier::ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r)
556{
557  return ::ncSA_1xy0x0yG(i, j, n, m, m_g, r);
558}
559
560poly CFormulaPowerMultiplier::ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
561{
562  return ::ncSA_1xyAx0y0(i, j, n, m, m_shiftCoef, r);
563}
564
565poly CFormulaPowerMultiplier::ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
566{
567  return ::ncSA_1xy0xBy0(i, j, n, m, m_shiftCoef, r);
568}
569#endif
Note: See TracBrowser for help on using the repository browser.