source: git/libfac/charset/charset.cc @ 405ebc

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