Changeset 5812c69 in git for Singular/fglmhom.cc


Ignore:
Timestamp:
Sep 24, 1998, 11:59:51 AM (26 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
56c52a7879fbc6de62cefba1da86b2edf2aadd4c
Parents:
073d2edeb03013a5c6ed0913687b024305a1427d
Message:
cosmetic changes


git-svn-id: file:///usr/local/Singular/svn/trunk@2517 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/fglmhom.cc

    r073d2e r5812c69  
    11// emacs edit mode for this file is -*- C++ -*-
    2 // $Id: fglmhom.cc,v 1.10 1998-06-15 14:30:06 Singular Exp $
     2// $Id: fglmhom.cc,v 1.11 1998-09-24 09:59:39 Singular Exp $
    33
    44/****************************************
    55*  Computer Algebra System SINGULAR     *
    66****************************************/
    7 /* 
     7/*
    88* ABSTRACT - The FGLM-Algorithm extended for homogeneous ideals.
    99*   Calculates via the hilbert-function a groebner basis.
    … …  
    1818#include "tok.h"
    1919#include "structs.h"
    20 #include "subexpr.h" 
     20#include "subexpr.h"
    2121#include "polys.h"
    2222#include "ideals.h"
    … …  
    6262    homogElem() : v(), dv(), basis(0), destbasis(0), inDest(FALSE) {}
    6363    homogElem( poly m, int b, BOOLEAN ind ) :
    64         basis(b), inDest(ind)
     64        basis(b), inDest(ind)
    6565    {
    66         mon.dm= m;
    67         mon.sm= NULL;
     66        mon.dm= m;
     67        mon.sm= NULL;
    6868    }
    6969};
    7070
    71 struct homogData 
     71struct homogData
    7272{
    7373    ideal sourceIdeal;
    … …  
    8383    int overall;  // nur zum testen.
    8484    int numberofdestbasismonoms;
    85 //     homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL), 
    86 //      numMonoms(0), basisSize(0) {}
     85//     homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL),
     86//         numMonoms(0), basisSize(0) {}
    8787};
    8888
    8989int
    90 hfglmNextdegree( intvec * source, ideal current, int & deg ) 
     90hfglmNextdegree( intvec * source, ideal current, int & deg )
    9191{
    9292    int numelems;
    9393    intvec * newhilb = hHstdSeries( current, NULL, currQuotient );
    9494
    95     loop 
     95    loop
    9696    {
    97         if ( deg < newhilb->length() )
    98         {
    99             if ( deg < source->length() )
    100                 numelems= (*newhilb)[deg]-(*source)[deg];
    101             else
    102                 numelems= (*newhilb)[deg];
    103         }
    104         else
    105         {
    106             if (deg < source->length())
    107                 numelems= -(*source)[deg];
    108             else
    109             {
    110                 deg= 0;
    111                 return 0;
    112             }
    113         }
    114         if (numelems != 0)
    115             return numelems;
    116         deg++;
    117     } 
     97        if ( deg < newhilb->length() )
     98        {
     99            if ( deg < source->length() )
     100                numelems= (*newhilb)[deg]-(*source)[deg];
     101            else
     102                numelems= (*newhilb)[deg];
     103        }
     104        else
     105        {
     106            if (deg < source->length())
     107                numelems= -(*source)[deg];
     108            else
     109            {
     110                deg= 0;
     111                return 0;
     112            }
     113        }
     114        if (numelems != 0)
     115            return numelems;
     116        deg++;
     117    }
    118118    delete newhilb;
    119119}
    120120
    121 void 
     121void
    122122generateMonoms( poly m, int var, int deg, homogData * dat )
    123123{
    124124    if ( var == pVariables ) {
    125         BOOLEAN inSource = FALSE;
    126         BOOLEAN inDest = FALSE;
    127         poly mon = pCopy( m );
    128         pSetExp( mon, var, deg );
    129         pSetm( mon );
    130         ++dat->overall;
    131         int i;
    132         for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
    133             if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
    134                 inSource= TRUE;
    135             }
    136         }
    137         for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
    138             if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
    139                 inDest= TRUE;
    140             }
    141         }
    142         if ( (!inSource) || (!inDest) ) {
    143             int basis = 0;
    144             if ( !inSource )
    145                 basis= ++(dat->basisSize);
    146             if ( !inDest )
    147                 ++dat->numberofdestbasismonoms;
    148             if ( dat->numMonoms == dat->monlistmax ) {
    149                 dat->monlist= (homogElem * )ReAlloc( dat->monlist, (dat->monlistmax)*sizeof( homogElem ), (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
    150                 int k;
    151                 for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
    152                     dat->monlist[k].homogElem();
    153                 dat->monlistmax+= dat->monlistblock;
    154             }
    155             dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
    156             dat->numMonoms++;
    157             if ( inSource && ! inDest ) PROT( "\\" );
    158             if ( ! inSource && inDest ) PROT( "/" );
    159             if ( ! inSource && ! inDest ) PROT( "." );
    160         }
    161         else {
    162             pDelete( & mon );
    163         }
    164         return;
     125        BOOLEAN inSource = FALSE;
     126        BOOLEAN inDest = FALSE;
     127        poly mon = pCopy( m );
     128        pSetExp( mon, var, deg );
     129        pSetm( mon );
     130        ++dat->overall;
     131        int i;
     132        for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
     133            if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
     134                inSource= TRUE;
     135            }
     136        }
     137        for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
     138            if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
     139                inDest= TRUE;
     140            }
     141        }
     142        if ( (!inSource) || (!inDest) ) {
     143            int basis = 0;
     144            if ( !inSource )
     145                basis= ++(dat->basisSize);
     146            if ( !inDest )
     147                ++dat->numberofdestbasismonoms;
     148            if ( dat->numMonoms == dat->monlistmax ) {
     149                dat->monlist= (homogElem * )ReAlloc( dat->monlist, (dat->monlistmax)*sizeof( homogElem ), (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
     150                int k;
     151                for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
     152                    dat->monlist[k].homogElem();
     153                dat->monlistmax+= dat->monlistblock;
     154            }
     155            dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
     156            dat->numMonoms++;
     157            if ( inSource && ! inDest ) PROT( "\\" );
     158            if ( ! inSource && inDest ) PROT( "/" );
     159            if ( ! inSource && ! inDest ) PROT( "." );
     160        }
     161        else {
     162            pDelete( & mon );
     163        }
     164        return;
    165165    }
    166166    else {
    167         poly newm = pCopy( m );
    168         while ( deg >= 0 ) {
    169             generateMonoms( newm, var+1, deg, dat );
    170             pIncrExp( newm, var );
    171             pSetm( newm );
    172             deg--;
    173         }
    174         pDelete( & newm );
     167        poly newm = pCopy( m );
     168        while ( deg >= 0 ) {
     169            generateMonoms( newm, var+1, deg, dat );
     170            pIncrExp( newm, var );
     171            pSetm( newm );
     172            deg--;
     173        }
     174        pDelete( & newm );
    175175    }
    176176    return;
    … …  
    178178
    179179void
    180 mapMonoms( ring oldRing, homogData & dat ) 
     180mapMonoms( ring oldRing, homogData & dat )
    181181{
    182182    int * vperm = (int *)Alloc( (currRing->N + 1)*sizeof(int) );
    … …  
    185185    int s;
    186186    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
    187 //      dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
     187//        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
    188188      // obachman: changed the folowing to reflect the new calling interface of
    189189      // pPermPoly -- Tim please check whether this is correct!
    190         dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );     
     190        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );
    191191    }
    192192}
    193193
    194194void
    195 getVectorRep( homogData & dat ) 
     195getVectorRep( homogData & dat )
    196196{
    197197    // Calculate the NormalForms
    198198    int s;
    199199    for ( s= 0;  s < dat.numMonoms; s++ ) {
    200         if ( dat.monlist[s].inDest == FALSE ) {
    201             fglmVector v;
    202             if ( dat.monlist[s].basis == 0 ) {
    203                 v= fglmVector( dat.basisSize );
    204                 // now the monom is in L(source)
    205                 PROT( "(" );
    206                 poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
    207                 PROT( ")" );
    208                 poly temp = nf;
    209                 while (temp != NULL ) {
    210                     int t;
    211                     for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
    212                         if ( dat.monlist[t].basis > 0 ) {
    213                             if ( pEqual( dat.monlist[t].mon.sm, temp ) ) {
    214                                 number coeff= nCopy( pGetCoeff( temp ) );
    215                                 v.setelem( dat.monlist[t].basis, coeff );
    216                             }
    217                         }
    218                     }
    219                     pIter(temp);
    220                 }
    221                 pDelete( & nf );
    222             }
    223             else {
    224                 PROT( "." );
    225                 v= fglmVector( dat.basisSize, dat.monlist[s].basis );
    226             }
    227             dat.monlist[s].v= v;
    228         }
     200        if ( dat.monlist[s].inDest == FALSE ) {
     201            fglmVector v;
     202            if ( dat.monlist[s].basis == 0 ) {
     203                v= fglmVector( dat.basisSize );
     204                // now the monom is in L(source)
     205                PROT( "(" );
     206                poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
     207                PROT( ")" );
     208                poly temp = nf;
     209                while (temp != NULL ) {
     210                    int t;
     211                    for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
     212                        if ( dat.monlist[t].basis > 0 ) {
     213                            if ( pEqual( dat.monlist[t].mon.sm, temp ) ) {
     214                                number coeff= nCopy( pGetCoeff( temp ) );
     215                                v.setelem( dat.monlist[t].basis, coeff );
     216                            }
     217                        }
     218                    }
     219                    pIter(temp);
     220                }
     221                pDelete( & nf );
     222            }
     223            else {
     224                PROT( "." );
     225                v= fglmVector( dat.basisSize, dat.monlist[s].basis );
     226            }
     227            dat.monlist[s].v= v;
     228        }
    229229    }
    230230}
    231231
    232232void
    233 remapVectors( ring oldring, homogData & dat ) 
     233remapVectors( ring oldring, homogData & dat )
    234234{
    235235    nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly );
    236236    int s;
    237237    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
    238         if ( dat.monlist[s].inDest == FALSE ) {
    239             int k;
    240             fglmVector newv( dat.basisSize );
    241             for ( k= dat.basisSize; k > 0; k-- ){
    242                 number newnum= nMap( dat.monlist[s].v.getelem( k ) );
    243                 newv.setelem( k, newnum );
    244             }
    245             dat.monlist[s].dv= newv;
    246         }
     238        if ( dat.monlist[s].inDest == FALSE ) {
     239            int k;
     240            fglmVector newv( dat.basisSize );
     241            for ( k= dat.basisSize; k > 0; k-- ){
     242                number newnum= nMap( dat.monlist[s].v.getelem( k ) );
     243                newv.setelem( k, newnum );
     244            }
     245            dat.monlist[s].dv= newv;
     246        }
    247247    }
    248248}
    249249
    250250void
    251 gaussreduce( homogData & dat, int maxnum, int BS ) 
     251gaussreduce( homogData & dat, int maxnum, int BS )
    252252{
    253253    int s;
    … …  
    256256    int destbasisSize = 0;
    257257    gaussReducer gauss( dat.basisSize );
    258    
     258
    259259    for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
    260         if ( dat.monlist[s].inDest == FALSE ) {
    261             if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
    262                 destbasisSize++;
    263                 dat.monlist[s].destbasis= destbasisSize;
    264                 gauss.store();
    265                 PROT( "." );
    266             }
    267             else {
    268                 fglmVector p= gauss.getDependence();
    269                 poly result = pCopy( dat.monlist[s].mon.dm );
    270                 pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
    271                 int l = 0;
    272                 int k;
    273                 for ( k= 1; k < p.size(); k++ ) {
    274                     if ( ! p.elemIsZero( k ) ) {
    275                         while ( dat.monlist[l].destbasis != k )
    276                             l++;
    277                         poly temp = pCopy( dat.monlist[l].mon.dm );
    278                         pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
    279                         result= pAdd( result, temp );
    280                     }
    281                 }
    282                 if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
    283 //              PROT2( "(%s)", pString( result ) );
    284                 PROT( "+" );
    285                 found++;
    286                 (dat.destIdeal->m)[dat.numDestPolys]= result;
    287                 dat.numDestPolys++;
    288                 if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
    289                     pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
    290                     IDELEMS( dat.destIdeal )+= BS;
    291                 }
    292                
    293             }
    294                
    295         }
     260        if ( dat.monlist[s].inDest == FALSE ) {
     261            if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
     262                destbasisSize++;
     263                dat.monlist[s].destbasis= destbasisSize;
     264                gauss.store();
     265                PROT( "." );
     266            }
     267            else {
     268                fglmVector p= gauss.getDependence();
     269                poly result = pCopy( dat.monlist[s].mon.dm );
     270                pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
     271                int l = 0;
     272                int k;
     273                for ( k= 1; k < p.size(); k++ ) {
     274                    if ( ! p.elemIsZero( k ) ) {
     275                        while ( dat.monlist[l].destbasis != k )
     276                            l++;
     277                        poly temp = pCopy( dat.monlist[l].mon.dm );
     278                        pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
     279                        result= pAdd( result, temp );
     280                    }
     281                }
     282                if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
     283//                PROT2( "(%s)", pString( result ) );
     284                PROT( "+" );
     285                found++;
     286                (dat.destIdeal->m)[dat.numDestPolys]= result;
     287                dat.numDestPolys++;
     288                if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
     289                    pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
     290                    IDELEMS( dat.destIdeal )+= BS;
     291                }
     292
     293            }
     294
     295        }
    296296    }
    297297    PROT2( "(%i", s );
    … …  
    313313    ring sourceRing = currRing;
    314314
    315     intvec * hilb = hHstdSeries( sourceIdeal, NULL, currQuotient ); 
     315    intvec * hilb = hHstdSeries( sourceIdeal, NULL, currQuotient );
    316316    int s;
    317317    dat.sourceIdeal= sourceIdeal;
    318318    dat.sourceHeads= (doublepoly *)Alloc( IDELEMS( sourceIdeal ) * sizeof( doublepoly ) );
    319319    for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) {
    320         dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
     320        dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
    321321    }
    322322    dat.numSourceHeads= IDELEMS( sourceIdeal );
    … …  
    329329    nSetMap( sourceRing->ch, sourceRing->parameter, sourceRing->P, sourceRing->minpoly );
    330330    for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) {
    331         dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
     331        dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
    332332    }
    333333
    … …  
    336336
    337337    while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 ) {
    338         int num = 0;  // the number of monoms of degree deg
    339         PROT2( "deg= %i ", deg );
    340         PROT2( "num= %i\ngen>", numGBelems );
    341         dat.monlistblock= 512;
    342         dat.monlistmax= dat.monlistblock;
    343         dat.monlist= (homogElem *)Alloc( dat.monlistmax*sizeof( homogElem ) );
    344         int j;
    345         for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
    346         dat.numMonoms= 0;
    347         dat.basisSize= 0;
    348         dat.overall= 0;
    349         dat.numberofdestbasismonoms= 0;
    350        
    351         poly start= pOne();
    352         generateMonoms( start, 1, deg, &dat );
    353         pDelete( & start );
    354 
    355         PROT2( "(%i/", dat.basisSize );
    356         PROT2( "%i)\nvec>", dat.overall );
    357         // switch to sourceRing and map monoms
    358         rSetHdl( sourceRingHdl, TRUE );
    359         mapMonoms( destRing, dat );
    360         getVectorRep( dat );
    361        
    362         // switch to destination Ring and remap the vectors
    363         rSetHdl( destRingHdl, TRUE );
    364         remapVectors( sourceRing, dat );
    365        
    366         PROT( "<\nred>" );
    367         // now do gaussian reduction
    368         gaussreduce( dat, numGBelems, groebnerBS );
    369 
    370         Free( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
    371         PROT( "<\n" );
     338        int num = 0;  // the number of monoms of degree deg
     339        PROT2( "deg= %i ", deg );
     340        PROT2( "num= %i\ngen>", numGBelems );
     341        dat.monlistblock= 512;
     342        dat.monlistmax= dat.monlistblock;
     343        dat.monlist= (homogElem *)Alloc( dat.monlistmax*sizeof( homogElem ) );
     344        int j;
     345        for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
     346        dat.numMonoms= 0;
     347        dat.basisSize= 0;
     348        dat.overall= 0;
     349        dat.numberofdestbasismonoms= 0;
     350
     351        poly start= pOne();
     352        generateMonoms( start, 1, deg, &dat );
     353        pDelete( & start );
     354
     355        PROT2( "(%i/", dat.basisSize );
     356        PROT2( "%i)\nvec>", dat.overall );
     357        // switch to sourceRing and map monoms
     358        rSetHdl( sourceRingHdl, TRUE );
     359        mapMonoms( destRing, dat );
     360        getVectorRep( dat );
     361
     362        // switch to destination Ring and remap the vectors
     363        rSetHdl( destRingHdl, TRUE );
     364        remapVectors( sourceRing, dat );
     365
     366        PROT( "<\nred>" );
     367        // now do gaussian reduction
     368        gaussreduce( dat, numGBelems, groebnerBS );
     369
     370        Free( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
     371        PROT( "<\n" );
    372372    }
    373373    PROT( "\n" );
Note: See TracChangeset for help on using the changeset viewer.