source: git/libpolys/polys/nc/sca.cc @ c7aad0

jengelh-datetimespielwiese
Last change on this file since c7aad0 was c7aad0, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
FIX: non-commutative (left) GB should now work. NOTE: SCA was disabled with this!
  • Property mode set to 100644
File size: 36.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    sca.cc
6 *  Purpose: supercommutative kernel procedures
7 *  Author:  motsak (Oleksandr Motsak)
8 *  Created: 2006/12/18
9 *  Version: $Id$
10 *******************************************************************/
11
12// set it here if needed.
13#define OUTPUT 0
14#define MYTEST 0
15
16#if MYTEST
17#define OM_CHECK 4
18#define OM_TRACK 5
19#endif
20
21// #define PDEBUG 2
22#include "config.h"
23#include <misc/auxiliary.h>
24
25#ifdef HAVE_PLURAL
26
27// for
28#define PLURAL_INTERNAL_DECLARATIONS
29
30#include <polys/nc/sca.h>
31#include <polys/nc/nc.h>
32// #include <polys/gring.h>
33
34
35#include <coeffs/numbers.h>
36#include <polys/coeffrings.h>
37
38
39// #include <polys/febase.h>
40#include <misc/options.h>
41
42#include <polys/monomials/p_polys.h>
43
44// #include <polys/kutil.h>
45#include <polys/simpleideals.h>
46#include <misc/intvec.h>
47
48#include <polys/monomials/ring.h>
49#include <polys/kbuckets.h>
50// #include <polys/kstd1.h>
51#include <polys/sbuckets.h>
52
53#include <polys/prCopy.h>
54
55#include <polys/operations/p_Mult_q.h>
56#include <polys/templates/p_MemAdd.h>
57
58// #include <polys/kutil.h>
59// #include <polys/kstd1.h>
60
61#include <polys/weight.h>
62
63
64// poly functions defined in p_Procs :
65
66// return pPoly * pMonom; preserve pPoly and pMonom.
67poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing, poly &);
68
69// return pMonom * pPoly; preserve pPoly and pMonom.
70static poly sca_mm_Mult_pp(const poly pMonom, const poly pPoly, const ring rRing);
71
72// return pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
73poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing);
74
75// return pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
76static poly sca_mm_Mult_p(const poly pMonom, poly pPoly, const ring rRing);
77
78
79// compute the spolynomial of p1 and p2
80poly sca_SPoly(const poly p1, const poly p2, const ring r);
81poly sca_ReduceSpoly(const poly p1, poly p2, const ring r);
82
83
84////////////////////////////////////////////////////////////////////////////////////////////////////
85// Super Commutative Algebra extension by Oleksandr
86////////////////////////////////////////////////////////////////////////////////////////////////////
87
88// returns the sign of: lm(pMonomM) * lm(pMonomMM),
89// preserves input, may return +/-1, 0
90static inline int sca_Sign_mm_Mult_mm( const poly pMonomM, const poly pMonomMM, const ring rRing )
91{
92#ifdef PDEBUG
93    p_Test(pMonomM,  rRing);
94    p_Test(pMonomMM, rRing);
95#endif
96
97    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
98    const unsigned int iLastAltVar  = scaLastAltVar(rRing);
99
100    register unsigned int tpower = 0;
101    register unsigned int cpower = 0;
102
103    for( register unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
104    {
105      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
106      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
107
108#ifdef PDEBUG
109      assume( iExpM <= 1);
110      assume( iExpMM <= 1);
111#endif
112
113      if( iExpMM != 0 ) // TODO: think about eliminating there if-s...
114      {
115        if( iExpM != 0 )
116        {
117          return 0; // lm(pMonomM) * lm(pMonomMM) == 0
118        }
119        tpower ^= cpower; // compute degree of (-1).
120      }
121      cpower ^= iExpM;
122    }
123
124#ifdef PDEBUG
125    assume(tpower <= 1);
126#endif
127
128    // 1 => -1  // degree is odd => negate coeff.
129    // 0 =>  1
130
131    return(1 - (tpower << 1) );
132}
133
134
135
136
137// returns and changes pMonomM: lt(pMonomM) = lt(pMonomM) * lt(pMonomMM),
138// preserves pMonomMM. may return NULL!
139// if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomM) * lt(pMonomMM)
140// if result == NULL => pMonomM MUST BE deleted manually!
141static inline poly sca_m_Mult_mm( poly pMonomM, const poly pMonomMM, const ring rRing )
142{
143#ifdef PDEBUG
144    p_Test(pMonomM,  rRing);
145    p_Test(pMonomMM, rRing);
146#endif
147
148    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
149    const unsigned int iLastAltVar = scaLastAltVar(rRing);
150
151    register unsigned int tpower = 0;
152    register unsigned int cpower = 0;
153
154    for( register unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
155    {
156      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
157      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
158
159#ifdef PDEBUG
160      assume( iExpM <= 1);
161      assume( iExpMM <= 1);
162#endif
163
164      if( iExpMM != 0 )
165      {
166        if( iExpM != 0 ) // result is zero!
167        {
168          return NULL; // we do nothing with pMonomM in this case!
169        }
170
171        tpower ^= cpower; // compute degree of (-1).
172      }
173
174      cpower ^= iExpM;
175    }
176
177#ifdef PDEBUG
178    assume(tpower <= 1);
179#endif
180
181    p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
182
183    number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
184
185    if( (tpower) != 0 ) // degree is odd => negate coeff.
186      nCoeffM = n_Neg(nCoeffM, rRing); // negate nCoeff (will destroy the original number)
187
188    const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
189
190    number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing); // new number!
191
192    p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
193
194#ifdef PDEBUG
195    p_LmTest(pMonomM, rRing);
196#endif
197
198    return(pMonomM);
199}
200
201// returns and changes pMonomM: lt(pMonomM) = lt(pMonomMM) * lt(pMonomM),
202// preserves pMonomMM. may return NULL!
203// if result != NULL => next(result) = next(pMonomM), lt(result) = lt(pMonomMM) * lt(pMonomM)
204// if result == NULL => pMonomM MUST BE deleted manually!
205static inline poly sca_mm_Mult_m( const poly pMonomMM, poly pMonomM, const ring rRing )
206{
207#ifdef PDEBUG
208    p_Test(pMonomM,  rRing);
209    p_Test(pMonomMM, rRing);
210#endif
211
212    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
213    const unsigned int iLastAltVar = scaLastAltVar(rRing);
214
215    register unsigned int tpower = 0;
216    register unsigned int cpower = 0;
217
218    for( register unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
219    {
220      const unsigned int iExpMM = p_GetExp(pMonomMM, j, rRing);
221      const unsigned int iExpM  = p_GetExp(pMonomM,  j, rRing);
222
223#ifdef PDEBUG
224      assume( iExpM <= 1);
225      assume( iExpMM <= 1);
226#endif
227
228      if( iExpM != 0 )
229      {
230        if( iExpMM != 0 ) // result is zero!
231        {
232          return NULL; // we do nothing with pMonomM in this case!
233        }
234
235        tpower ^= cpower; // compute degree of (-1).
236      }
237
238      cpower ^= iExpMM;
239    }
240
241#ifdef PDEBUG
242    assume(tpower <= 1);
243#endif
244
245    p_ExpVectorAdd(pMonomM, pMonomMM, rRing); // "exponents" are additive!!!
246
247    number nCoeffM = p_GetCoeff(pMonomM, rRing); // no new copy! should be deleted!
248
249    if( (tpower) != 0 ) // degree is odd => negate coeff.
250      nCoeffM = n_Neg(nCoeffM, rRing); // negate nCoeff (will destroy the original number), creates new number!
251
252    const number nCoeffMM = p_GetCoeff(pMonomMM, rRing); // no new copy!
253
254    number nCoeff = n_Mult(nCoeffM, nCoeffMM, rRing); // new number!
255
256    p_SetCoeff(pMonomM, nCoeff, rRing); // delete lc(pMonomM) and set lc(pMonomM) = nCoeff
257
258#ifdef PDEBUG
259    p_LmTest(pMonomM, rRing);
260#endif
261
262    return(pMonomM);
263}
264
265
266
267// returns: result = lt(pMonom1) * lt(pMonom2),
268// preserves pMonom1, pMonom2. may return NULL!
269// if result != NULL => next(result) = NULL, lt(result) = lt(pMonom1) * lt(pMonom2)
270static inline poly sca_mm_Mult_mm( poly pMonom1, const poly pMonom2, const ring rRing )
271{
272#ifdef PDEBUG
273    p_Test(pMonom1, rRing);
274    p_Test(pMonom2, rRing);
275#endif
276
277    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
278    const unsigned int iLastAltVar = scaLastAltVar(rRing);
279
280    register unsigned int tpower = 0;
281    register unsigned int cpower = 0;
282
283    for( register unsigned int j = iLastAltVar; j >= iFirstAltVar; j-- )
284    {
285      const unsigned int iExp1 = p_GetExp(pMonom1, j, rRing);
286      const unsigned int iExp2 = p_GetExp(pMonom2, j, rRing);
287
288#ifdef PDEBUG
289      assume( iExp1 <= 1);
290      assume( iExp2 <= 1);
291#endif
292
293      if( iExp2 != 0 )
294      {
295        if( iExp1 != 0 ) // result is zero!
296        {
297          return NULL;
298        }
299        tpower ^= cpower; // compute degree of (-1).
300      }
301      cpower ^= iExp1;
302    }
303
304#ifdef PDEBUG
305    assume(cpower <= 1);
306#endif
307
308    poly pResult;
309    p_AllocBin(pResult,rRing->PolyBin,rRing);
310
311    p_ExpVectorSum(pResult, pMonom1, pMonom2, rRing); // "exponents" are additive!!!
312
313    pNext(pResult) = NULL;
314
315    const number nCoeff1 = p_GetCoeff(pMonom1, rRing); // no new copy!
316    const number nCoeff2 = p_GetCoeff(pMonom2, rRing); // no new copy!
317
318    number nCoeff = n_Mult(nCoeff1, nCoeff2, rRing); // new number!
319
320    if( (tpower) != 0 ) // degree is odd => negate coeff.
321      nCoeff = n_Neg(nCoeff, rRing); // negate nCoeff (will destroy the original number)
322
323    p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
324
325#ifdef PDEBUG
326    p_LmTest(pResult, rRing);
327#endif
328
329    return(pResult);
330}
331
332// returns: result =  x_i * lt(pMonom),
333// preserves pMonom. may return NULL!
334// if result != NULL => next(result) = NULL, lt(result) = x_i * lt(pMonom)
335static inline poly sca_xi_Mult_mm(unsigned int i, const poly pMonom, const ring rRing)
336{
337#ifdef PDEBUG
338    p_Test(pMonom, rRing);
339#endif
340
341    assume( i <= scaLastAltVar(rRing));
342    assume( scaFirstAltVar(rRing) <= i );
343
344    if( p_GetExp(pMonom, i, rRing) != 0 ) // => result is zero!
345      return NULL;
346
347    const unsigned int iFirstAltVar = scaFirstAltVar(rRing);
348
349    register unsigned int cpower = 0;
350
351    for( register unsigned int j = iFirstAltVar; j < i ; j++ )
352      cpower ^= p_GetExp(pMonom, j, rRing);
353
354#ifdef PDEBUG
355    assume(cpower <= 1);
356#endif
357
358    poly pResult = p_LmInit(pMonom, rRing);
359
360    p_SetExp(pResult, i, 1, rRing); // pResult*=X_i &&
361    p_Setm(pResult, rRing);         // addjust degree after previous step!
362
363    number nCoeff = n_Copy(p_GetCoeff(pMonom, rRing), rRing); // new number!
364
365    if( cpower != 0 ) // degree is odd => negate coeff.
366      nCoeff = n_Neg(nCoeff, rRing); // negate nCoeff (will destroy the original number)
367
368    p_SetCoeff0(pResult, nCoeff, rRing); // set lc(pResult) = nCoeff, no destruction!
369
370#ifdef PDEBUG
371    p_LmTest(pResult, rRing);
372#endif
373
374    return(pResult);
375}
376
377//-----------------------------------------------------------------------------------//
378
379// return poly = pPoly * pMonom; preserve pMonom, destroy or reuse pPoly.
380poly sca_p_Mult_mm(poly pPoly, const poly pMonom, const ring rRing)
381{
382  assume( rIsSCA(rRing) );
383
384#ifdef PDEBUG
385//  Print("sca_p_Mult_mm\n"); // !
386
387  p_Test(pPoly, rRing);
388  p_Test(pMonom, rRing);
389#endif
390
391  if( pPoly == NULL )
392    return NULL;
393
394  assume(pMonom !=NULL);
395  //if( pMonom == NULL )
396  //{
397  //  // pPoly != NULL =>
398  //  p_Delete( &pPoly, rRing );
399  //  return NULL;
400  //}
401
402  const int iComponentMonomM = p_GetComp(pMonom, rRing);
403
404  poly p = pPoly; poly* ppPrev = &pPoly;
405
406  loop
407  {
408#ifdef PDEBUG
409    p_Test(p, rRing);
410#endif
411    const int iComponent = p_GetComp(p, rRing);
412
413    if( iComponent!=0 )
414    {
415      if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
416      {
417        // REPORT_ERROR
418        Werror("sca_p_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
419        // what should we do further?!?
420
421        p_Delete( &pPoly, rRing); // delete the result AND rest
422        return NULL;
423      }
424#ifdef PDEBUG
425      if(iComponentMonomM==0 )
426      {
427        dReportError("sca_p_Mult_mm: Multiplication in the left module from the right");
428      }
429#endif
430    }
431
432    // terms will be in the same order because of quasi-ordering!
433    poly v = sca_m_Mult_mm(p, pMonom, rRing);
434
435    if( v != NULL )
436    {
437      ppPrev = &pNext(p); // fixed!
438
439    // *p is changed if v != NULL ( p == v )
440      pIter(p);
441
442      if( p == NULL )
443        break;
444    }
445    else
446    { // Upps! Zero!!! we must kill this term!
447
448      //
449      p = p_LmDeleteAndNext(p, rRing);
450
451      *ppPrev = p;
452
453      if( p == NULL )
454        break;
455    }
456  }
457
458#ifdef PDEBUG
459  p_Test(pPoly,rRing);
460#endif
461
462  return(pPoly);
463}
464
465// return new poly = pPoly * pMonom; preserve pPoly and pMonom.
466poly sca_pp_Mult_mm(const poly pPoly, const poly pMonom, const ring rRing, poly &)
467{
468  assume( rIsSCA(rRing) );
469
470#ifdef PDEBUG
471//  Print("sca_pp_Mult_mm\n"); // !
472
473  p_Test(pPoly, rRing);
474  p_Test(pMonom, rRing);
475#endif
476
477  if( ( pPoly == NULL ) /*|| ( pMonom == NULL )*/ )
478    return NULL;
479
480  const int iComponentMonomM = p_GetComp(pMonom, rRing);
481
482  poly pResult = NULL;
483  poly* ppPrev = &pResult;
484
485  for( poly p = pPoly; p!= NULL; pIter(p) )
486  {
487#ifdef PDEBUG
488    p_Test(p, rRing);
489#endif
490    const int iComponent = p_GetComp(p, rRing);
491
492    if( iComponent!=0 )
493    {
494      if( iComponentMonomM!=0 ) // TODO: make global if on iComponentMonomM =?= 0
495      {
496        // REPORT_ERROR
497        Werror("sca_pp_Mult_mm: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
498        // what should we do further?!?
499
500        p_Delete( &pResult, rRing); // delete the result
501        return NULL;
502      }
503
504#ifdef PDEBUG
505      if(iComponentMonomM==0 )
506      {
507        dReportError("sca_pp_Mult_mm: Multiplication in the left module from the right");
508      }
509#endif
510    }
511
512    // terms will be in the same order because of quasi-ordering!
513    poly v = sca_mm_Mult_mm(p, pMonom, rRing);
514
515    if( v != NULL )
516    {
517      *ppPrev = v;
518      ppPrev = &pNext(v);
519    }
520  }
521
522#ifdef PDEBUG
523  p_Test(pResult,rRing);
524#endif
525
526  return(pResult);
527}
528
529//-----------------------------------------------------------------------------------//
530
531// return x_i * pPoly; preserve pPoly.
532static inline poly sca_xi_Mult_pp(unsigned int i, const poly pPoly, const ring rRing)
533{
534  assume( rIsSCA(rRing) );
535
536#ifdef PDEBUG
537  p_Test(pPoly, rRing);
538#endif
539
540  assume(i <= scaLastAltVar(rRing));
541  assume(scaFirstAltVar(rRing) <= i);
542
543  if( pPoly == NULL )
544    return NULL;
545
546  poly pResult = NULL;
547  poly* ppPrev = &pResult;
548
549  for( poly p = pPoly; p!= NULL; pIter(p) )
550  {
551
552    // terms will be in the same order because of quasi-ordering!
553    poly v = sca_xi_Mult_mm(i, p, rRing);
554
555#ifdef PDEBUG
556    p_Test(v, rRing);
557#endif
558
559    if( v != NULL )
560    {
561      *ppPrev = v;
562      ppPrev = &pNext(*ppPrev);
563    }
564  }
565
566
567#ifdef PDEBUG
568  p_Test(pResult, rRing);
569#endif
570
571  return(pResult);
572}
573
574
575// return new poly = pMonom * pPoly; preserve pPoly and pMonom.
576static poly sca_mm_Mult_pp(const poly pMonom, const poly pPoly, const ring rRing)
577{
578  assume( rIsSCA(rRing) );
579
580#ifdef PDEBUG
581//  Print("sca_mm_Mult_pp\n"); // !
582
583  p_Test(pPoly, rRing);
584  p_Test(pMonom, rRing);
585#endif
586
587  if ((pPoly==NULL) || (pMonom==NULL)) return NULL;
588
589  assume( (pPoly != NULL) && (pMonom !=NULL));
590
591  const int iComponentMonomM = p_GetComp(pMonom, rRing);
592
593  poly pResult = NULL;
594  poly* ppPrev = &pResult;
595
596  for( poly p = pPoly; p!= NULL; pIter(p) )
597  {
598#ifdef PDEBUG
599    p_Test(p, rRing);
600#endif
601    const int iComponent = p_GetComp(p, rRing);
602
603    if( iComponentMonomM!=0 )
604    {
605      if( iComponent!=0 ) // TODO: make global if on iComponentMonomM =?= 0
606      {
607        // REPORT_ERROR
608        Werror("sca_mm_Mult_pp: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
609        // what should we do further?!?
610
611        p_Delete( &pResult, rRing); // delete the result
612        return NULL;
613      }
614#ifdef PDEBUG
615      if(iComponent==0 )
616      {
617        dReportError("sca_mm_Mult_pp: Multiplication in the left module from the right!");
618//        PrintS("mm = "); p_Write(pMonom, rRing);
619//        PrintS("pp = "); p_Write(pPoly, rRing);
620//        assume(iComponent!=0);
621      }
622#endif
623    }
624
625    // terms will be in the same order because of quasi-ordering!
626    poly v = sca_mm_Mult_mm(pMonom, p, rRing);
627
628    if( v != NULL )
629    {
630      *ppPrev = v;
631      ppPrev = &pNext(*ppPrev); // nice line ;-)
632    }
633  }
634
635#ifdef PDEBUG
636  p_Test(pResult,rRing);
637#endif
638
639  return(pResult);
640}
641
642
643// return poly = pMonom * pPoly; preserve pMonom, destroy or reuse pPoly.
644static poly sca_mm_Mult_p(const poly pMonom, poly pPoly, const ring rRing) // !!!!! the MOST used procedure !!!!!
645{
646  assume( rIsSCA(rRing) );
647
648#ifdef PDEBUG
649  p_Test(pPoly, rRing);
650  p_Test(pMonom, rRing);
651#endif
652
653  if( pPoly == NULL )
654    return NULL;
655
656  assume(pMonom!=NULL);
657  //if( pMonom == NULL )
658  //{
659  //  // pPoly != NULL =>
660  //  p_Delete( &pPoly, rRing );
661  //  return NULL;
662  //}
663
664  const int iComponentMonomM = p_GetComp(pMonom, rRing);
665
666  poly p = pPoly; poly* ppPrev = &pPoly;
667
668  loop
669  {
670#ifdef PDEBUG
671    if( !p_Test(p, rRing) )
672    {
673      PrintS("p is wrong!");
674      p_Write(p,rRing);
675    }
676#endif
677
678    const int iComponent = p_GetComp(p, rRing);
679
680    if( iComponentMonomM!=0 )
681    {
682      if( iComponent!=0 )
683      {
684        // REPORT_ERROR
685        Werror("sca_mm_Mult_p: exponent mismatch %d and %d\n", iComponent, iComponentMonomM);
686        // what should we do further?!?
687
688        p_Delete( &pPoly, rRing); // delete the result
689        return NULL;
690      }
691#ifdef PDEBUG
692      if(iComponent==0)
693      {
694        dReportError("sca_mm_Mult_p: Multiplication in the left module from the right!");
695//        PrintS("mm = "); p_Write(pMonom, rRing);
696//        PrintS("p = "); p_Write(pPoly, rRing);
697//        assume(iComponent!=0);
698      }
699#endif
700    }
701
702    // terms will be in the same order because of quasi-ordering!
703    poly v = sca_mm_Mult_m(pMonom, p, rRing);
704
705    if( v != NULL )
706    {
707      ppPrev = &pNext(p);
708
709    // *p is changed if v != NULL ( p == v )
710      pIter(p);
711
712      if( p == NULL )
713        break;
714    }
715    else
716    { // Upps! Zero!!! we must kill this term!
717      p = p_LmDeleteAndNext(p, rRing);
718
719      *ppPrev = p;
720
721      if( p == NULL )
722        break;
723    }
724  }
725
726#ifdef PDEBUG
727    if( !p_Test(pPoly, rRing) )
728    {
729      PrintS("pPoly is wrong!");
730      p_Write(pPoly, rRing);
731    }
732#endif
733
734  return(pPoly);
735}
736
737//-----------------------------------------------------------------------------------//
738
739#ifdef PDEBUG
740#endif
741
742
743
744
745//-----------------------------------------------------------------------------------//
746
747// GB computation routines:
748
749/*4
750* creates the S-polynomial of p1 and p2
751* does not destroy p1 and p2
752*/
753poly sca_SPoly( const poly p1, const poly p2, const ring r )
754{
755  assume( rIsSCA(r) );
756
757  const long lCompP1 = p_GetComp(p1,r);
758  const long lCompP2 = p_GetComp(p2,r);
759
760  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
761  {
762#ifdef PDEBUG
763    dReportError("sca_SPoly: different non-zero components!\n");
764#endif
765    return(NULL);
766  }
767
768  poly pL = p_Lcm(p1, p2, si_max(lCompP1, lCompP2), r);       // pL = lcm( lm(p1), lm(p2) )
769
770  poly m1 = p_One( r);
771  p_ExpVectorDiff(m1, pL, p1, r);                  // m1 = pL / lm(p1)
772
773  //p_SetComp(m1,0,r);
774  //p_Setm(m1,r);
775#ifdef PDEBUG
776  p_Test(m1,r);
777#endif
778
779
780  poly m2 = p_One( r);
781  p_ExpVectorDiff (m2, pL, p2, r);                  // m2 = pL / lm(p2)
782
783  //p_SetComp(m2,0,r);
784  //p_Setm(m2,r);
785#ifdef PDEBUG
786  p_Test(m2,r);
787#endif
788
789  p_Delete(&pL,r);
790
791  number C1  = n_Copy(p_GetCoeff(p1,r),r);      // C1 = lc(p1)
792  number C2  = n_Copy(p_GetCoeff(p2,r),r);      // C2 = lc(p2)
793
794  number C = n_Gcd(C1,C2,r);                     // C = gcd(C1, C2)
795
796  if (!n_IsOne(C, r))                              // if C != 1
797  {
798    C1=n_Div(C1, C, r);                              // C1 = C1 / C
799    C2=n_Div(C2, C, r);                              // C2 = C2 / C
800  }
801
802  n_Delete(&C,r); // destroy the number C
803
804  const int iSignSum = sca_Sign_mm_Mult_mm (m1, p1, r) + sca_Sign_mm_Mult_mm (m2, p2, r);
805  // zero if different signs
806
807  assume( (iSignSum*iSignSum == 0) || (iSignSum*iSignSum == 4) );
808
809  if( iSignSum != 0 ) // the same sign!
810    C2=n_Neg (C2, r);
811
812  p_SetCoeff(m1, C2, r);                           // lc(m1) = C2!!!
813  p_SetCoeff(m2, C1, r);                           // lc(m2) = C1!!!
814
815  poly tmp1 = nc_mm_Mult_pp (m1, pNext(p1), r);    // tmp1 = m1 * tail(p1),
816  p_Delete(&m1,r);  //  => n_Delete(&C1,r);
817
818  poly tmp2 = nc_mm_Mult_pp (m2, pNext(p2), r);    // tmp1 = m2 * tail(p2),
819  p_Delete(&m2,r);  //  => n_Delete(&C1,r);
820
821  poly spoly = p_Add_q (tmp1, tmp2, r); // spoly = spoly(lt(p1), lt(p2)) + m1 * tail(p1), delete tmp1,2
822
823  if (spoly!=NULL) p_Cleardenom (spoly, r);
824//  if (spoly!=NULL) p_Content (spoly); // r?
825
826#ifdef PDEBUG
827  p_Test (spoly, r);
828#endif
829
830  return(spoly);
831}
832
833
834
835
836/*2
837* reduction of p2 with p1
838* do not destroy p1, but p2
839* p1 divides p2 -> for use in NF algorithm
840*/
841poly sca_ReduceSpoly(const poly p1, poly p2, const ring r)
842{
843  assume( rIsSCA(r) );
844
845  assume( p1 != NULL );
846
847  const long lCompP1 = p_GetComp (p1, r);
848  const long lCompP2 = p_GetComp (p2, r);
849
850  if ((lCompP1!=lCompP2) && (lCompP1!=0) && (lCompP2!=0))
851  {
852#ifdef PDEBUG
853    dReportError("sca_ReduceSpoly: different non-zero components!");
854#endif
855    return(NULL);
856  }
857
858  poly m = p_ISet (1, r);
859  p_ExpVectorDiff (m, p2, p1, r);                      // m = lm(p2) / lm(p1)
860  //p_Setm(m,r);
861#ifdef PDEBUG
862  p_Test (m,r);
863#endif
864
865  number C1 = n_Copy( p_GetCoeff(p1, r), r);
866  number C2 = n_Copy( p_GetCoeff(p2, r), r);
867
868  /* GCD stuff */
869  number C = n_Gcd(C1, C2, r);
870
871  if (!n_IsOne(C, r))
872  {
873    C1 = n_Div(C1, C, r);
874    C2 = n_Div(C2, C, r);
875  }
876  n_Delete(&C,r);
877
878  const int iSign = sca_Sign_mm_Mult_mm( m, p1, r );
879
880  if(iSign == 1)
881    C2 = n_Neg(C2,r);
882
883  p_SetCoeff(m, C2, r);
884
885#ifdef PDEBUG
886  p_Test(m,r);
887#endif
888
889  p2 = p_LmDeleteAndNext( p2, r );
890
891  p2 = p_Mult_nn(p2, C1, r); // p2 !!!
892  n_Delete(&C1,r);
893
894  poly T = nc_mm_Mult_pp(m, pNext(p1), r);
895  p_Delete(&m, r);
896
897  p2 = p_Add_q(p2, T, r);
898
899  if ( p2!=NULL ) p_Content(p2,r);
900
901#ifdef PDEBUG
902  p_Test(p2,r);
903#endif
904
905  return(p2);
906}
907
908// should be used only inside nc_SetupQuotient!
909// Check whether this our case:
910//  1. rG is  a commutative polynomial ring \otimes anticommutative algebra
911//  2. factor ideal rGR->qideal contains squares of all alternating variables.
912//
913// if yes, make rGR a super-commutative algebra!
914// NOTE: Factors of SuperCommutative Algebras are supported this way!
915//
916//  rG == NULL means that there is no separate base G-algebra in this case take rGR == rG
917
918// special case: bCopy == true (default value: false)
919// meaning: rGR copies structure from rG
920// (maybe with some minor changes, which don't change the type!)
921bool sca_SetupQuotient(ring rGR, ring rG, bool bCopy)
922{
923
924  //////////////////////////////////////////////////////////////////////////
925  // checks...
926  //////////////////////////////////////////////////////////////////////////
927  if( rG == NULL )
928    rG = rGR;
929
930  assume(rGR != NULL);
931  assume(rG  != NULL);
932  assume(rIsPluralRing(rG));
933
934#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
935  PrintS("sca_SetupQuotient(rGR, rG, bCopy)");
936
937  {
938    PrintS("\nrG: \n"); rWrite(rG);
939    PrintS("\nrGR: \n"); rWrite(rGR);
940    PrintLn();
941  }
942#endif
943
944
945  if(bCopy)
946  {
947    if(rIsSCA(rG) && (rG != rGR))
948      return sca_Force(rGR, scaFirstAltVar(rG), scaLastAltVar(rG));
949    else
950      return false;
951  }
952
953  assume(!bCopy);
954
955  const int N = rG->N;
956
957  if(N < 2)
958    return false;
959
960
961//  if( (ncRingType(rG) != nc_skew) || (ncRingType(rG) != nc_comm) )
962//    return false;
963
964#if OUTPUT
965  PrintS("sca_SetupQuotient: qring?\n");
966#endif
967
968  if(rGR->qideal == NULL) // there should be a factor!
969    return false;
970
971#if OUTPUT
972  PrintS("sca_SetupQuotient: qideal!!!\n");
973#endif
974
975//  if((rG->qideal != NULL) && (rG != rGR) ) // we cannot change from factor to factor at the time, sorry!
976//    return false;
977
978
979  int iAltVarEnd = -1;
980  int iAltVarStart   = N+1;
981
982  const nc_struct* NC = rG->GetNC();
983  const ring rBase = rG; //NC->basering;
984  const matrix C   = NC->C; // live in rBase!
985  const matrix D   = NC->D; // live in rBase!
986
987#if OUTPUT
988  PrintS("sca_SetupQuotient: AltVars?!\n");
989#endif
990
991  for(int i = 1; i < N; i++)
992  {
993    for(int j = i + 1; j <= N; j++)
994    {
995      if( MATELEM(D,i,j) != NULL) // !!!???
996      {
997#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
998        Print("Nonzero D[%d, %d]\n", i, j);
999#endif
1000        return false;
1001      }
1002
1003
1004      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1005      number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1006
1007      if( n_IsMOne(c, rBase) ) // !!!???
1008      {
1009        if( i < iAltVarStart)
1010          iAltVarStart = i;
1011
1012        if( j > iAltVarEnd)
1013          iAltVarEnd = j;
1014      } else
1015      {
1016        if( !n_IsOne(c, rBase) )
1017        {
1018#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1019          Print("Wrong Coeff at: [%d, %d]\n", i, j);
1020#endif
1021          return false;
1022        }
1023      }
1024    }
1025  }
1026
1027#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1028  Print("AltVars?1: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1029#endif
1030
1031
1032  if( (iAltVarEnd == -1) || (iAltVarStart == (N+1)) )
1033    return false; // either no alternating varables, or a single one => we are in commutative case!
1034
1035
1036  for(int i = 1; i < N; i++)
1037  {
1038    for(int j = i + 1; j <= N; j++)
1039    {
1040      assume(MATELEM(C,i,j) != NULL); // after CallPlural!
1041      number c = p_GetCoeff(MATELEM(C,i,j), rBase);
1042
1043      if( (iAltVarStart <= i) && (j <= iAltVarEnd) ) // S <= i < j <= E
1044      { // anticommutative part
1045        if( !n_IsMOne(c, rBase) )
1046        {
1047#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1048           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1049#endif
1050          return false;
1051        }
1052      } else
1053      { // should commute
1054        if( !n_IsOne(c, rBase) )
1055        {
1056#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1057           Print("Wrong Coeff at: [%d, %d]\n", i, j);
1058#endif
1059          return false;
1060        }
1061      }
1062    }
1063  }
1064
1065#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1066  Print("AltVars!?: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1067#endif
1068
1069  assume( 1            <= iAltVarStart );
1070  assume( iAltVarStart < iAltVarEnd   );
1071  assume( iAltVarEnd   <= N            );
1072
1073
1074//  ring rSaveRing = assureCurrentRing(rG);
1075
1076
1077  assume(rGR->qideal != NULL);
1078  assume(rGR->N == rG->N);
1079//  assume(rG->qideal == NULL); // ?
1080
1081  const ideal idQuotient = rGR->qideal;
1082
1083
1084#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1085  PrintS("Analyzing quotient ideal:\n");
1086  idPrint(idQuotient); // in rG!!!
1087#endif
1088
1089
1090  // check for
1091  // y_{iAltVarStart}^2, y_{iAltVarStart+1}^2, \ldots, y_{iAltVarEnd}^2  (iAltVarEnd > iAltVarStart)
1092  // to be within quotient ideal.
1093
1094  bool bSCA = true;
1095
1096  int b = N+1;
1097  int e = -1;
1098
1099  if(rIsSCA(rG))
1100  {
1101    b = si_min(b, scaFirstAltVar(rG));
1102    e = si_max(e, scaLastAltVar(rG));
1103
1104#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1105    Print("AltVars!?: [%d, %d]\n", b, e);
1106#endif
1107  }
1108
1109  for ( int i = iAltVarStart; (i <= iAltVarEnd) && bSCA; i++ )
1110    if( (i < b) || (i > e) ) // otherwise it's ok since rG is an SCA!
1111  {
1112    poly square = p_One( rG);
1113    p_SetExp(square, i, 2, rG); // square = var(i)^2.
1114    p_Setm(square, rG);
1115
1116    // square = NF( var(i)^2 | Q )
1117    // NOTE: rSaveRing == currRing now!
1118    // NOTE: there is no better way to check this in general!
1119    //
1120//////    TODO!!!!!
1121//////    extern poly kNF(ideal I, ideal Q, poly f, int a, int b, const ring r);
1122//////    square = kNF(idQuotient, NULL, square, 0, 1, rG); // must ran in currRing == rG!
1123
1124    if( square != NULL ) // var(i)^2 is not in Q?
1125    {
1126      p_Delete(&square, rG);
1127      bSCA = false;
1128      break;
1129    }
1130  }
1131
1132//  assureCurrentRing(rSaveRing);
1133
1134  if(!bSCA)
1135  {
1136    WarnS("sca_SetupQuotient: Sorry can not use NF: unfixed implementational issue :(");
1137    return false;
1138  }
1139   
1140
1141
1142#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1143  Print("ScaVars!: [%d, %d]\n", iAltVarStart, iAltVarEnd);
1144#endif
1145
1146
1147  //////////////////////////////////////////////////////////////////////////
1148  // ok... here we go. let's setup it!!!
1149  //////////////////////////////////////////////////////////////////////////
1150  ideal tempQ = id_KillSquares(idQuotient, iAltVarStart, iAltVarEnd, rG); // in rG!!!
1151
1152
1153#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1154  PrintS("Quotient: \n");
1155  iiWriteMatrix((matrix)idQuotient,"__",1);
1156  PrintS("tempSCAQuotient: \n");
1157  iiWriteMatrix((matrix)tempQ,"__",1);
1158#endif
1159
1160  idSkipZeroes( tempQ );
1161
1162  ncRingType( rGR, nc_exterior );
1163
1164  scaFirstAltVar( rGR, iAltVarStart );
1165  scaLastAltVar( rGR, iAltVarEnd );
1166
1167  if( idIs0(tempQ) )
1168    rGR->GetNC()->SCAQuotient() = NULL;
1169  else
1170    rGR->GetNC()->SCAQuotient() = idrMoveR(tempQ, rG, rGR); // deletes tempQ!
1171
1172  nc_p_ProcsSet(rGR, rGR->p_Procs); // !!!!!!!!!!!!!!!!!
1173
1174
1175#if ((defined(PDEBUG) && OUTPUT) || MYTEST)
1176  PrintS("SCAQuotient: \n");
1177  if(tempQ != NULL)
1178    iiWriteMatrix((matrix)tempQ,"__",1);
1179  else
1180    PrintS("(NULL)\n");
1181#endif
1182
1183  return true;
1184}
1185
1186
1187bool sca_Force(ring rGR, int b, int e)
1188{
1189  assume(rGR != NULL);
1190  assume(rIsPluralRing(rGR));
1191  assume(!rIsSCA(rGR));
1192
1193  const int N = rGR->N;
1194
1195//  ring rSaveRing = currRing;
1196//  if(rSaveRing != rGR)
1197//    rChangeCurrRing(rGR);
1198
1199  const ideal idQuotient = rGR->qideal;
1200
1201  ideal tempQ = idQuotient;
1202
1203  if( b <= N && e >= 1 )
1204    tempQ = id_KillSquares(idQuotient, b, e, rGR);
1205
1206  idSkipZeroes( tempQ );
1207
1208  ncRingType( rGR, nc_exterior );
1209
1210  if( idIs0(tempQ) )
1211    rGR->GetNC()->SCAQuotient() = NULL;
1212  else
1213    rGR->GetNC()->SCAQuotient() = tempQ;
1214
1215
1216  scaFirstAltVar( rGR, b );
1217  scaLastAltVar( rGR, e );
1218
1219
1220  nc_p_ProcsSet(rGR, rGR->p_Procs);
1221
1222//  if(rSaveRing != rGR)
1223//    rChangeCurrRing(rSaveRing);
1224
1225  return true;
1226}
1227
1228// return x_i * pPoly; preserve pPoly.
1229poly sca_pp_Mult_xi_pp(unsigned int i, const poly pPoly, const ring rRing)
1230{
1231  assume(1 <= i);
1232  assume(i <= (unsigned int)rRing->N);
1233
1234  if(rIsSCA(rRing))
1235    return sca_xi_Mult_pp(i, pPoly, rRing);
1236
1237
1238
1239  poly xi =  p_One( rRing);
1240  p_SetExp(xi, i, 1, rRing);
1241  p_Setm(xi, rRing);
1242
1243  poly pResult = pp_Mult_qq(xi, pPoly, rRing);
1244
1245  p_Delete( &xi, rRing);
1246
1247  return pResult;
1248
1249}
1250
1251void sca_p_ProcsSet(ring rGR, p_Procs_s* p_Procs)
1252{
1253
1254  // "commutative" procedures:
1255  rGR->p_Procs->p_Mult_mm     = sca_p_Mult_mm;
1256  rGR->p_Procs->pp_Mult_mm    = sca_pp_Mult_mm;
1257
1258  p_Procs->p_Mult_mm          = sca_p_Mult_mm;
1259  p_Procs->pp_Mult_mm         = sca_pp_Mult_mm;
1260
1261  // non-commutaitve
1262  rGR->GetNC()->p_Procs.mm_Mult_p   = sca_mm_Mult_p;
1263  rGR->GetNC()->p_Procs.mm_Mult_pp  = sca_mm_Mult_pp;
1264
1265//   rGR->GetNC()->p_Procs.SPoly         = sca_SPoly;
1266//   rGR->GetNC()->p_Procs.ReduceSPoly   = sca_ReduceSpoly;
1267
1268#if 0
1269
1270        // Multiplication procedures:
1271
1272        p_Procs->p_Mult_mm   = sca_p_Mult_mm;
1273        _p_procs->p_Mult_mm  = sca_p_Mult_mm;
1274
1275        p_Procs->pp_Mult_mm  = sca_pp_Mult_mm;
1276        _p_procs->pp_Mult_mm = sca_pp_Mult_mm;
1277
1278        r->GetNC()->mmMultP()     = sca_mm_Mult_p;
1279        r->GetNC()->mmMultPP()    = sca_mm_Mult_pp;
1280
1281/*
1282        // ??????????????????????????????????????? coefficients swell...
1283        r->GetNC()->SPoly()         = sca_SPoly;
1284        r->GetNC()->ReduceSPoly()   = sca_ReduceSpoly;
1285*/
1286//         r->GetNC()->BucketPolyRed() = gnc_kBucketPolyRed;
1287//         r->GetNC()->BucketPolyRed_Z()= gnc_kBucketPolyRed_Z;
1288#endif
1289}
1290
1291
1292// bi-Degree (x, y) of monomial "m"
1293// x-es and y-s are weighted by wx and wy resp.
1294// [optional] components have weights by wCx and wCy.
1295static inline void m_GetBiDegree(const poly m,
1296  const intvec *wx, const intvec *wy,
1297  const intvec *wCx, const intvec *wCy,
1298  int& dx, int& dy, const ring r)
1299{
1300  const unsigned int N  = r->N;
1301
1302  p_Test(m, r);
1303
1304  assume( wx != NULL );
1305  assume( wy != NULL );
1306
1307  assume( wx->cols() == 1 );
1308  assume( wy->cols() == 1 );
1309
1310  assume( (unsigned int)wx->rows() >= N );
1311  assume( (unsigned int)wy->rows() >= N );
1312
1313  int x = 0;
1314  int y = 0;
1315
1316  for(int i = N; i > 0; i--)
1317  {
1318    const int d = p_GetExp(m, i, r);
1319    x += d * (*wx)[i-1];
1320    y += d * (*wy)[i-1];
1321  }
1322
1323  if( (wCx != NULL) && (wCy != NULL) )
1324  {
1325    const int c = p_GetComp(m, r);
1326
1327    if( wCx->range(c) )
1328      x += (*wCx)[c];
1329
1330    if( wCy->range(c) )
1331      x += (*wCy)[c];
1332  }
1333
1334  dx = x;
1335  dy = y;
1336}
1337
1338// returns true if polynom p is bi-homogenous with respect to the given weights
1339// simultaneously sets bi-Degree
1340bool p_IsBiHomogeneous(const poly p,
1341  const intvec *wx, const intvec *wy,
1342  const intvec *wCx, const intvec *wCy,
1343  int &dx, int &dy,
1344  const ring r)
1345{
1346  if( p == NULL )
1347  {
1348    dx = 0;
1349    dy = 0;
1350    return true;
1351  }
1352
1353  poly q = p;
1354
1355
1356  int ddx, ddy;
1357
1358  m_GetBiDegree( q, wx, wy, wCx, wCy, ddx, ddy, r); // get bi degree of lm(p)
1359
1360  pIter(q);
1361
1362  for(; q != NULL; pIter(q) )
1363  {
1364    int x, y;
1365
1366    m_GetBiDegree( q, wx, wy, wCx, wCy, x, y, r); // get bi degree of q
1367
1368    if ( (x != ddx) || (y != ddy) ) return false;
1369  }
1370
1371  dx = ddx;
1372  dy = ddy;
1373
1374  return true;
1375}
1376
1377
1378// returns true if id is bi-homogenous without respect to the given weights
1379bool id_IsBiHomogeneous(const ideal id,
1380  const intvec *wx, const intvec *wy,
1381  const intvec *wCx, const intvec *wCy,
1382  const ring r)
1383{
1384  if (id == NULL) return true; // zero ideal
1385
1386  const int iSize = IDELEMS(id);
1387
1388  if (iSize == 0) return true;
1389
1390  bool b = true;
1391  int x, y;
1392
1393  for(int i = iSize - 1; (i >= 0 ) && b; i--)
1394    b = p_IsBiHomogeneous(id->m[i], wx, wy, wCx, wCy, x, y, r);
1395
1396  return b;
1397}
1398
1399
1400// returns an intvector with [nvars(r)] integers [1/0]
1401// 1 - for commutative variables
1402// 0 - for anticommutative variables
1403intvec *ivGetSCAXVarWeights(const ring r)
1404{
1405  const unsigned int N  = r->N;
1406
1407  const int CommutativeVariable = 0; // bug correction!
1408  const int AntiCommutativeVariable = 0;
1409
1410  intvec* w = new intvec(N, 1, CommutativeVariable);
1411
1412  if(AntiCommutativeVariable != CommutativeVariable)
1413  if( rIsSCA(r) )
1414  {
1415    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
1416    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
1417
1418    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
1419    {
1420      (*w)[i-1] = AntiCommutativeVariable;
1421    }
1422  }
1423
1424  return w;
1425}
1426
1427
1428// returns an intvector with [nvars(r)] integers [1/0]
1429// 0 - for commutative variables
1430// 1 - for anticommutative variables
1431intvec *ivGetSCAYVarWeights(const ring r)
1432{
1433  const unsigned int N  = r->N;
1434
1435  const int CommutativeVariable = 0;
1436  const int AntiCommutativeVariable = 1;
1437
1438  intvec* w = new intvec(N, 1, CommutativeVariable);
1439
1440  if(AntiCommutativeVariable != CommutativeVariable)
1441  if( rIsSCA(r) )
1442  {
1443    const unsigned int m_iFirstAltVar = scaFirstAltVar(r);
1444    const unsigned int m_iLastAltVar  = scaLastAltVar(r);
1445
1446    for (unsigned int i = m_iFirstAltVar; i<= m_iLastAltVar; i++)
1447    {
1448      (*w)[i-1] = AntiCommutativeVariable;
1449    }
1450  }
1451  return w;
1452}
1453
1454
1455
1456
1457// reduce term lt(m) modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar:
1458// either create a copy of m or return NULL
1459static inline poly m_KillSquares(const poly m,
1460  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
1461  const ring r)
1462{
1463#ifdef PDEBUG
1464  p_Test(m, r);
1465  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
1466
1467#if 0
1468  Print("m_KillSquares, m = "); // !
1469  p_Write(m, r);
1470#endif
1471#endif
1472
1473  assume( m != NULL );
1474
1475  for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
1476    if( p_GetExp(m, k, r) > 1 )
1477      return NULL;
1478
1479  return p_Head(m, r);
1480}
1481
1482
1483// reduce polynomial p modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
1484// returns a new poly!
1485poly p_KillSquares(const poly p,
1486  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
1487  const ring r)
1488{
1489#ifdef PDEBUG
1490  p_Test(p, r);
1491
1492  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
1493
1494#if 0
1495  Print("p_KillSquares, p = "); // !
1496  p_Write(p, r);
1497#endif
1498#endif
1499
1500
1501  if( p == NULL )
1502    return NULL;
1503
1504  poly pResult = NULL;
1505  poly* ppPrev = &pResult;
1506
1507  for( poly q = p; q!= NULL; pIter(q) )
1508  {
1509#ifdef PDEBUG
1510    p_Test(q, r);
1511#endif
1512
1513    // terms will be in the same order because of quasi-ordering!
1514    poly v = m_KillSquares(q, iFirstAltVar, iLastAltVar, r);
1515
1516    if( v != NULL )
1517    {
1518      *ppPrev = v;
1519      ppPrev = &pNext(v);
1520    }
1521
1522  }
1523
1524#ifdef PDEBUG
1525  p_Test(pResult, r);
1526#if 0
1527  Print("p_KillSquares => "); // !
1528  p_Write(pResult, r);
1529#endif
1530#endif
1531
1532  return(pResult);
1533}
1534
1535
1536
1537
1538// reduces ideal id modulo <y_i^2> , i = iFirstAltVar .. iLastAltVar
1539// returns the reduced ideal or zero ideal.
1540ideal id_KillSquares(const ideal id,
1541  const unsigned int iFirstAltVar, const unsigned int iLastAltVar,
1542  const ring r, const bool bSkipZeroes)
1543{
1544  if (id == NULL) return id; // zero ideal
1545
1546  assume( (iFirstAltVar >= 1) && (iLastAltVar <= r->N) && (iFirstAltVar <= iLastAltVar) );
1547
1548  const int iSize = IDELEMS(id);
1549
1550  if (iSize == 0) return id;
1551
1552  ideal temp = idInit(iSize, id->rank);
1553
1554#if 0
1555   PrintS("<id_KillSquares>\n");
1556  {
1557    PrintS("ideal id: \n");
1558    for (int i = 0; i < IDELEMS(id); i++)
1559    {
1560      Print("; id[%d] = ", i+1);
1561      p_Write(id->m[i], r);
1562    }
1563    PrintS(";\n");
1564    PrintLn();
1565  }
1566#endif
1567
1568
1569  for (int j = 0; j < iSize; j++)
1570    temp->m[j] = p_KillSquares(id->m[j], iFirstAltVar, iLastAltVar, r);
1571
1572  if( bSkipZeroes )
1573    idSkipZeroes(temp);
1574
1575#if 0
1576   PrintS("<id_KillSquares>\n");
1577  {
1578    PrintS("ideal temp: \n");
1579    for (int i = 0; i < IDELEMS(temp); i++)
1580    {
1581      Print("; temp[%d] = ", i+1);
1582      p_Write(temp->m[i], r);
1583    }
1584    PrintS(";\n");
1585    PrintLn();
1586  }
1587   PrintS("</id_KillSquares>\n");
1588#endif
1589
1590//  temp->rank = idRankFreeModule(temp, r);
1591
1592  return temp;
1593}
1594
1595
1596
1597
1598#endif
Note: See TracBrowser for help on using the repository browser.