source: git/kernel/ncSAFormula.cc @ 098f98f

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