source: git/kernel/ncSAFormula.cc @ 61d32c

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