source: git/libpolys/polys/nc/ncSAFormula.cc @ d2f720

spielwiese
Last change on this file since d2f720 was d2f720, checked in by Niels Ranosch <ranosch@…>, 13 years ago
CHG: unused variables
  • 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); // not used for now
120         if( p_LmIsConstantComp(d, r) ) // Weyl
121          return _ncSA_1xy0x0yG;
122
123        const int k = p_IsPurePower(d, r); // k if not pure power
124
125        if( k > 0 )
126          if( p_GetExp(d, k, r) == 1 )
127          {
128            if(k == i)
129              return _ncSA_1xyAx0y0;
130
131            if(k == j)
132              return _ncSA_1xy0xBy0;             
133          }
134      }
135    }
136  }
137
138  return _ncSA_notImplemented;
139}
140
141
142
143CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring r): m_BaseRing(r), m_NVars(r->N)
144{
145#if OUTPUT 
146  Print("CFormulaPowerMultiplier::CFormulaPowerMultiplier(ring)!");
147  PrintLn();
148#endif
149
150  m_SAPairTypes = (Enum_ncSAType*)omAlloc0( ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) );
151
152  for( int i = 1; i < NVars(); i++ )
153    for( int j = i + 1; j <= NVars(); j++ )
154      GetPair(i, j) = AnalyzePairType(GetBasering(), i, j);
155}
156
157
158
159
160CFormulaPowerMultiplier::~CFormulaPowerMultiplier()
161{
162#if OUTPUT 
163  Print("CFormulaPowerMultiplier::~CFormulaPowerMultiplier()!");
164  PrintLn();
165#endif
166
167  omFreeSize((ADDRESS)m_SAPairTypes, ((NVars() * (NVars()-1)) / 2) * sizeof(Enum_ncSAType) );
168}
169
170
171
172
173///////////////////////////////////////////////////////////////////////////////////////////
174static inline void CorrectPolyWRTOrdering(poly &pResult, const ring r)
175{
176  if( pNext(pResult) != NULL )
177  {
178    const int cmp = p_LmCmp(pResult, pNext(pResult), r);
179    assume( cmp != 0 ); // must not be equal!!!
180    if( cmp != 1 ) // Wrong order!!!
181      pResult = pReverse(pResult); // Reverse!!!
182  }
183  p_Test(pResult, r);
184}
185///////////////////////////////////////////////////////////////////////////////////////////
186static inline poly ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r)
187{
188#if OUTPUT 
189  Print("ncSA_1xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); 
190  PrintLn();
191#endif
192
193  poly p = p_One( r);
194  p_SetExp(p, j, m, r);
195  p_SetExp(p, i, n, r);
196  p_Setm(p, r);
197
198  p_Test(p, r);
199
200  return p;
201
202} //    return ncSA_1xy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
203///////////////////////////////////////////////////////////////////////////////////////////
204static inline poly ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r)
205{
206#if OUTPUT 
207  Print("ncSA_{M = -1}xy0x0y0(var(%d)^{%d}, var(%d)^{%d}, r)!", j, m, i, n); 
208  PrintLn();
209#endif
210
211  const int  sign = 1 - ((n & (m & 1)) << 1);
212  poly p = p_ISet(sign, r);
213  p_SetExp(p, j, m, r);
214  p_SetExp(p, i, n, r);
215  p_Setm(p, r);
216
217
218  p_Test(p, r);
219
220  return p;
221
222} //    return ncSA_Mxy0x0y0(GetI(), GetJ(), expRight, expLeft, r);
223///////////////////////////////////////////////////////////////////////////////////////////
224static inline poly ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r)
225{
226#if OUTPUT 
227  Print("ncSA_Qxy0x0y0(var(%d)^{%d}, var(%d)^{%d}, Q, r)!", j, m, i, n); 
228  PrintLn();
229#endif
230
231  int min, max;
232
233  if( n < m )
234  {
235    min = n;
236    max = m;
237  } else
238  {
239    min = m;
240    max = n;
241  }
242
243  number qN;
244
245  if( max == 1 )
246    qN = n_Copy(m_q, r);
247  else
248  {
249    number t;
250    n_Power(m_q, max, &t, r);
251
252    if( min > 1 )
253    {
254      n_Power(t, min, &qN, r);
255      n_Delete(&t, r);     
256    }
257    else
258      qN = t;
259  }
260
261  poly p = p_NSet(qN, r);
262  p_SetExp(p, j, m, r);
263  p_SetExp(p, i, n, r);
264  p_Setm(p, r);
265
266
267  p_Test(p, r);
268
269  return p;
270
271} //    return ncSA_Qxy0x0y0(GetI(), GetJ(), expRight, expLeft, m_q, r);
272///////////////////////////////////////////////////////////////////////////////////////////
273static inline poly ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r)
274{
275#if OUTPUT 
276  Print("ncSA_1xy0x0yG(var(%d)^{%d}, var(%d)^{%d}, G, r)!", j, m, i, n); 
277  PrintLn();
278  number t = n_Copy(m_g, r);
279  PrintS("Parameter G: "); n_Write(t, r);
280  n_Delete(&t, r);
281#endif
282
283  int kn = n;
284  int km = m;
285
286  number c = n_Init(1, r);
287
288  poly p = p_One( r);
289
290  p_SetExp(p, j, km--, r); // y ^ (m-k)
291  p_SetExp(p, i, kn--, r); // x ^ (n-k)
292
293  p_Setm(p, r); // pResult = x^n * y^m
294
295
296  poly pResult = p;
297  poly pLast = p;
298
299  int min = si_min(m, n);
300
301  int k = 1;
302
303  for(; k < min; k++ )
304  {
305    number t = n_Init(km + 1, r);
306    n_InpMult(t, m_g, r); // t = ((m - k) + 1) * gamma
307    n_InpMult(c, t, r);   // c = c'* ((m - k) + 1) * gamma
308    n_Delete(&t, r);
309
310    t = n_Init(kn + 1, r);
311    n_InpMult(c, t, r);   // c = (c'* ((m - k) + 1) * gamma) * ((n - k) + 1)
312    n_Delete(&t, r);
313
314    t = n_Init(k, r);
315    c = n_Div(c, t, r);
316    n_Delete(&t, r);
317
318//    n_Normalize(c, r);
319
320    t = n_Copy(c, r); // not the last!
321
322    p = p_NSet(t, r);
323
324    p_SetExp(p, j, km--, r); // y ^ (m-k)
325    p_SetExp(p, i, kn--, r); // x ^ (n-k)
326
327    p_Setm(p, r); // pResult = x^n * y^m
328
329    pNext(pLast) = p;
330    pLast = p;
331  }
332
333  assume(k == min);
334  assume((km == 0) || (kn == 0) );
335
336  {
337    n_InpMult(c, m_g, r);   // c = c'* gamma
338
339    if( km > 0 )
340    {
341      number t = n_Init(km + 1, r);
342      n_InpMult(c, t, r);   // c = (c'* gamma) * (m - k + 1)
343      n_Delete(&t, r);
344    }
345
346    if( kn > 0 )
347    {
348      number t = n_Init(kn + 1, r);
349      n_InpMult(c, t, r);   // c = (c'* gamma) * (n - k + 1)
350      n_Delete(&t, r);
351    }
352
353    number t = n_Init(k, r); // c = ((c'* gamma) * ((n - k + 1) * (m - k + 1))) / k;
354    c = n_Div(c, t, r);
355    n_Delete(&t, r);
356  }
357
358  p = p_NSet(c, r);
359
360  p_SetExp(p, j, km, r); // y ^ (m-k)
361  p_SetExp(p, i, kn, r); // x ^ (n-k)
362
363  p_Setm(p, r); //
364
365  pNext(pLast) = p;
366
367  CorrectPolyWRTOrdering(pResult, r);
368
369  return pResult;
370}
371///////////////////////////////////////////////////////////////////////////////////////////
372///////////////////////////////////////////////////////////////////////////////////////////
373static inline poly ncSA_ShiftAx(int i, int j, int n, int m, const number m_shiftCoef, const ring r)
374{
375  // Char == 0, otherwise - problem!
376
377  int k = m; // to 0
378
379  number c = n_Init(1, r); // k = m, C_k = 1
380  poly p = p_One( r);
381
382  p_SetExp(p, j, k, r); // Y^{k}
383  p_SetExp(p, i, n, r); 
384
385  p_Setm(p, r); // pResult = C_k * x^n * y^k, k == m
386
387
388  poly pResult = p;
389  poly pLast = p;
390
391  number nn =  n_Init(n, r); // number(n)!
392  n_InpMult(nn, m_shiftCoef, r); // nn = (alpha*n)
393
394  --k;
395
396  int mk = 1; // mk = (m - k)
397
398  for(; k > 0; k-- )
399  {
400    number t = n_Init(k + 1, r);  // t = k+1
401    n_InpMult(c, t, r);           // c = c' * (k+1)
402    n_InpMult(c, nn, r);          // c = (c' * (k+1)) * (alpha * n)
403
404    n_Delete(&t, r);
405    t = n_Init(mk++, r);         
406    c = n_Div(c, t, r);           // c = ((c' * (k+1))  * (alpha * n)) / (m-k);
407    n_Delete(&t, r);
408
409//    n_Normalize(c, r);
410
411    t = n_Copy(c, r); // not the last!
412
413    p = p_NSet(t, r);
414
415    p_SetExp(p, j, k, r); // y^k
416    p_SetExp(p, i, n, r); // x^n
417
418    p_Setm(p, r); // pResult = x^n * y^m
419
420    pNext(pLast) = p;
421    pLast = p;
422  }
423
424  assume(k == 0);
425
426  {
427    n_InpMult(c, nn, r);          // c = (c' * (0+1)) * (alpha * n)
428
429    number t = n_Init(m, r);         
430    c = n_Div(c, t, r);           // c = ((c' * (0+1))  * (alpha * n)) / (m-0);
431    n_Delete(&t, r);
432  }
433
434  n_Delete(&nn, r);
435
436  p = p_NSet(c, r);
437
438  p_SetExp(p, j, k, r); // y^k
439  p_SetExp(p, i, n, r); // x^n
440
441  p_Setm(p, r); //
442
443  pNext(pLast) = p;
444
445  CorrectPolyWRTOrdering(pResult, r);
446
447  return pResult;
448}
449///////////////////////////////////////////////////////////////////////////////////////////
450static inline poly ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
451{
452#if OUTPUT 
453  Print("ncSA_1xyAx0y0(var(%d)^{%d}, var(%d)^{%d}, A, r)!", j, m, i, n); 
454  PrintLn();
455  number t = n_Copy(m_shiftCoef, r);
456  PrintS("Parameter A: "); n_Write(t, r);
457  n_Delete(&t, r); 
458#endif
459
460  return ncSA_ShiftAx(i, j, n, m, m_shiftCoef, r);
461}
462///////////////////////////////////////////////////////////////////////////////////////////
463static inline poly ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
464{
465#if OUTPUT 
466  Print("ncSA_1xy0xBy0(var(%d)^{%d}, var(%d)^{%d}, B, r)!", j, m, i, n); 
467  PrintLn();
468  number t = n_Copy(m_shiftCoef, r);
469  PrintS("Parameter B: "); n_Write(t, r);
470  n_Delete(&t, r); 
471#endif
472
473  return ncSA_ShiftAx(j, i, m, n, m_shiftCoef, r);
474}
475///////////////////////////////////////////////////////////////////////////////////////////
476///////////////////////////////////////////////////////////////////////////////////////////
477///////////////////////////////////////////////////////////////////////////////////////////
478///////////////////////////////////////////////////////////////////////////////////////////
479
480
481static inline poly ncSA_Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r)
482{
483#if OUTPUT 
484  Print("ncSA_Multiply(type: %d, ring, (var(%d)^{%d} * var(%d)^{%d}, r)!", (int)type, j, m, i, n); 
485  PrintLn();
486#endif
487
488  assume( type != _ncSA_notImplemented );
489  assume( (n > 0) && (m > 0) );
490
491  if( type == _ncSA_1xy0x0y0 )
492    return ::ncSA_1xy0x0y0(i, j, n, m, r);
493
494  if( type == _ncSA_Mxy0x0y0 )
495    return ::ncSA_Mxy0x0y0(i, j, n, m, r);
496
497  if( type == _ncSA_Qxy0x0y0 )
498  {
499    const number q = p_GetCoeff(GetC(r, i, j), r);
500    return ::ncSA_Qxy0x0y0(i, j, n, m, q, r);
501  }
502
503  const number g = p_GetCoeff(GetD(r, i, j), r);
504
505  if( type == _ncSA_1xy0x0yG ) // Weyl
506    return ::ncSA_1xy0x0yG(i, j, n, m, g, r);
507
508  if( type == _ncSA_1xyAx0y0 ) // Shift 1
509    return ::ncSA_1xyAx0y0(i, j, n, m, g, r);
510
511  if( type == _ncSA_1xy0xBy0 ) // Shift 2
512    return ::ncSA_1xy0xBy0(i, j, n, m, g, r);
513
514  assume( type == _ncSA_notImplemented );
515
516  return NULL;
517}
518
519
520poly CFormulaPowerMultiplier::Multiply( Enum_ncSAType type, const int i, const int j, const int n, const int m, const ring r)
521{
522  return ncSA_Multiply( type, i, j, n, m, r);
523}
524
525
526Enum_ncSAType CFormulaPowerMultiplier::AnalyzePair(const ring r, int i, int j)
527{
528  return ::AnalyzePairType( r, i, j);
529}
530
531poly CFormulaPowerMultiplier::Multiply( int i, int j, const int n, const int m)
532{
533  return ncSA_Multiply( GetPair(i, j), i, j, n, m, GetBasering());
534}
535
536
537
538
539poly CFormulaPowerMultiplier::ncSA_1xy0x0y0(const int i, const int j, const int n, const int m, const ring r)
540{
541  return ::ncSA_1xy0x0y0(i, j, n, m, r);
542}
543
544poly CFormulaPowerMultiplier::ncSA_Mxy0x0y0(const int i, const int j, const int n, const int m, const ring r)
545{
546  return ::ncSA_Mxy0x0y0(i, j, n, m, r);
547}
548
549poly CFormulaPowerMultiplier::ncSA_Qxy0x0y0(const int i, const int j, const int n, const int m, const number m_q, const ring r)
550{
551  return ::ncSA_Qxy0x0y0(i, j, n, m, m_q, r);
552}
553
554poly CFormulaPowerMultiplier::ncSA_1xy0x0yG(const int i, const int j, const int n, const int m, const number m_g, const ring r)
555{
556  return ::ncSA_1xy0x0yG(i, j, n, m, m_g, r);
557}
558
559poly CFormulaPowerMultiplier::ncSA_1xyAx0y0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
560{
561  return ::ncSA_1xyAx0y0(i, j, n, m, m_shiftCoef, r);
562}
563
564poly CFormulaPowerMultiplier::ncSA_1xy0xBy0(const int i, const int j, const int n, const int m, const number m_shiftCoef, const ring r)
565{
566  return ::ncSA_1xy0xBy0(i, j, n, m, m_shiftCoef, r);
567}
568#endif
Note: See TracBrowser for help on using the repository browser.