source: git/libfac/charset/charset.cc @ f152c5

jengelh-datetimespielwiese
Last change on this file since f152c5 was f152c5, checked in by Hans Schoenemann <hannes@…>, 12 years ago
removed Factorize2 git-svn-id: file:///usr/local/Singular/svn/trunk@14303 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.8 KB
Line 
1////////////////////////////////////////////////////////////
2// emacs edit mode for this file is -*- C++ -*-
3/* $Id$ */
4/////////////////////////////////////////////////////////////
5// FACTORY - Includes
6#include <factory.h>
7#ifdef HAVE_IOSTREAM
8#define CIN std::cin
9#else
10#define CIN cin
11#endif
12// Factor - Includes
13#include <SqrFree.h>
14#include <Factor.h>
15#include <interrupt.h>
16// Charset - Includes
17#include "csutil.h"
18#include "algfactor.h"
19#include "alg_factor.h"
20// Some CC's need this:
21#include "charset.h"
22
23#ifdef BASICSETDEBUG
24#  define DEBUGOUTPUT
25#else
26#  undef DEBUGOUTPUT
27#endif
28
29#include "debug.h"
30#include "timing.h"
31TIMING_DEFINE_PRINT(subfactorize_time);
32
33
34// forward declarations:
35static CFList     irras(CFList & AS, int &ja, CanonicalForm & reducible);
36
37#ifdef BASICSETDEBUG
38#  define DEBUGOUTPUT
39#else
40#  undef DEBUGOUTPUT
41#endif
42#include "debug.h"
43
44// the next computes a characteristic set (a basic set in Wang's sense)
45CFList
46BasicSet( const CFList &PS )
47{
48    CFList QS = PS, BS, RS;
49    CanonicalForm b;
50    int cb;
51
52    DEBOUTLN(CERR, "BasicSet: called with ps= ", PS);
53    if ( PS.length() < 2 ) return PS;
54    while ( ! QS.isEmpty() ) {
55        b = lowestRank( QS );
56        cb = rank( b );
57        DEBOUTLN(CERR, "BasicSet: choose b  = ", b);
58        DEBOUTLN(CERR, "BasicSet: it's rank = ", cb);
59        BS=Union(CFList(b),BS);//BS.append( b );
60        if ( rank( b ) == 0 )
61            return Union(PS, CFList(b)) ; // b should be the first elem!
62        else {
63            RS = CFList();
64            // QS:= {q \in QS -{B} | q is reduced wrt b}
65            // We can process whole QS, because b is not reduced wrt. b
66            for ( CFListIterator i = QS; i.hasItem(); ++i )
67                if ( degree( i.getItem(), cb ) < degree( b ) )
68                    //RS.append( i.getItem() );
69                    RS = Union(CFList(i.getItem()),RS);
70            QS = RS;
71        }
72    }
73    DEBOUTLN(CERR, "BasicSet: returning bs= ", BS);
74    return BS;
75}
76
77int
78checkok( const CFList & PS, CFList & FS2){
79  CanonicalForm elem;
80
81  for ( CFListIterator i=PS; i.hasItem(); i++){
82    elem= i.getItem();
83    for (CFListIterator j=FS2; j.hasItem(); j++){
84      if (elem == j.getItem()){
85         //        FS2= Difference(FS2,CFList(elem));
86         return 0;
87      }
88    }
89  }
90  return 1;
91}
92
93#ifdef MCHARSETNDEBUG
94#  define DEBUGOUTPUT
95#else
96#  undef DEBUGOUTPUT
97#endif
98#include "debug.h"
99
100// The modified CharSet (an extended characteristic set with certain factors
101// canceled; this is a characteristic set in Wang's sense)
102CFList
103MCharSetN( const CFList &PS, PremForm & Remembern ){
104  CFList QS = PS, RS = PS, CSet, OLDCS;
105
106  DEBOUTLN(CERR, "MCharSetN: called with ps= ", PS);
107  while ( ! RS.isEmpty() ) {
108    CSet = BasicSet( QS );
109    OLDCS=CSet;
110    DEBOUTLN(CERR, "MCharSetN: CS= ", CSet);
111//     if ( getNumVars(CSet.getFirst()) > 1 ){
112//       //CSet = removecontent(CSet, Remembern);
113// #ifdef MCHARSETNDEBUG
114//       CERR << "MCharSetN: CSet= " << CSet << "\n";
115// #endif
116//     }
117    Remembern.FS1 = Union(Remembern.FS1, initalset1(CSet));
118    DEBOUTLN(CERR, "MCharSetN: Remembern.FS1= ", Remembern.FS1);
119    DEBOUTLN(CERR, "MCharSetN: Remembern.FS2= ", Remembern.FS2);
120    RS = CFList();
121    if ( rank( CSet.getFirst() ) != 0 ) {
122      CFList D = Difference( QS, CSet );
123      DEBOUT(CERR, "MCharSetN: Difference( ", QS);
124      DEBOUT(CERR, " , ", CSet);
125      DEBOUTLN(CERR, " ) = ", D);
126//CERR << "MCharSetN: Difference( " << QS << " , " << CSet << " ) = " << D << "\n";
127      //PremForm Oldremember=Remembern;
128      //PremForm Newremember=Remembern;
129      for ( CFListIterator i = D; i.hasItem(); ++i ) {
130        CanonicalForm r = Prem( i.getItem(), CSet );
131        DEBOUT(CERR,"MCharSetN: Prem(", i.getItem()  );
132        DEBOUT(CERR, ",", CSet);
133        DEBOUTLN(CERR,") = ", r);
134//CERR << "MCharSetN: Prem("<< i.getItem() << "," << CSet << ") = " << r << "\n";
135        if ( r != 0 ){
136          //removefactor( r, Newremember );
137          removefactor( r, Remembern );
138          //Remembern.FS2 = Union(Remembern.FS2, Newremember.FS2);
139          //Newremember = Oldremember;
140          //if ( cls(r) > 0 )
141            //RS=Union(CFList(r),RS);//RS.append( r );
142            RS=Union(RS,CFList(r));
143        }
144      }
145      if ( ! checkok(RS,Remembern.FS2)) return CFList(CanonicalForm(1));
146      DEBOUTLN(CERR, "MCharSetN: RS= ", RS);
147      //QS = Union( QS, RS );
148      QS = Union(OLDCS,RS);
149      DEBOUTLN(CERR, "MCharSetN: QS= Union(QS,RS)= ", QS);
150    }
151    else{ return CFList(CanonicalForm(1)); }
152  }
153  DEBOUTLN(CERR, "MCharSetN: Removed factors: ", Remembern.FS2);
154  DEBOUTLN(CERR, "MCharSetN: Remembern.FS1: ", Remembern.FS1);
155
156  return CSet;
157}
158
159
160CFList
161mcharset( const CFList &PS, PremForm & Remembern ){
162  CFList cs= MCharSetN(PS, Remembern );
163  CFList rs= remsetb(Difference(PS,cs),cs);
164
165  DEBOUTLN(CERR, "mcharset: cs= ", cs);
166  DEBOUTLN(CERR, "mcharset: rs= ", rs);
167  if ( rs.length() > 0 )
168    cs= mcharset(Union(PS,Union(cs,rs)), Remembern);
169  return cs;
170}
171
172// the "original" extended characteristic set
173CFList
174CharSet( const CFList &PS ){
175  CFList QS = PS, RS = PS, CSet;
176
177  while ( ! RS.isEmpty() ) {
178    CSet = BasicSet( QS );
179    DEBOUTLN(CERR, "CharSet: CSet= ", CSet);
180    RS = CFList();
181    if ( rank( CSet.getFirst() ) != 0 ) {
182      CFList D = Difference( QS, CSet );
183      for ( CFListIterator i = D; i.hasItem(); ++i ) {
184        CanonicalForm r = Prem( i.getItem(), CSet );
185        if ( r != 0 )  RS=Union(CFList(r),RS);//RS.append( r );
186      }
187      QS = Union( QS, RS );
188    }
189  }
190  return CSet;
191}
192
193static CFList
194charseta( const CFList & PS ){
195  CFList QS = PS, RS = PS, CSet;
196
197  while ( ! RS.isEmpty() ) {
198    CSet = CharSet( QS );
199    RS = CFList();
200    if ( rank( CSet.getFirst() ) != 0 ) {
201      CFList D = Difference( QS, CSet );
202      for ( CFListIterator i = D; i.hasItem(); ++i ) {
203        CanonicalForm r = Prem( i.getItem(), CSet );
204        if ( r != 0 )  RS=Union(CFList(r),RS);//RS.append( r );
205      }
206      QS = Union(CSet,Union( QS, RS ));
207    }
208    else return CFList(CanonicalForm(1));
209  }
210  return CSet;
211}
212
213static bool
214contractsub( const CFList & cs1, const CFList & cs2 ){
215  CFListIterator i;
216
217  for ( i=cs1; i.hasItem(); i++ )
218    if ( Prem(i.getItem(),cs2 ) != 0 ) return false;
219  CFList is=initalset1(cs1);
220  for ( i=is; i.hasItem(); i++ )
221    if ( Prem(i.getItem(),cs2 ) == 0 ) return false;
222  return true;
223}
224
225static ListCFList
226contract( const ListCFList & cs){
227  ListCFList mem,ts;
228  CFList iitem,jitem;
229
230  if ( cs.length() < 2 ) return cs;
231
232  for ( ListCFListIterator i=cs; i.hasItem(); i++ ){
233    iitem=i.getItem();
234    if ( ! member(iitem, mem))
235      for ( ListCFListIterator j=i; j.hasItem(); j++){
236        jitem=j.getItem();
237        if ( ! same( iitem, jitem ) )
238          if ( ! member(jitem, mem))
239            if ( contractsub(iitem, jitem) ){
240              ts.append(jitem); mem.append(jitem);
241            }
242            else
243              if ( contractsub(jitem, iitem) ){
244                ts.append(iitem);
245              }
246      }
247  }
248  return Minus(cs,ts);
249}
250
251static ListCFList
252adjoin(const CFList & is, const CFList & qs, const ListCFList & qh ){
253  ListCFList iss,qhi;
254  ListCFListIterator j;
255  CFList iscopy,itt;
256  CFListIterator i;
257  CanonicalForm elem;
258  int ind, length;
259
260  for ( i=is ; i.hasItem(); i++ ){
261    elem=i.getItem();
262    if ( cls(elem) > 0 ) iscopy=Union(CFList(elem),iscopy);
263  }
264  if ( iscopy.isEmpty() ) return iss;
265  qhi = MyDifference(qh,qs);
266  length = qhi.length();
267  for ( i =iscopy; i.hasItem(); i++){
268    itt = Union(qs,CFList(i.getItem()));
269    ind = 0;
270    if ( length > 0 )
271      for ( j=qhi; j.hasItem(); j++ )
272        if ( subset(j.getItem(),itt )) ind=1;
273    if ( ind == 0 ) iss.append(itt);
274  }
275  return iss;
276}
277
278static ListCFList
279adjoinb(const CFList & is, const CFList & qs, const ListCFList & qh ,const CFList & cs){
280  ListCFList iss,qhi;
281  ListCFListIterator j;
282  CFList iscopy,itt;
283  CFListIterator i;
284  CanonicalForm elem;
285  int ind, length;
286
287  for ( i=is ; i.hasItem(); i++ ){
288    elem=i.getItem();
289    if ( cls(elem) > 0 ) iscopy=Union(CFList(elem),iscopy);
290  }
291  if ( iscopy.isEmpty() ) return iss;
292  qhi = MyDifference(qh,qs);
293  length = qhi.length();
294  for ( i =iscopy; i.hasItem(); i++){
295    itt = Union(Union(qs,CFList(i.getItem())),cs);
296    ind = 0;
297    if ( length > 0 )
298      for ( j=qhi; j.hasItem(); j++ )
299        if ( subset(j.getItem(),itt )) ind=1;
300    if ( ind == 0 ) {iss.append(itt);}
301  }
302  return iss;
303}
304
305static ListCFList
306sort( const ListCFList & list_to_sort ){
307  ListCFList output,copy=list_to_sort;
308  CFList l,qs1,elem;
309
310  l = copy.getLast(); copy.removeLast();
311  if ( copy.length() == 0 ){ return ListCFList(l); }
312  for ( ListCFListIterator i=copy ; i.hasItem(); i++ ){
313    elem = i.getItem();
314    if ( elem.length() > l.length() ) {
315      output = MyUnion( ListCFList(l), output);
316      l= elem;
317    }
318    else{ output = MyUnion(output, ListCFList(elem) ); }
319  }
320  //output = MyUnion(ListCFList(l),sort(output));
321  output = MyUnion(ListCFList(l),output);
322  return output;
323}
324
325#ifdef IRRCHARSERIESDEBUG
326#  define DEBUGOUTPUT
327#else
328#  undef DEBUGOUTPUT
329#endif
330#include "debug.h"
331
332ListCFList
333IrrCharSeries( const CFList &PS, int opt ){
334  CanonicalForm reducible,reducible2;
335  CFList qs,cs,factorset,is,ts;
336  ListCFList pi,ppi,qqi,qsi,iss,qhi= ListCFList(PS);
337  int nr_of_iteration=0,ts2,highestlevel=0;
338
339//  CERR << getCharacteristic() << "\n";
340  for ( CFListIterator Ps=PS; Ps.hasItem(); Ps++ )
341    if ( level(Ps.getItem() ) > highestlevel ) highestlevel = level(Ps.getItem()) ;
342//  for ( int xx=1; xx <= highestlevel; xx++)
343//   CERR << Variable(xx) ;
344//  CERR << "\n";
345//  for ( CFListIterator Ps=PS; Ps.hasItem(); Ps++ )
346//    CERR << Ps.getItem() << ", " ;//<< "\n";
347//  CERR <<  "\n";
348  while ( ! qhi.isEmpty() ) {
349    qhi=sort(qhi);
350    DEBOUTLN(CERR, "qhi is: ", qhi);
351    qs=qhi.getFirst();
352    DEBOUTLN(CERR, "qs  is: ", qs);
353    DEBOUTLN(CERR, "ppi is: ", ppi);
354    ListCFList ppi1,ppi2;
355    select(ppi,qs.length(),ppi1,ppi2);
356    DEBOUTLN(CERR, "ppi1 is: ", ppi1);
357    DEBOUTLN(CERR, "ppi2 is: ", ppi2);
358    qqi = MyUnion(ppi2,qqi);
359    DEBOUTLN(CERR, "qqi is: ", qqi);
360    if ( nr_of_iteration == 0 ){ nr_of_iteration += 1; ppi = ListCFList(); }
361    else{ nr_of_iteration += 1; ppi = MyUnion(ListCFList(qs),ppi1); }
362    DEBOUTLN(CERR,"ppi is: ", ppi);
363    PremForm Remembern;
364    cs = MCharSetN(qs,Remembern);
365    DEBOUTLN(CERR, "cs is: ", cs);
366    DEBOUTLN(CERR, "factorset is: ", Remembern.FS2);
367    cs = removecontent(cs,Remembern);
368    factorset=Remembern.FS2;
369    DEBOUTLN(CERR, "cs (after removecontent) is: ", cs);
370    DEBOUTLN(CERR, "factorset is: ", factorset);
371    // Hier: removecontent einfuegen!!!!
372    if ( cls(cs.getFirst()) > 0 ){
373      ts = irras(cs,ts2, reducible);
374
375      // INTERRUPTHANDLER
376      if ( interrupt_handle() ) return ListCFList() ;
377      // INTERRUPTHANDLER
378
379      DEBOUTLN(CERR, "ts is: ", ts);
380      DEBOUTLN(CERR, "ts2 is: ", ts2);
381      // next is preliminary: should be ==0
382      if ( ts2 <= 0 ){ //irreducible
383        if ( ! subset(cs,qs) ){
384          DEBOUTMSG(CERR, "cs is not a subset of qs");
385          cs = charseta(Union(qs,cs));
386          DEBOUTLN(CERR, "new cs is: ", cs);
387        }
388        if ( ! member(cs,pi) ){
389          pi = MyUnion(pi, ListCFList(cs));
390          DEBOUTMSG(CERR, "cs is not a member of pi");
391          DEBOUTLN(CERR, "pi is: ", pi);
392          if ( cls(cs.getFirst()) > 0 ){
393            ts = irras(cs,ts2,reducible);
394
395            // INTERRUPTHANDLER
396            if ( interrupt_handle() ) return ListCFList() ;
397            // INTERRUPTHANDLER
398
399            DEBOUTLN(CERR, "ts is: ", ts);
400            DEBOUTLN(CERR, "ts2 is: ", ts2);
401            // next is preliminary: should be ==0
402            if ( ts2 <= 0 ){ //irreducible
403              qsi = MyUnion(qsi,ListCFList(cs));
404              DEBOUTLN(CERR, "qsi is: ", qsi);
405              if ( cs.length() == highestlevel ){
406                DEBOUTLN(CERR, "cs.length() == nops(ord) :", cs.length());
407                is = factorps(factorset);
408              }
409              else{
410                DEBOUT(CERR,"cs.length() != nops(ord) :", cs.length());
411                DEBOUTLN(CERR, "  nops(ord)= ", highestlevel);
412                is = Union(initalset1(cs),factorps(factorset));
413              }
414              DEBOUTLN(CERR, "is is: ", is);
415              iss = adjoin(is,qs,qqi);
416              DEBOUTLN(CERR, "iss is: ", iss);
417            }
418          }
419          else{ iss = adjoin(factorps(factorset),qs,qqi); }
420        }
421        else{
422          DEBOUTMSG(CERR, "cs is a member of pi");
423          iss = adjoin(factorps(factorset),qs,qqi); }
424        DEBOUTLN(CERR, "iss is: ", iss);
425        DEBOUTLN(CERR, "   factorps(factorset)= ", factorps(factorset));
426        DEBOUTLN(CERR, "   qs= ", qs);
427        DEBOUTLN(CERR, "   qqi= ", qqi);
428      }
429      // next is preliminary: should be !=0
430      if ( ts2 > 0 ){
431        is = factorps(factorset);
432        DEBOUTLN(CERR, "is is: ", is);
433        if ( ts2 > 1 ){
434          // setup cst: need it later for adjoinb
435          CFList cst;
436          for ( CFListIterator i=cs ; i.hasItem(); i++){
437            if ( i.getItem() == reducible ) { break; }
438            else { cst.append(i.getItem()); }
439          }
440          is = Union(initalset1(cst), is);
441          iss = MyUnion(adjoin(is,qs,qqi), adjoinb(ts,qs,qqi,cst));
442        }
443        else{ iss = adjoin(Union(is,ts),qs,qqi); }
444        DEBOUTLN(CERR, "iss is: ", iss);
445      }
446    }
447    else{
448      iss = adjoin(factorps(factorset),qs,qqi);
449      DEBOUTMSG(CERR, "case: cs is a constant.");
450      DEBOUTLN(CERR, "  qs = ", qs);
451      DEBOUTLN(CERR, "  qqi = ", qqi);
452      DEBOUTLN(CERR, "  iss = adjoin(factorps(factorset),qs,qqi) = ",iss);
453    }
454    if ( qhi.length() > 1 ){ qhi.removeFirst(); qhi = MyUnion(qhi,iss); }
455    else{ qhi = iss; }
456    DEBOUTLN(CERR, "iss is: ", iss);
457  }
458  if ( ! qsi.isEmpty() ){
459    DEBOUTLN(CERR, "qsi before contract= ", qsi);
460    if ( opt == 0 ){
461       return contract( qsi );
462    }
463    else { return qsi; }
464  }
465  else{ return ListCFList() ; }
466}
467
468// tests for characteristic sets
469//////////////////////////////////
470#ifdef IRRASDEBUG
471#  define DEBUGOUTPUT
472#else
473#  undef DEBUGOUTPUT
474#endif
475
476#include "debug.h"
477
478static CFList
479irras( CFList & AS, int & ja, CanonicalForm & reducible)
480{
481  CFFList qs;
482  CFList ts,as;
483  CanonicalForm elem;
484  int ind=1,nr=0, success=-1;
485  CFListIterator i;
486
487  ja = 0;
488  DEBOUTLN(CERR, "irras: called with: AS= ", AS);
489  for ( i=AS; i.hasItem(); i++ )
490  {
491    elem = i.getItem();
492    nr += 1;
493    DEBOUT(CERR, "irras: factoring: ", elem);
494    if ( degree(elem) > 1 ) // linear poly's are irreduzible
495    {
496      qs = factorize(elem);
497      // remove a constant
498      if (qs.getFirst().factor().degree()==0) qs.removeFirst();
499    }
500    else
501    {
502      qs=(CFFactor(elem,1));
503    }
504    DEBOUTLN(CERR, "  = ", qs);
505    // INTERRUPTHANDLER
506    if ( interrupt_handle() ) return CFList() ;
507    // INTERRUPTHANDLER
508
509    if ( (qs.length() >= 2 ) || (qs.getFirst().exp() > 1))
510    {
511      DEBOUTLN(CERR, "irras: Setting ind=0, ja= ", nr);
512      ja=nr; ind=0; reducible= elem;
513      break;
514    }
515    //    else{ as.append(elem) ; }
516  }
517  //  CERR << "ind= " << ind << "\n";
518  if ( (ind == 1) ) //&& ( as.length() > 1) )
519  {
520    if ( irreducible(AS) )
521    { // as quasilinear? => irreducible!
522      ja = 0;
523      DEBOUTLN(CERR, "as is irreducible. as= ", AS);
524    }
525    else
526    {
527      i=AS;
528      for ( nr=1; nr< AS.length(); nr++)
529      {
530        as.append(i.getItem());
531        i++;
532        if ( degree(i.getItem()) > 1 )
533        {  // search for a non linear elem
534          elem=i.getItem();
535          //if (as.length()==1)
536          //  qs = Factorize2(elem,as.getFirst());
537          //else
538            qs= newfactoras(elem,as,success);
539          if ( qs.length() > 1 || qs.getFirst().exp() > 1 )
540          { //found elem is reducible
541            reducible=elem;
542            ja=nr+1;
543            break;
544          }
545        }
546      }
547    }
548  }
549  for ( CFFListIterator k=qs; k.hasItem();k++)
550    ts.append(k.getItem().factor());
551  return ts;
552}
553
554///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.