Changeset e212bd in git


Ignore:
Timestamp:
Nov 20, 2013, 4:07:38 PM (9 years ago)
Author:
Adrian <adi_popescum@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '48f1dd268d0ff74ef2f7dccbf02545425002ddcc')
Children:
fd1101c51f85808be5f4f270758360b46a425a4f
Parents:
867d59cf46048ef5ad1c520db05b0b6d5f71d305
Message:
First Version: implemented some a new function that computes hilb, but there are some errors
File:
1 edited

Legend:

Unmodified
Added
Removed
  • kernel/hilb.cc

    r867d59 re212bd  
    1717#include <kernel/hutil.h>
    1818#include <kernel/stairc.h>
     19//ADICHANGES:
     20#include <kernel/ideals.h>
     21#include <kernel/kstd1.h>
     22#include<gmp.h>
     23#include<vector>
     24
    1925
    2026static int  **Qpol;
    2127static int  *Q0, *Ql;
    2228static int  hLength;
     29//ADICHANGES
     30static int NNN;
     31static mpz_t ec;
     32static mpz_ptr hilbertcoef = (mpz_ptr)omAlloc(NNN*sizeof(mpz_t));
     33std::vector<int> hilbertpowers;
     34std::vector<int> hilbertpowerssorted;
     35   
    2336
    2437static int hMinModulweight(intvec *modulweight)
     
    213226  }
    214227}
    215 
     228// ---------------------------------- ADICHANGES ---------------------------------------------
     229//!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
     230void PutInVector(int degvalue, int oldposition, int where)
     231{
     232    hilbertpowerssorted.resize(hilbertpowerssorted.size()+2);
     233    int i;
     234    for(i=hilbertpowerssorted.size()-1;i>=where+2;i--)
     235    {
     236        hilbertpowerssorted[i] = hilbertpowerssorted[i-2];
     237    }
     238    hilbertpowerssorted[where] = degvalue;
     239    hilbertpowerssorted[where+1] = oldposition;
     240}
     241
     242void SortPowerVec()
     243{
     244    int i,j;
     245    int test;
     246    bool flag;
     247    hilbertpowerssorted.resize(2);
     248    hilbertpowerssorted[0] = hilbertpowers[0];
     249    hilbertpowerssorted[1] = 0;
     250    for(i=1;i<hilbertpowers.size();i++)
     251    {
     252        flag=TRUE;
     253        if(hilbertpowers[i]<=hilbertpowerssorted[0])
     254        {
     255            flag=FALSE;
     256            PutInVector(hilbertpowers[i], i, 0);
     257        }
     258        if(hilbertpowers[i]>=hilbertpowerssorted[hilbertpowerssorted.size()-2])
     259        {
     260            flag=FALSE;
     261            PutInVector(hilbertpowers[i], i, hilbertpowerssorted.size()-2);
     262        }
     263        if(flag==TRUE)
     264        {
     265            for(j=2;(j<=hilbertpowerssorted.size()-4)&&(flag);j=j+2)
     266            {
     267                if(hilbertpowers[i]<=hilbertpowerssorted[j])
     268                {
     269                    flag=FALSE;
     270                    PutInVector(hilbertpowers[i], i, j);
     271                }
     272            }
     273        }
     274        //printf("\n----------------TEST----------------\n");
     275        //for(test=0;test<hilbertpowerssorted.size();test++)
     276        //    {
     277        //        printf(" %d ", hilbertpowerssorted[test]);
     278        //    }   
     279       
     280    }
     281}
     282
     283int DegMon(poly p, ring tailRing)
     284{
     285    int i,deg;
     286    deg = 0;
     287    for(i=0;i<=tailRing->N;i++)
     288    {
     289        deg = deg + p_GetExp(p, i, tailRing);
     290    }
     291    return(deg);
     292}
     293
     294poly SearchP(ideal I, ring tailRing)
     295{
     296    int i,j,exp;
     297    poly res;
     298    for(i=IDELEMS(I)-1;i>=0;i--)
     299    {
     300        if(DegMon(I->m[i], tailRing)>=2)
     301        {
     302            res = p_ISet(1, tailRing);
     303            res = p_Copy(I->m[i], tailRing);
     304            for(j=0;j<=tailRing->N;j++)
     305            {
     306                exp = p_GetExp(I->m[i], j, tailRing);
     307                if(exp > 0)
     308                {
     309                    p_SetExp(res, j, exp - 1, tailRing);
     310                    return(res);
     311                }
     312            }
     313        }
     314    }
     315    return(NULL);
     316}
     317
     318poly ChoosePVar (ideal I, ring tailRing)
     319{
     320    bool flag=TRUE;
     321    int i,j;
     322    poly res;
     323    for(i=1;i<=tailRing->N;i++)
     324    {
     325        flag=TRUE;
     326        for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
     327        {
     328            if(p_GetExp(I->m[j], i, tailRing)>0)
     329            {
     330                flag=FALSE;
     331            }
     332        }
     333       
     334        if(flag == TRUE)
     335        {
     336            res = p_ISet(1, tailRing);
     337            p_SetExp(res, i, 1, tailRing);
     338            //idPrint(I);
     339            //pWrite(res);
     340            return(res);
     341        }
     342    }
     343    return(NULL); //i.e. it is the maximal ideal
     344}
     345
     346ideal idSimplify(ideal I, ring tailRing)
     347{
     348    int i,j;
     349    bool flag;
     350    /*std::vector<int>  var;
     351    for(i=0;i<=tailRing->N;i++)
     352    {
     353        if(p_GetExp(I->[IDELEMS(I)-1], tailRing)>0)
     354        {
     355            var.resize(var.size()+1);
     356            var[var.size()-1] = i;
     357        }
     358    }
     359    */
     360    for(i=IDELEMS(I)-2;i>=0;i--)
     361    {
     362        flag=TRUE;
     363        for(j=0;(j<=tailRing->N)&&(flag);j++)
     364        {
     365            if(p_GetExp(I->m[i], j, tailRing) < p_GetExp(I->m[IDELEMS(I)-1], j, tailRing))
     366            {
     367                flag=FALSE;
     368            }
     369        }
     370        if(flag)
     371        {
     372            I->m[i]=NULL;
     373        }
     374    }
     375    idSkipZeroes(I);
     376    return(I);
     377}
     378
     379
     380void eulerchar (ideal I, ring tailRing)
     381{
     382    //gmp_printf("\nEuler char: %Zd\n", ec);
     383    mpz_t dummy;
     384    if( idElem(I) == 0)
     385    {
     386        mpz_init(dummy);
     387        //change: era 1 in loc de 0
     388        mpz_set_si(dummy, 1);
     389        mpz_sub(ec, ec, dummy);
     390        return;
     391    }
     392    if( idElem(I) == 1)
     393    {
     394        if(!p_IsOne(I->m[0], tailRing))
     395        {
     396            //aici nu era nimic
     397            //mpz_set_si(dummy, 1);
     398            //mpz_add(ec, ec, dummy);
     399            return;
     400        }
     401        else
     402        {
     403            mpz_init(dummy);
     404            mpz_set_si(dummy, 1);
     405            mpz_sub(ec, ec, dummy);
     406            return;
     407        }
     408    }
     409    ideal p = idInit(1,1);
     410    p->m[0] = SearchP(I, tailRing);
     411    //idPrint(I);
     412    //printf("\nSearchP founded: \n");pWrite(p->m[0]);
     413    if(p->m[0] == NULL)
     414    {
     415        mpz_init(dummy);
     416        mpz_set_si(dummy, IDELEMS(I)-1);
     417        mpz_add(ec, ec, dummy);
     418        return;
     419    }
     420    ideal Ip = idQuot(I,p,TRUE,TRUE);
     421    ideal Iplusp = id_Add(I,p,tailRing);
     422    Iplusp=idSimplify(Iplusp, tailRing);
     423    //mpz_t i,j;
     424    //i = eulerchar(Ip, ec, tailRing);
     425    //j = eulerchar(Iplusp, ec, tailRing);
     426    //mpz_add(ec, i,j);
     427    eulerchar(Ip, tailRing);
     428    eulerchar(Iplusp, tailRing);
     429    return;
     430}
     431
     432poly SqFree (ideal I, ring tailRing)
     433{
     434    int i,j;
     435    bool flag=TRUE;
     436    poly notsqrfree = NULL;
     437    for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
     438    {
     439        for(j=0;(j<=tailRing->N)&&(flag);j++)
     440        {
     441            if(p_GetExp(I->m[i],j,tailRing)>1)
     442            {
     443                flag=FALSE;
     444                notsqrfree = p_ISet(1,tailRing);
     445                p_SetExp(notsqrfree,j,1,tailRing);
     446            }
     447        }
     448    }
     449    return(notsqrfree);
     450}
     451
     452void rouneslice(ideal I, ideal S, poly q, ring tailRing, poly x)
     453{
     454    int i,j;
     455    int dummy;
     456    mpz_t dummympz;
     457    poly m;
     458    ideal p, koszsimp;
     459    I = idMinBase(I);
     460   
     461    //Work on it
     462    //----------- PRUNING OF S ---------------
     463    S = idMinBase(S);
     464    for(i=IDELEMS(S)-1;i>=0;i--)
     465    {
     466        if(kNF(I,NULL,S->m[i],NULL,NULL)==NULL)
     467        {
     468            S->m[i]=NULL;
     469        }
     470    }
     471    idSkipZeroes(S);
     472    //----------------------------------------
     473    for(i=IDELEMS(I)-1;i>=0;i--)
     474    {
     475        m = p_Copy(I->m[i],tailRing);
     476        for(j=0;j<=tailRing->N;j++)
     477        {
     478            dummy = p_GetExp(m,j,tailRing);
     479            if(dummy > 0)
     480            {
     481                p_SetExp(m,j,dummy-1,tailRing);
     482            }       
     483        }
     484        if(kNF(S,NULL,m,NULL,NULL)==NULL)
     485        {
     486            I->m[i]=NULL;
     487            //printf("\n Check if m is in S: \n");pWrite(m);idPrint(S);
     488            //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
     489        }
     490    }
     491    idSkipZeroes(I);
     492   
     493    //----------- MORE PRUNING OF S ------------
     494    //strictly divide???   
     495    //------------------------------------------
     496
     497    //printf("Ideal I: \n");idPrint(I);
     498    //printf("Ideal S: \n");idPrint(S);
     499    //printf("Poly  q: \n");pWrite(q);
     500    if(idElem(I)==0)
     501    {
     502        //I = 0
     503        //printf("\nx does not divide lcm(I)");
     504        //pWrite(x);pWrite(m);idPrint(I);
     505        printf("\nEmpty Set: ");pWrite(q);
     506        return;
     507    }
     508    m = p_ISet(1,tailRing);
     509    for(i=0;i<=tailRing->N;i++)
     510    {
     511        dummy=0;
     512        for(j=IDELEMS(I)-1;j>=0;j--)
     513        {
     514            if(p_GetExp(I->m[j],i,tailRing) > dummy)
     515            {
     516                dummy = p_GetExp(I->m[j],i,tailRing);
     517            }
     518        }
     519        p_SetExp(m,i,dummy,tailRing);
     520    }
     521    p_Setm(m,tailRing);
     522    if(p_DivisibleBy(x,m,tailRing)==FALSE)
     523    {
     524        //printf("\nx does not divide lcm(I)");
     525        //pWrite(x);pWrite(m);idPrint(I);
     526        //printf("\nEmpty set");pWrite(q);
     527        return;
     528    }
     529    m = SqFree(I, tailRing);
     530    koszsimp = idInit(IDELEMS(I),1);
     531    if(m==NULL)
     532    {
     533        printf("\n      Corner: ");
     534        pWrite(q);
     535        printf("\n      With the facets of the simplex:\n");
     536        for(i=IDELEMS(I)-1;i>=0;i--)
     537        {
     538            m = p_ISet(1,tailRing);
     539            //m = p_Divide(x,I->m[i],tailRing);
     540            for(j=tailRing->N; j>=0; j--)
     541            {
     542                p_SetExp(m,j, p_GetExp(x,j,tailRing)- p_GetExp(I->m[i],j,tailRing),tailRing);
     543            }
     544            printf("    ");
     545            p_Setm(m, tailRing);
     546            p_Write(m, tailRing);
     547            koszsimp->m[i] = p_Copy(m, tailRing);
     548        }
     549        printf("\n Euler Characteristic = ");
     550        mpz_init(dummympz);
     551        eulerchar(koszsimp, tailRing);
     552        //gmp_printf("dummy %Zd \n", dummympz);
     553        gmp_printf("ec %Zd \n", ec);
     554        mpz_set(dummympz, ec);
     555        mpz_sub(ec, ec, ec);
     556        hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
     557        mpz_init(&hilbertcoef[NNN]);
     558        mpz_set(&hilbertcoef[NNN], dummympz);
     559        hilbertpowers.resize(NNN+1);
     560        hilbertpowers[NNN] = DegMon(q, tailRing);
     561        NNN++;
     562        return;
     563    }
     564    if(IDELEMS(S)!=1)
     565    {
     566        m = SearchP(S, tailRing);
     567        if(m == NULL)
     568        {
     569            m = ChoosePVar(S, tailRing);
     570        }
     571    }
     572    p = idInit(1,1);
     573    p->m[0] = m;
     574    printf("Poly  p: \n");pWrite(m);
     575    ideal Ip = idQuot(I,p,TRUE,TRUE);
     576    ideal Sp = idQuot(S,p,TRUE,TRUE);
     577    ideal Splusp = id_Add(S,p,tailRing);
     578    Splusp = idSimplify(Splusp, tailRing);
     579    poly pq = pp_Mult_mm(q,m,tailRing);
     580    rouneslice(Ip, Sp, pq, tailRing, x);
     581    rouneslice(I, Splusp, q, tailRing, x);
     582    return;
     583}
     584
     585void slicehilb(ideal I, ring tailRing)
     586{
     587    printf("Adi changes are here: \n");
     588    int i;
     589    ideal S = idInit(0,1);
     590    poly q = p_ISet(1,tailRing);
     591    ideal X = idInit(1,1);
     592    X->m[0]=p_One(tailRing);
     593    for(i=1;i<=tailRing->N;i++)
     594    {
     595            p_SetExp(X->m[0],i,1,tailRing);   
     596    }
     597    p_Setm(X->m[0],tailRing);
     598    I = id_Mult(I,X,tailRing);
     599    printf("\n-------------RouneSlice--------------\n");
     600    rouneslice(I,S,q,tailRing,X->m[0]);
     601    for(i=0;i<NNN;i++)
     602    {
     603        gmp_printf(" %Zd ", &hilbertcoef[i]);
     604    }
     605    printf("\n");
     606    for(i=0;i<NNN;i++)
     607    {
     608        printf(" %d ", hilbertpowers[i]);
     609    }
     610    SortPowerVec();
     611    printf("\n");
     612    for(i=0;i<hilbertpowerssorted.size();i++)
     613    {
     614        printf(" %d ", hilbertpowerssorted[i]);
     615    }
     616    mpz_t coefhilb;
     617    mpz_t dummy;
     618    mpz_init(coefhilb);
     619    mpz_init(dummy);
     620    printf("\n//        1 t^0");
     621    for(i=0;i<hilbertpowerssorted.size();i=i+2)
     622    {
     623        mpz_set(dummy, &hilbertcoef[hilbertpowerssorted[i+1]]);
     624        mpz_add(coefhilb, coefhilb, dummy);
     625        while((hilbertpowerssorted[i]==hilbertpowerssorted[i+2]) && (i<=hilbertpowerssorted.size()-2))
     626        {
     627            i=i+2;
     628            mpz_add(coefhilb, coefhilb, &hilbertcoef[hilbertpowerssorted[i+1]]);
     629        }
     630        if(mpz_sgn(coefhilb)!=0)
     631        {
     632            gmp_printf("\n//        %Zd t^%d",coefhilb,hilbertpowerssorted[i]);
     633        }
     634        mpz_sub(coefhilb, coefhilb, coefhilb);
     635    }
     636    omFreeSize(hilbertcoef, NNN*sizeof(mpz_t));
     637    printf("\n-------------------------------------\n");
     638}
     639
     640// -------------------------------- END OF CHANGES -------------------------------------------
    216641static intvec * hSeries(ideal S, intvec *modulweight,
    217642                int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
    218643{
     644  slicehilb(S,tailRing); // ADICHANGES
    219645  intvec *work, *hseries1=NULL;
    220646  int  mc;
     
    472898}
    473899
     900
     901
Note: See TracChangeset for help on using the changeset viewer.