Changeset ebecd1b in git


Ignore:
Timestamp:
Oct 4, 2018, 2:40:42 PM (6 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
a0d057c19c0772366f00887ef3f52e956cd70656
Parents:
3a66ff4125e3b112bc7befbe551f30a1342923db3db01637348691ba2d4080e4c14f5e4f134c5be5
Message:
Merge branch 'letterplace_kernel_multiplication' of https://github.com/kabouzeid/Singular into kabouzeid-letterplace_kernel_multiplication
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • Singular/extra.cc

    r3a66ff rebecd1b  
    5757#include "polys/weight.h"
    5858
     59#include "polys/shiftop.h"
     60
    5961#include "coeffs/bigintmat.h"
    6062#include "kernel/fast_mult.h"
    … …  
    6668#include "kernel/GBEngine/kutil.h"
    6769
    68 #include "kernel/GBEngine/shiftgb.h"
    6970#include "kernel/linear_algebra/linearAlgebra.h"
    7071
    … …  
    11811182          return TRUE;
    11821183        }
    1183         res->data = p_LPshift(p,sh,currRing);
     1184        p_LPshift(p,sh,currRing);
     1185        res->data = p;
    11841186        res->rtyp = POLY_CMD;
    11851187        return FALSE;

  • kernel/GBEngine/kspoly.cc

    r3a66ff rebecd1b  
    2121#endif
    2222#include "kernel/GBEngine/shiftgb.h"
     23#include "polys/shiftop.h"
    2324
    2425#ifdef KDEBUG

  • kernel/GBEngine/shiftgb.cc

    r3a66ff rebecd1b  
    3131#include "kernel/GBEngine/shiftgb.h"
    3232#include "polys/nc/sca.h"
    33 
    34 
    35 #define freeT(A,v) omFreeSize((ADDRESS)A,(v+1)*sizeof(int))
    36 
    37 poly p_LPshift(poly p, int sh, const ring r)
    38 {
    39   if (sh == 0 || p == NULL) return(p);
    40 
    41   poly q  = NULL;
    42   poly pp = p;
    43   while (pp!=NULL)
    44   {
    45     poly h=pp;
    46     pIter(pp);
    47     pNext(h)=NULL;
    48     h=p_mLPshift(h,sh,r);
    49     q = p_Add_q(q, h,r);
    50   }
    51   return(q);
    52 }
    53 
    54 poly p_mLPshift(poly p, int sh, const ring r)
    55 {
    56   if (sh == 0 || p == NULL || p_LmIsConstantComp(p,r)) return(p);
    57 
    58   int lV = r->isLPring;
    59 
    60   int L = p_mLastVblock(p,r);
    61   assume(L+sh>=1);
    62   assume(L+sh<=r->N/lV);
    63 
    64   int *e=(int *)omAlloc0((r->N+1)*sizeof(int));
    65   int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
    66   p_GetExpV(p,e,r);
    67 
    68   int j;
    69   //  for (j=1; j<=r->N; j++)
    70   // L*lV gives the last position of the last block
    71   for (j=1; j<= L*lV ; j++)
    72   {
    73     assume(e[j]<=1);
    74     if (e[j]==1)
    75     {
    76       assume(j + (sh*lV)<=r->N);
    77       assume(j + (sh*lV)>=1);
    78       s[j + (sh*lV)] = e[j]; /* actually 1 */
    79     }
    80   }
    81   p_SetExpV(p,s,r);
    82   freeT(e, r->N);
    83   freeT(s, r->N);
    84   /*  pSetm(m); */ /* done in the pSetExpV */
    85   /* think on the component and coefficient */
    86   //  number c = pGetCoeff(p);
    87   //  p_SetCoeff0(m,p_GetCoeff(p,r),r);
    88   return(p);
    89 }
     33#include "polys/shiftop.h"
    9034
    9135poly p_LPCopyAndShiftLM(poly p, int sh, const ring r)
    9236{
    9337  if (sh == 0 || p == NULL) return p;
    94   poly q = p_mLPshift(p_Head(p, r), sh, r);
     38
     39  poly q = p_Head(p, r);
     40  p_mLPshift(q, sh, r);
    9541  pNext(q) = pNext(p);
    9642  return q;
    97 }
    98 
    99 /* returns the number of maximal block */
    100 /* appearing among the monomials of p */
    101 /* the 0th block is the 1st one */
    102 int p_LastVblock(poly p, const ring r)
    103 {
    104   poly q = p;
    105   int ans = 0;
    106   while (q!=NULL)
    107   {
    108     int ansnew = p_mLastVblock(q, r);
    109     ans    = si_max(ans,ansnew);
    110     pIter(q);
    111   }
    112   return(ans);
    113 }
    114 
    115 /* for a monomial p, returns the number of the last block */
    116 /* where a nonzero exponent is sitting */
    117 int p_mLastVblock(poly p, const ring r)
    118 {
    119   if (p == NULL || p_LmIsConstantComp(p,r))
    120   {
    121     return(0);
    122   }
    123 
    124   int lV = r->isLPring;
    125   int *e=(int *)omAlloc0((r->N+1)*sizeof(int));
    126   p_GetExpV(p,e,r);
    127   int j,b;
    128   j = r->N;
    129   while ( (!e[j]) && (j>=1) ) j--;
    130   freeT(e, r->N);
    131   assume(j>0);
    132   b = (int)((j+lV-1)/lV); /* the number of the block, >=1 */
    133   return (b);
    134 }
    135 
    136 /* returns the number of maximal block */
    137 /* appearing among the monomials of p */
    138 /* the 0th block is the 1st one */
    139 int p_FirstVblock(poly p, const ring r)
    140 {
    141   if (p == NULL) {
    142     return 0;
    143   }
    144 
    145   poly q = p;
    146   int ans = p_mFirstVblock(q, r);
    147   while (q!=NULL)
    148   {
    149     int ansnew = p_mFirstVblock(q, r);
    150     if (ansnew > 0) { // don't count constants
    151       ans = si_min(ans,ansnew);
    152     }
    153     pIter(q);
    154   }
    155   /* do not need to delete q */
    156   return(ans);
    157 }
    158 
    159 /* for a monomial p, returns the number of the first block */
    160 /* where a nonzero exponent is sitting */
    161 int p_mFirstVblock(poly p, const ring r)
    162 {
    163   if (p == NULL || p_LmIsConstantComp(p,r))
    164   {
    165     return(0);
    166   }
    167 
    168   int lV = r->isLPring;
    169   int *e=(int *)omAlloc0((r->N+1)*sizeof(int));
    170   p_GetExpV(p,e,r);
    171   int j,b;
    172   j = 1;
    173   while ( (!e[j]) && (j<=r->N-1) ) j++;
    174   freeT(e, r->N);
    175   assume(j <= r->N);
    176   b = (int)(j+lV-1)/lV; /* the number of the block, 1<= b <= r->N  */
    177   return (b);
    17843}
    17944
    … …  
    21277  }
    21378  /* do not need e anymore */
    214   freeT(e, currRing->N);
     79  omFreeSize((ADDRESS) e, (currRing->N+1)*sizeof(int));
    21580
    21681  if (j==0) goto ret_true;
    … …  
    22590    if (B[j]!=1)
    22691    {
    227       freeT(B, b);
     92      omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
    22893      return(0);
    22994    }
    23095  }
    23196 ret_true:
    232   freeT(B, b);
     97  omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
    23398  return(1);
    23499}
    … …  
    297162  m1 = p_GetExp_k_n(m1, hole, r->N, r);
    298163
    299   p_mLPshift(m2, 1 - p_mFirstVblock(m2, r), r);
     164  p_mLPunshift(m2, r);
    300165  p_SetCoeff(m1, m1Coeff, r);
     166
     167  assume(pFirstVblock(m1,r) <= 1);
     168  assume(pFirstVblock(m2,r) <= 1);
    301169}
    302170

  • kernel/GBEngine/shiftgb.h

    r3a66ff rebecd1b  
    1212#include "polys/nc/nc.h"
    1313
    14 poly p_LPshift(poly p, int sh, const ring r);
    15 poly p_mLPshift(poly p, int sh, const ring r);
    1614poly p_LPCopyAndShiftLM(poly p, int sh, const ring r);
    17 
    18 #define pLPshift(p, sh) p_LPshift(p, sh, currRing)
    19 #define pmLPshift(p, sh) p_mLPshift(p, sh, currRing)
    2015#define pLPCopyAndShiftLM(p, sh) p_LPCopyAndShiftLM(p, sh, currRing)
    21 
    22 int p_LastVblock(poly p, const ring r);
    23 int p_mLastVblock(poly p, const ring r);
    24 
    25 #define pLastVblock(p) p_LastVblock(p,currRing)
    26 #define pmLastVblock(p) p_mLastVblock(p,currRing)
    27 
    28 int p_FirstVblock(poly p, const ring r);
    29 int p_mFirstVblock(poly p, const ring r);
    30 
    31 #define pFirstVblock(p) p_FirstVblock(p,currRing)
    32 #define pmFirstVblock(p) p_mFirstVblock(p,currRing)
    3316
    3417int isInV(poly p, const ring r);

  • libpolys/polys/shiftop.cc

    r3a66ff rebecd1b  
    11#include "shiftop.h"
    22#include "templates/p_MemCopy.h"
     3#include "monomials/p_polys.h"
    34
    45/* #define SHIFT_MULT_DEBUG */
    … …  
    3132#ifdef SHIFT_MULT_COMPAT_MODE
    3233  _m = p_Copy(_m, ri);
    33   p_mLPUnShift(_m, ri);
     34  p_mLPunshift(_m, ri);
    3435  p = p_Copy(p, ri);
    3536  poly pCopyHead = p; // used to delete p later
    36   p_LPUnShift(p, ri);
     37  p_LPunshift(p, ri);
    3738#else
    3839  assume(p_mFirstVblock(_m, ri) <= 1);
    … …  
    4849  pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
    4950
    50   int mLength = p_mLastVblock(_m, ri) * lV;
    5151  int *mExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    5252  p_GetExpV(_m,mExpV,ri);
     53  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
    5354  do
    5455  {
    … …  
    5960    int *pExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    6061    p_GetExpV(p, pExpV, ri);
    61     p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, ri) * lV, mLength, ri);
     62    p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
    6263    p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
    6364    p_SetExpV(q, pExpV, ri);
    … …  
    9697#ifdef SHIFT_MULT_COMPAT_MODE
    9798  _m = p_Copy(_m, ri);
    98   p_mLPUnShift(_m, ri);
    99   p_LPUnShift(p, ri);
     99  p_mLPunshift(_m, ri);
     100  p_LPunshift(p, ri);
    100101#else
    101102  assume(p_mFirstVblock(_m, ri) <= 1);
    … …  
    109110  pAssume(!n_IsZero(mCoeff, ri->cf));
    110111
    111   int mLength = p_mLastVblock(_m, ri) * lV;
    112112  int *mExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    113113  p_GetExpV(_m,mExpV,ri);
     114  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
    114115  while (p != NULL)
    115116  {
    … …  
    120121    int *pExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    121122    p_GetExpV(p,pExpV,ri);
    122     p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, ri) * lV, mLength, ri);
     123    p_LPExpVappend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
    123124    p_SetExpV(p, pExpV, ri);
    124125    omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
    … …  
    154155#ifdef SHIFT_MULT_COMPAT_MODE
    155156  _m = p_Copy(_m, ri);
    156   p_mLPUnShift(_m, ri);
     157  p_mLPunshift(_m, ri);
    157158  p = p_Copy(p, ri);
    158159  poly pCopyHead = p; // used to delete p later
    159   p_LPUnShift(p, ri);
     160  p_LPunshift(p, ri);
    160161#else
    161162  assume(p_mFirstVblock(_m, ri) <= 1);
    … …  
    171172  pAssume1(p_GetComp(m, ri) == 0 || p_MaxComp(p, ri) == 0);
    172173
    173   int mLength = p_mLastVblock(_m, ri) * lV;
    174174  int *mExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    175175  p_GetExpV(_m,mExpV,ri);
     176  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
    176177  do
    177178  {
    … …  
    182183    int *pExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    183184    p_GetExpV(p, pExpV, ri);
    184     p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, ri) * lV, mLength, ri);
     185    p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
    185186    p_MemCopy_LengthGeneral(q->exp, p->exp, ri->ExpL_Size); // otherwise q is not initialized correctly
    186187    p_SetExpV(q, pExpV, ri);
    … …  
    219220#ifdef SHIFT_MULT_COMPAT_MODE
    220221  _m = p_Copy(_m, ri);
    221   p_mLPUnShift(_m, ri);
    222   p_LPUnShift(p, ri);
     222  p_mLPunshift(_m, ri);
     223  p_LPunshift(p, ri);
    223224#else
    224225  assume(p_mFirstVblock(_m, ri) <= 1);
    … …  
    232233  pAssume(!n_IsZero(mCoeff, ri->cf));
    233234
    234   int mLength = p_mLastVblock(_m, ri) * lV;
    235235  int *mExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    236236  p_GetExpV(_m,mExpV,ri);
     237  int mLength = p_mLastVblock(_m, mExpV, ri) * lV;
    237238  while (p != NULL)
    238239  {
    … …  
    243244    int *pExpV = (int *) omAlloc0((ri->N+1)*sizeof(int));
    244245    p_GetExpV(p,pExpV,ri);
    245     p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, ri) * lV, mLength, ri);
     246    p_LPExpVprepend(pExpV, mExpV, p_mLastVblock(p, pExpV, ri) * lV, mLength, ri);
    246247    p_SetExpV(p, pExpV, ri);
    247248    omFreeSize((ADDRESS) pExpV, (ri->N+1)*sizeof(int));
    … …  
    311312
    312313// unshifts the monomial m
    313 void p_mLPUnShift(poly m, const ring ri)
     314void p_mLPunshift(poly m, const ring ri)
    314315{
    315316  if (m == NULL || p_LmIsConstantComp(m,ri)) return;
    … …  
    320321
    321322  if (shift == 0) return;
    322 
    323   int L = p_mLastVblock(m, ri);
    324323
    325324  int *e=(int *)omAlloc0((ri->N+1)*sizeof(int));
    … …  
    328327
    329328  int expVoffset = shift*lV;
    330   for (int i = 1 + expVoffset; i <= L*lV; ++i)
     329  for (int i = 1 + expVoffset; i <= ri->N; i++)
    331330  {
    332331    assume(e[i] <= 1);
    … …  
    339338
    340339// unshifts the polynomial p, note: the ordering can be destroyed if the shifts for the monomials are not equal
    341 void p_LPUnShift(poly p, const ring ri)
     340void p_LPunshift(poly p, const ring ri)
    342341{
    343342  while (p!=NULL)
    344343  {
    345     p_mLPUnShift(p, ri);
     344    p_mLPunshift(p, ri);
    346345    pIter(p);
    347346  }
    348   p_Test(p, ri); // check if ordering was destroyed
     347}
     348
     349void p_mLPshift(poly m, int sh, const ring ri)
     350{
     351  if (sh == 0 || m == NULL || p_LmIsConstantComp(m,ri)) return;
     352
     353  int lV = ri->isLPring;
     354
     355  assume(p_mFirstVblock(m,ri) + sh >= 1);
     356  assume(p_mLastVblock(m,ri) + sh <= ri->N/lV);
     357
     358  int *e=(int *)omAlloc0((ri->N+1)*sizeof(int));
     359  int *s=(int *)omAlloc0((ri->N+1)*sizeof(int));
     360  p_GetExpV(m,e,ri);
     361
     362  for (int i = 1; i <= ri->N; i++)
     363  {
     364    assume(e[i]<=1);
     365    if (e[i]==1)
     366    {
     367      assume(i + (sh*lV) <= r->N);
     368      assume(i + (sh*lV) >= 1);
     369      s[i + (sh*lV)] = e[i]; /* actually 1 */
     370    }
     371  }
     372  p_SetExpV(m,s,ri);
     373  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
     374  omFreeSize((ADDRESS) s, (ri->N+1)*sizeof(int));
     375}
     376
     377void p_LPshift(poly p, int sh, const ring ri)
     378{
     379  if (sh == 0) return;
     380
     381  while (p!=NULL)
     382  {
     383    p_mLPshift(p, sh, ri);
     384    pIter(p);
     385  }
     386}
     387
     388/* returns the number of maximal block */
     389/* appearing among the monomials of p */
     390/* the 0th block is the 1st one */
     391int p_LastVblock(poly p, const ring r)
     392{
     393  poly q = p;
     394  int ans = 0;
     395  while (q!=NULL)
     396  {
     397    int ansnew = p_mLastVblock(q, r);
     398    ans    = si_max(ans,ansnew);
     399    pIter(q);
     400  }
     401  return(ans);
     402}
     403
     404/* for a monomial p, returns the number of the last block */
     405/* where a nonzero exponent is sitting */
     406int p_mLastVblock(poly p, const ring ri)
     407{
     408  if (p == NULL || p_LmIsConstantComp(p,ri))
     409  {
     410    return(0);
     411  }
     412
     413  int *e=(int *)omAlloc0((ri->N+1)*sizeof(int));
     414  p_GetExpV(p,e,ri);
     415  int b = p_mLastVblock(p, e, ri);
     416  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
     417  return b;
     418}
     419
     420/* for a monomial p with exponent vector expV, returns the number of the last block */
     421/* where a nonzero exponent is sitting */
     422int p_mLastVblock(poly p, int *expV, const ring ri)
     423{
     424  if (p == NULL || p_LmIsConstantComp(p,ri))
     425  {
     426    return(0);
     427  }
     428
     429  int lV = ri->isLPring;
     430  int j,b;
     431  j = ri->N;
     432  while ( (!expV[j]) && (j>=1) ) j--;
     433  assume(j>0);
     434  b = (int)((j+lV-1)/lV); /* the number of the block, >=1 */
     435  return b;
     436}
     437
     438/* returns the number of maximal block */
     439/* appearing among the monomials of p */
     440/* the 0th block is the 1st one */
     441int p_FirstVblock(poly p, const ring r)
     442{
     443  if (p == NULL) {
     444    return 0;
     445  }
     446
     447  poly q = p;
     448  int ans = p_mFirstVblock(q, r);
     449  while (q!=NULL)
     450  {
     451    int ansnew = p_mFirstVblock(q, r);
     452    if (ansnew > 0) { // don't count constants
     453      ans = si_min(ans,ansnew);
     454    }
     455    pIter(q);
     456  }
     457  /* do not need to delete q */
     458  return(ans);
     459}
     460
     461/* for a monomial p, returns the number of the first block */
     462/* where a nonzero exponent is sitting */
     463int p_mFirstVblock(poly p, const ring ri)
     464{
     465  if (p == NULL || p_LmIsConstantComp(p,ri))
     466  {
     467    return(0);
     468  }
     469
     470  int *e=(int *)omAlloc0((ri->N+1)*sizeof(int));
     471  p_GetExpV(p,e,ri);
     472  int b = p_mFirstVblock(p, e, ri);
     473  omFreeSize((ADDRESS) e, (ri->N+1)*sizeof(int));
     474  return b;
     475}
     476
     477/* for a monomial p with exponent vector expV, returns the number of the first block */
     478/* where a nonzero exponent is sitting */
     479int p_mFirstVblock(poly p, int *expV, const ring ri)
     480{
     481  if (p == NULL || p_LmIsConstantComp(p,ri))
     482  {
     483    return(0);
     484  }
     485
     486  int lV = ri->isLPring;
     487  int j,b;
     488  j = 1;
     489  while ( (!expV[j]) && (j<=ri->N-1) ) j++;
     490  assume(j <= r->N);
     491  b = (int)(j+lV-1)/lV; /* the number of the block, 1<= b <= r->N  */
     492  return b;
    349493}
    350494
    … …  
    356500  PrintLn(); WriteLPExpV(m2ExpV, ri);
    357501#endif
    358   assume(m1Length + m2Length <= ri->N);
     502  assume(m1Length + m2Length <= ri->N); // always throw an error?
    359503  for (int i = 1 + m1Length; i < 1 + m1Length + m2Length; ++i)
    360504  {
    … …  
    377521  PrintLn(); WriteLPExpV(m2ExpV, ri);
    378522#endif
    379   assume(m1Length + m2Length <= ri->N);
     523  assume(m1Length + m2Length <= ri->N); // always throw an error?
    380524
    381525  // shift m1 by m2Length

  • libpolys/polys/shiftop.h

    r3a66ff rebecd1b  
    11#include "misc/auxiliary.h"
    22#include "monomials/ring.h"
    3 #include "../kernel/GBEngine/shiftgb.h"
    43
    54poly shift_pp_Mult_mm(poly p, const poly m, const ring r);
    … …  
    1110poly shift_pp_Mult_Coeff_mm_DivSelectMult_STUB(poly p,const poly m, const poly a, const poly b, int &shorter,const ring r);
    1211poly shift_pp_Mult_Coeff_mm_DivSelect_STUB(poly p, const poly m, int &shorter, const ring r);
    13 void p_mLPUnShift(poly m, const ring ri);
    14 void p_LPUnShift(poly p, const ring ri);
     12
     13void p_mLPunshift(poly m, const ring ri);
     14void p_LPunshift(poly p, const ring ri);
     15void p_mLPshift(poly p, int sh, const ring r);
     16void p_LPshift(poly p, int sh, const ring r);
     17#define pLPunshift(p, sh) p_LPunshift(p, currRing)
     18#define pmLPunshift(p, sh) p_mLPunshift(p, currRing)
     19#define pLPshift(p, sh) p_LPshift(p, sh, currRing)
     20#define pmLPshift(p, sh) p_mLPshift(p, sh, currRing)
     21
     22int p_LastVblock(poly p, const ring r);
     23int p_mLastVblock(poly p, const ring r);
     24int p_mLastVblock(poly p, int *expV, const ring r);
     25int p_FirstVblock(poly p, const ring r);
     26int p_mFirstVblock(poly p, const ring r);
     27int p_mFirstVblock(poly p, int *expV, const ring r);
     28#define pLastVblock(p) p_LastVblock(p,currRing)
     29#define pmLastVblock(p) p_mLastVblock(p,currRing)
     30#define pFirstVblock(p) p_FirstVblock(p,currRing)
     31#define pmFirstVblock(p) p_mFirstVblock(p,currRing)
     32
    1533void p_LPExpVappend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri);
    1634void p_LPExpVprepend(int *m1ExpV, int *m2ExpV, int m1Length, int m2Length, const ring ri);
     35
    1736void WriteLPExpV(int *expV, ring ri);
Note: See TracChangeset for help on using the changeset viewer.