Changeset 0b0bc3 in git


Ignore:
Timestamp:
Feb 8, 2013, 6:03:33 PM (10 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
c462b60ab56fe1781c7f591d40ccd2799acd6c92
Parents:
67242b64d178faf65f7a6059740248dcfbf15c53
Message:
chg: n_Farey, p_Farey, p_ChineseRemainder

from master
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • Singular/iparith.cc

    r67242b r0b0bc3  
    17281728  ideal result;
    17291729  ideal *x=(ideal *)omAlloc(rl*sizeof(ideal));
     1730  number *xx=NULL;
    17301731  int i;
    17311732  int return_type=c->m[0].Typ();
    … …  
    17341735  && (return_type!=MATRIX_CMD))
    17351736  {
    1736     WerrorS("ideal/module/matrix expected");
    1737     omFree(x); // delete c
    1738     return TRUE;
    1739   }
    1740   for(i=rl-1;i>=0;i--)
    1741   {
    1742     if (c->m[i].Typ()!=return_type)
    1743     {
    1744       Werror("%s expected at pos %d",Tok2Cmdname(return_type),i+1);
     1737    if((return_type!=BIGINT_CMD)&&(return_type!=INT_CMD))
     1738    {
     1739      WerrorS("ideal/module/matrix expected");
    17451740      omFree(x); // delete c
    17461741      return TRUE;
    17471742    }
    1748     x[i]=((ideal)c->m[i].Data());
     1743    else
     1744      return_type=BIGINT_CMD;
     1745  }
     1746  if (return_type!=BIGINT_CMD)
     1747  {
     1748    for(i=rl-1;i>=0;i--)
     1749    {
     1750      if (c->m[i].Typ()!=return_type)
     1751      {
     1752        Werror("%s expected at pos %d",Tok2Cmdname(return_type),i+1);
     1753        omFree(x); // delete c
     1754        return TRUE;
     1755      }
     1756      x[i]=((ideal)c->m[i].Data());
     1757    }
     1758  }
     1759  else
     1760  {
     1761    xx=(number *)omAlloc(rl*sizeof(number));
     1762    for(i=rl-1;i>=0;i--)
     1763    {
     1764      if (c->m[i].Typ()==INT_CMD)
     1765      {
     1766        xx[i]=n_Init(((int)(long)c->m[i].Data()),coeffs_BIGINT);
     1767      }
     1768      else if (c->m[i].Typ()==BIGINT_CMD)
     1769      {
     1770        xx[i]=(number)c->m[i].Data();
     1771      }
     1772      else
     1773      {
     1774        Werror("bigint expected at pos %d",i+1);
     1775        omFree(x); // delete c
     1776        omFree(xx); // delete c
     1777        return TRUE;
     1778      }
     1779    }
    17491780  }
    17501781  number *q=(number *)omAlloc(rl*sizeof(number));
    … …  
    17621793      if (pl->m[i].Typ()==INT_CMD)
    17631794      {
    1764         q[i]=n_Init((int)(long)pl->m[i].Data(),currRing->cf);
     1795        if (return_type==BIGINT_CMD)
     1796          q[i]=n_Init((int)(long)pl->m[i].Data(),coeffs_BIGINT);
     1797        else
     1798          q[i]=n_Init((int)(long)pl->m[i].Data(),currRing->cf);
    17651799      }
    17661800      else if (pl->m[i].Typ()==BIGINT_CMD)
    17671801      {
    1768         q[i]=n_Init_bigint((number)(pl->m[i].Data()),coeffs_BIGINT,currRing->cf);
     1802        if (return_type==BIGINT_CMD)
     1803          q[i]=n_Copy((number)(pl->m[i].Data()),coeffs_BIGINT);
     1804        else
     1805          q[i]=n_Init_bigint((number)(pl->m[i].Data()),coeffs_BIGINT,currRing->cf);
    17691806      }
    17701807      else
    17711808      {
    17721809        Werror("bigint expected at pos %d",i+1);
     1810        if (return_type==BIGINT_CMD)
     1811        for(i++;i<rl;i++)
     1812        {
     1813          n_Delete(&(q[i]),coeffs_BIGINT);
     1814        }
     1815        else
    17731816        for(i++;i<rl;i++)
    17741817        {
    … …  
    17811824    }
    17821825  }
    1783   result=id_ChineseRemainder(x,q,rl,currRing);
     1826  if (return_type==BIGINT_CMD)
     1827  {
     1828    number n=n_ChineseRemainder(xx,q,rl,coeffs_BIGINT);
     1829    res->data=(char *)n;
     1830  }
     1831  else
     1832  {
     1833    result=id_ChineseRemainder(x,q,rl,currRing);
     1834    // deletes also x
     1835    res->data=(char *)result;
     1836  }
     1837  if (return_type==BIGINT_CMD)
    17841838  for(i=rl-1;i>=0;i--)
    17851839  {
     1840    n_Delete(&(q[i]),coeffs_BIGINT);
     1841  }
     1842  else
     1843  for(i=rl-1;i>=0;i--)
     1844  {
    17861845    n_Delete(&(q[i]),currRing->cf);
    17871846  }
    17881847  omFree(q);
    1789   res->data=(char *)result;
    17901848  res->rtyp=return_type;
    17911849  return FALSE;
    … …  
    21432201static BOOLEAN jjFAREY_ID(leftv res, leftv u, leftv v)
    21442202{
    2145   if (rField_is_Q(currRing))
    2146   {
    2147     ideal uu=(ideal)u->Data();
    2148     number vv=(number)v->Data();
    2149     res->data=(void*)id_Farey(uu,vv,currRing);
    2150     res->rtyp=u->Typ();
    2151     return FALSE;
    2152   }
    2153   else return TRUE;
     2203  ideal uu=(ideal)u->Data();
     2204  number vv=(number)v->Data();
     2205  res->data=(void*)id_Farey(uu,vv,currRing);
     2206  res->rtyp=u->Typ();
     2207  return FALSE;
    21542208}
    21552209static BOOLEAN jjFETCH(leftv res, leftv u, leftv v)

  • kernel/ideals.cc

    r67242b r0b0bc3  
    24002400#endif
    24012401
     2402#ifdef HAVE_FACTORY
    24022403/*2
    24032404* xx,q: arrays of length 0..rl-1
    … …  
    24072408* destroys xx
    24082409*/
    2409 #ifdef HAVE_FACTORY
    24102410ideal id_ChineseRemainder(ideal *xx, number *q, int rl, const ring R)
    24112411{
    … …  
    24152415  result->ncols=xx[0]->ncols; // for lifting matrices
    24162416  int i,j;
    2417   poly r,h,hh,res_p;
    24182417  number *x=(number *)omAlloc(rl*sizeof(number));
     2418  poly *p=(poly *)omAlloc(rl*sizeof(poly));
    24192419  for(i=cnt-1;i>=0;i--)
    24202420  {
    2421     res_p=NULL;
    2422     loop
    2423     {
    2424       r=NULL;
    2425       for(j=rl-1;j>=0;j--)
    2426       {
    2427         h=xx[j]->m[i];
    2428         if ((h!=NULL)
    2429         &&((r==NULL)||(p_LmCmp(r,h,R)==-1)))
    2430           r=h;
    2431       }
    2432       if (r==NULL) break;
    2433       h=p_Head(r,R);
    2434       for(j=rl-1;j>=0;j--)
    2435       {
    2436         hh=xx[j]->m[i];
    2437         if ((hh!=NULL) && (p_LmCmp(r,hh,R)==0))
    2438         {
    2439           x[j]=pGetCoeff(hh);
    2440           hh=p_LmFreeAndNext(hh,R);
    2441           xx[j]->m[i]=hh;
    2442         }
    2443         else
    2444           x[j]=n_Init(0, R->cf);
    2445       }
    2446       number n=n_ChineseRemainder(x,q,rl,R->cf);
    2447       for(j=rl-1;j>=0;j--)
    2448       {
    2449         x[j]=NULL; // nlInit(0...) takes no memory
    2450       }
    2451       if (n_IsZero(n,R->cf)) p_Delete(&h,R);
    2452       else
    2453       {
    2454         p_SetCoeff(h,n,R);
    2455         //Print("new mon:");pWrite(h);
    2456         res_p=p_Add_q(res_p,h,R);
    2457       }
    2458     }
    2459     result->m[i]=res_p;
    2460   }
    2461   omFree(x);
     2421    for(j=rl-1;j>=0;j--)
     2422    {
     2423      p[j]=xx[j]->m[i];
     2424    }
     2425    result->m[i]=p_ChineseRemainder(p,x,q,rl,R);
     2426    for(j=rl-1;j>=0;j--)
     2427    {
     2428      xx[j]->m[i]=p[j];
     2429    }
     2430  }
     2431  omFreeSize(p,rl*sizeof(poly));
     2432  omFreeSize(x,rl*sizeof(number));
    24622433  for(i=rl-1;i>=0;i--) id_Delete(&(xx[i]),R);
    2463   omFree(xx);
     2434  omFreeSize(xx,rl*sizeof(ideal));
    24642435  return result;
    24652436}
    … …  
    25612532  for(i=cnt-1;i>=0;i--)
    25622533  {
    2563     poly h=p_Copy(x->m[i],r);
    2564     result->m[i]=h;
    2565     while(h!=NULL)
    2566     {
    2567       number c=pGetCoeff(h);
    2568       pSetCoeff0(h,n_Farey(c,N,r->cf));
    2569       n_Delete(&c,r->cf);
    2570       pIter(h);
    2571     }
    2572     while((result->m[i]!=NULL)&&(n_IsZero(pGetCoeff(result->m[i]),r->cf)))
    2573     {
    2574       p_LmDelete(&(result->m[i]),r);
    2575     }
    2576     h=result->m[i];
    2577     while((h!=NULL) && (pNext(h)!=NULL))
    2578     {
    2579       if(n_IsZero(pGetCoeff(pNext(h)),r->cf))
    2580       {
    2581         p_LmDelete(&pNext(h),r);
    2582       }
    2583       else pIter(h);
    2584     }
     2534    result->m[i]=p_Farey(x->m[i],N,r);
    25852535  }
    25862536  return result;

  • libpolys/polys/monomials/p_polys.cc

    r67242b r0b0bc3  
    5959#endif
    6060
     61/*
     62 * lift ideal with coeffs over Z (mod N) to Q via Farey
     63 */
     64poly p_Farey(poly p, number N, const ring r)
     65{
     66  poly h=p_Copy(p,r);
     67  poly hh=h;
     68  while(h!=NULL)
     69  {
     70    number c=pGetCoeff(h);
     71    pSetCoeff0(h,n_Farey(c,N,r->cf));
     72    n_Delete(&c,r->cf);
     73    pIter(h);
     74  }
     75  while((hh!=NULL)&&(n_IsZero(pGetCoeff(hh),r->cf)))
     76  {
     77    p_LmDelete(&hh,r);
     78  }
     79  h=hh;
     80  while((h!=NULL) && (pNext(h)!=NULL))
     81  {
     82    if(n_IsZero(pGetCoeff(pNext(h)),r->cf))
     83    {
     84      p_LmDelete(&pNext(h),r);
     85    }
     86    else pIter(h);
     87  }
     88  return hh;
     89}
     90/*2
     91* xx,q: arrays of length 0..rl-1
     92* xx[i]: SB mod q[i]
     93* assume: char=0
     94* assume: q[i]!=0
     95* destroys xx
     96*/
     97poly p_ChineseRemainder(poly *xx, number *x,number *q, int rl, const ring R)
     98{
     99  poly r,h,hh;
     100  int j;
     101  poly  res_p=NULL;
     102  loop
     103  {
     104    /* search the lead term */
     105    r=NULL;
     106    for(j=rl-1;j>=0;j--)
     107    {
     108      h=xx[j];
     109      if ((h!=NULL)
     110      &&((r==NULL)||(p_LmCmp(r,h,R)==-1)))
     111        r=h;
     112    }
     113    /* nothing found -> return */
     114    if (r==NULL) break;
     115    /* create the monomial in h */
     116    h=p_Head(r,R);
     117    /* collect the coeffs in x[..]*/
     118    for(j=rl-1;j>=0;j--)
     119    {
     120      hh=xx[j];
     121      if ((hh!=NULL) && (p_LmCmp(r,hh,R)==0))
     122      {
     123        x[j]=pGetCoeff(hh);
     124        hh=p_LmFreeAndNext(hh,R);
     125        xx[j]=hh;
     126      }
     127      else
     128        x[j]=n_Init(0, R);
     129    }
     130    number n=n_ChineseRemainder(x,q,rl,R->cf);
     131    for(j=rl-1;j>=0;j--)
     132    {
     133      x[j]=NULL; // nlInit(0...) takes no memory
     134    }
     135    if (n_IsZero(n,R)) p_Delete(&h,R);
     136    else
     137    {
     138      p_SetCoeff(h,n,R);
     139      //Print("new mon:");pWrite(h);
     140      pNext(h)=res_p;
     141      res_p=h; // building res_p in reverse order!
     142    }
     143  }
     144  return pReverse(res_p);
     145}
    61146/***************************************************************
    62147 *

  • libpolys/polys/monomials/p_polys.h

    r67242b r0b0bc3  
    3636#endif
    3737
     38poly p_Farey(poly p, number N, const ring r);
     39/*
     40* xx,q: arrays of length 0..rl-1
     41* xx[i]: SB mod q[i]
     42* assume: char=0
     43* assume: q[i]!=0
     44* destroys xx
     45*/
     46poly p_ChineseRemainder(poly *xx, number *x,number *q, int rl, const ring R);
    3847/***************************************************************
    3948 *
Note: See TracChangeset for help on using the changeset viewer.