source: git/libfac/charset/charset.cc @ 18500b

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