Changeset 58603f in git


Ignore:
Timestamp:
Aug 4, 2019, 9:34:39 PM (4 years ago)
Author:
Murray Heymann <heymann.murray@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
Children:
2f9e3f171a90f561e7a3c81734eb72998dc8fe82
Parents:
112c79d4f91036791f523bcfbb261d1ba63c2d8f557b878dd8c840afb5ec95de122ba27a0ec99926
Message:
Merge branch 'spielwiese' into machine_learning
Files:
17 added
22 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    r557b87 r58603f  
    117117ylwrap
    118118/resources/singular_resourcesconfig.h
     119__pycache__
     120keywords.txt
     121helpfiles
     122*.npy
     123.coverage
     124*.pyc

  • Singular/LIB/ellipticcovers.lib

    r112c79 r58603f  
    6767
    6868
    69 proc mod_init()
    70 {
    71 newstruct("graph","list vertices, list edges");
    72 newstruct("Net","list rows");
    73 
    74 system("install","graph","print",printGraph,1);
    75 system("install","Net","print",printNet,1);
    76 system("install","Net","+",catNet,2);
    77 
    78 }
    79 
     69static proc mod_init()
     70{
     71  newstruct("graph","list vertices, list edges");
     72  newstruct("Net","list rows");
     73
     74  system("install","graph","print",printGraph,1);
     75  system("install","Net","print",printNet,1);
     76  system("install","Net","+",catNet,2);
     77}
    8078
    8179static proc catNet(Net N, Net M)
    8280{
    83 list L;
    84 list LN=N.rows;
    85 list LM=M.rows;
    86 int widthN=size(LN[1]);
    87 int widthM=size(LM[1]);
    88 int nm=max(size(LN),size(LM));
    89 for (int j=1; j<=nm; j++)
    90 {
     81  list L;
     82  list LN=N.rows;
     83  list LM=M.rows;
     84  int widthN=size(LN[1]);
     85  int widthM=size(LM[1]);
     86  int nm=max(size(LN),size(LM));
     87  for (int j=1; j<=nm; j++)
     88  {
    9189    if (j>size(LN)){LN[j]=emptyString(widthN);}
    9290    if (j>size(LM)){LM[j]=emptyString(widthM);}
    9391    L[j]=LN[j]+LM[j];
    94 }
    95 Net NM;
    96 NM.rows=L;
    97 return(NM);}
     92  }
     93  Net NM;
     94  NM.rows=L;
     95  return(NM);
     96}
    9897
    9998
    … …  
    111110static proc printNet(Net N)
    112111{
    113 list L = N.rows;
    114 for (int j=1; j<=size(L); j++)
    115 {
    116    print(L[j]);
    117 }
    118 }
    119 
    120 static proc net(def M){
    121   if (typeof(M)=="list"){
     112  list L = N.rows;
     113  for (int j=1; j<=size(L); j++)
     114  {
     115    print(L[j]);
     116  }
     117}
     118
     119static proc net(def M)
     120{
     121  if (typeof(M)=="list")
     122  {
    122123    return(netList(M));
    123124  }
    … …  
    126127  L[1]=string(M);
    127128  N.rows=L;
    128 return(N);}
     129  return(N);
     130}
    129131
    130132
    … …  
    147149  G;
    148150}
    149 
    150 
    151151
    152152proc makeGraph(list v, list e)
    … …  
    170170  G;
    171171}
    172 
    173172
    174173proc propagator(def xy, def d)
    … …  
    195194"
    196195{
    197   if ((typeof(xy)=="list")||(typeof(d)=="int")) {
     196  if ((typeof(xy)=="list")||(typeof(d)=="int"))
     197  {
    198198    number x = xy[1];
    199199    number y = xy[2];
    … …  
    201201    if (d==0) {return(x^2*y^2/(x^2-y^2)^2);}
    202202    number p=0;
    203     for (int j=1; j<=d; j++){
     203    for (int j=1; j<=d; j++)
     204    {
    204205       if (d%j==0){p=p+(j*x^(4*j)+j*y^(4*j))/(x*y)^(2*j);}
    205206    }
    206207    return(p);
    207208  }
    208   if ((typeof(xy)=="graph")||(typeof(d)=="list"))  {
     209  if ((typeof(xy)=="graph")||(typeof(d)=="list"))
     210  {
    209211    list xl = ringlist(basering)[1][2];
    210212    list ed = xy.edges;
    211213    number f=1;
    212     for (int j=1; j<=size(ed); j++){
     214    for (int j=1; j<=size(ed); j++)
     215    {
    213216       execute("number xx1 = "+xl[ed[j][1]]);
    214217       execute("number xx2 = "+xl[ed[j][2]]);
    … …  
    219222    return(f);
    220223  }
    221   if ((typeof(xy)=="graph")||(typeof(d)=="int"))  {
    222   }
    223 ERROR("wrong input type");}
     224  if ((typeof(xy)=="graph")||(typeof(d)=="int"))
     225  {
     226  }
     227  ERROR("wrong input type");
     228}
    224229example
    225230{ "EXAMPLE:"; echo=2;
    … …  
    230235  propagator(G,list(1,1,1,0,0,0));
    231236}
    232 
    233 
    234 
    235 
    236237
    237238proc computeConstant(number f,number xx)
    … …  
    240241RETURN:  number, the constant coefficient of the Laurent series of f in the variable x.
    241242THEORY:  Computes the constant coefficient of the Laurent series by iterative differentiation.
    242 
    243243KEYWORDS: Laurent series
    244244EXAMPLE:  example computeConstant; shows an example
    … …  
    250250  int j;
    251251  poly de;
    252   while (tst==0){
    253      ff=f*xx^k;
    254      for (j=1; j<=k; j++){
    255         ff=diff(ff,xx)/j;
    256      }
    257   de = subst(denominator(ff),xx,0);
    258   if (de!=0){
    259      poly nu = subst(numerator(ff),xx,0);
    260      return(number(nu/de));
    261   }
    262   k=k+1;
    263   }
    264 ERROR("error in computeConstant");}
     252  while (tst==0)
     253  {
     254    ff=f*xx^k;
     255    for (j=1; j<=k; j++)
     256    {
     257      ff=diff(ff,xx)/j;
     258    }
     259    de = subst(denominator(ff),xx,0);
     260    if (de!=0)
     261    {
     262      poly nu = subst(numerator(ff),xx,0);
     263      return(number(nu/de));
     264    }
     265    k++;
     266  }
     267  ERROR("error in computeConstant");
     268}
    265269example
    266270{ "EXAMPLE:"; echo=2;
    … …  
    270274  computeConstant(P,x2);
    271275}
    272 
    273 
    274 
    275276
    276277proc evaluateIntegral(number P, list xL)
    … …  
    289290{
    290291  number p = P;
    291   for(int j=1; j<=size(xL); j++){
     292  for(int j=1; j<=size(xL); j++)
     293  {
    292294     p=computeConstant(p,xL[j]);
    293295  }
    294 return(p);}
     296  return(p);
     297}
    295298example
    296299{ "EXAMPLE:"; echo=2;
    … …  
    300303  evaluateIntegral(p,list(x1,x3,x4,x2));
    301304}
    302 
    303305
    304306proc permute (list N)
    … …  
    314316"
    315317{
    316    int i,j,k;
    317    list L,L1;
    318    if (size(N)==1){
    319      return(list(N));
    320    } else {
    321      k=1;
    322      for (i=1; i<=size(N); i++){
    323        L=permute(delete(N,i));
    324        for (j=1; j<=size(L); j++){
    325           L1[k]=L[j]+list(N[i]);
    326           k=k+1;
    327        }
    328      }
    329    }
    330 return(L1);}
     318  int i,j,k;
     319  list L,L1;
     320  if (size(N)==1)
     321  {
     322    return(list(N));
     323  }
     324  else
     325  {
     326    k=1;
     327    for (i=1; i<=size(N); i++)
     328    {
     329      L=permute(delete(N,i));
     330      for (j=1; j<=size(L); j++)
     331      {
     332        L1[k]=L[j]+list(N[i]);
     333        k=k+1;
     334      }
     335    }
     336  }
     337  return(L1);
     338}
    331339example
    332340{ "EXAMPLE:"; echo=2;
    … …  
    349357"
    350358{
    351 ring R = 2,(x(1..n)),dp;
    352 ideal I = maxideal(a);
    353 list L;
    354 for (int j=1;j<=size(I);j++){
    355   L[j]=leadexp(I[j]);
    356 }
    357 return(L);}
     359  ring R = 2,(x(1..n)),dp;
     360  ideal I = maxideal(a);
     361  list L;
     362  for (int j=1;j<=size(I);j++)
     363  {
     364    L[j]=leadexp(I[j]);
     365  }
     366  return(L);
     367}
    358368example
    359369{ "EXAMPLE:"; echo=2;
    360370  partitions(3,7);
    361371}
    362 
    363 
    364372
    365373proc gromovWitten(def P,list #)
    … …  
    385393"
    386394{
    387   if (typeof(P)=="number") {
    388   list xl = ringlist(basering)[1][2];
    389   int j;
    390   for(j=1; j<=size(xl); j++){
    391      execute("number n= "+xl[j]);
    392      xl[j]=n;
    393      kill n;
    394   }
    395   list pxl = permute(xl);
    396   number p = 0;
    397   for(j=1; j<=size(pxl); j++){
    398      p=p+evaluateIntegral(P,pxl[j]);
    399   }
    400   return(p);
    401   }
    402   if (typeof(P)=="graph"){
    403      if (size(#)>1){
    404        return(gromovWitten(propagator(P,#)));
    405      } else {
     395  if (typeof(P)=="number")
     396  {
     397    list xl = ringlist(basering)[1][2];
     398    int j;
     399    for(j=1; j<=size(xl); j++)
     400    {
     401      execute("number n= "+xl[j]);
     402      xl[j]=n;
     403      kill n;
     404    }
     405    list pxl = permute(xl);
     406    number p = 0;
     407    for(j=1; j<=size(pxl); j++)
     408    {
     409      p=p+evaluateIntegral(P,pxl[j]);
     410    }
     411    return(p);
     412  }
     413  if (typeof(P)=="graph")
     414  {
     415    if (size(#)>1)
     416    {
     417      return(gromovWitten(propagator(P,#)));
     418    }
     419    else
     420    {
    406421      int d =#[1];
    407422      list pa = partitions(size(P.edges),d);
    408423      list re;
    409424      int ti;
    410       for (int j=1; j<=size(pa); j++) {
    411        ti=timer;
    412        re[j]=gromovWitten(propagator(P,pa[j]));
    413        ti=timer-ti;
    414        //print(string(j)+" / "+string(size(pa))+"    "+string(pa[j])+"     "+string(re[j])+"      "+string(sum(re))+"     "+string(ti));
     425      for (int j=1; j<=size(pa); j++)
     426      {
     427        ti=timer;
     428        re[j]=gromovWitten(propagator(P,pa[j]));
     429        ti=timer-ti;
     430        //print(string(j)+" / "+string(size(pa))+"    "+string(pa[j])+"     "+string(re[j])+"      "+string(sum(re))+"     "+string(ti));
    415431      }
    416      return(lsum(re));
    417      }
     432      return(lsum(re));
     433    }
    418434  }
    419435}
    … …  
    427443  gromovWitten(G,2);
    428444}
    429 
    430 
    431445
    432446proc computeGromovWitten(graph P,int d, int st, int en, list #)
    … …  
    443457"
    444458{
    445      number s =0;
    446      list pararg;
    447      list re;
    448      list pa = partitions(size(P.edges),d);
    449      int vb=0;
    450      if (size(#)>0){vb=#[1];}
    451      int ti;
    452      if (vb>0){print(size(pa));}
    453      for (int j=1; j<=size(pa); j++) {
    454       if ((j>=st)&(j<=en)){
    455        ti=timer;
    456        //pararg[j]=list(propagator(G,pa[j]));
    457        re[j]=gromovWitten(propagator(P,pa[j]));
    458        ti=timer-ti;
    459        if (vb>0){print(string(j)+" / "+string(size(pa))+"    "+string(pa[j])+"     "+string(re[j])+"      "+string(lsum(re))+"     "+string(ti));}
    460       } else {re[j]=s;}
    461      }
    462      //list re = parallelWaitAll("gromovWitten", pararg, list(list(list(2))));
    463      return(re);
     459  number s =0;
     460  list pararg;
     461  list re;
     462  list pa = partitions(size(P.edges),d);
     463  int vb=0;
     464  if (size(#)>0){vb=#[1];}
     465  int ti;
     466  if (vb>0){print(size(pa));}
     467  for (int j=1; j<=size(pa); j++)
     468  {
     469    if ((j>=st)&(j<=en))
     470    {
     471      ti=timer;
     472      //pararg[j]=list(propagator(G,pa[j]));
     473      re[j]=gromovWitten(propagator(P,pa[j]));
     474      ti=timer-ti;
     475      if (vb>0){print(string(j)+" / "+string(size(pa))+"    "+string(pa[j])+"     "+string(re[j])+"      "+string(lsum(re))+"     "+string(ti));}
     476    }
     477    else
     478    {re[j]=s;}
     479  }
     480  //list re = parallelWaitAll("gromovWitten", pararg, list(list(list(2))));
     481  return(re);
    464482}
    465483example
    … …  
    471489  computeGromovWitten(G,2,3,7,1);
    472490}
    473 
    474491
    475492proc lsum(list L)
    … …  
    485502{
    486503  execute(typeof(L[1])+" s");
    487   for(int j=1; j<=size(L); j++){
    488      s=s+L[j];
    489   }
    490 return(s);}
     504  for(int j=1; j<=size(L); j++)
     505  {
     506    s=s+L[j];
     507  }
     508  return(s);
     509}
    491510example
    492511{ "EXAMPLE:"; echo=2;
    … …  
    494513  lsum(L);
    495514}
    496 
    497 
    498515
    499516proc generatingFunction(graph G, int d)
    … …  
    512529  int j,jj;
    513530  list pa,L;
    514   for (j=1; j<=d; j++){
     531  for (j=1; j<=d; j++)
     532  {
    515533    pa = partitions(size(G.edges),j);
    516534    L = computeGromovWitten(G,j,1,size(pa));
    517     for (jj=1; jj<=size(pa); jj++) {
    518        s=s+L[jj]*monomial(pa[jj]);
    519     }
    520   }
    521 return(s);}
     535    for (jj=1; jj<=size(pa); jj++)
     536    {
     537      s=s+L[jj]*monomial(pa[jj]);
     538    }
     539  }
     540  return(s);
     541}
    522542example
    523543{ "EXAMPLE:"; echo=2;

  • Singular/dyn_modules/python/Makefile.am

    r112c79 r58603f  
    5151python_moduledir =  ${datadir}/singular/LIB
    5252python_module_DATA = cart.py  interpreter.py  perf.py  symm.py  util.py
     53EXTRA_DIST=cart.py  interpreter.py  perf.py  symm.py  util.py

  • Singular/dyn_modules/systhreads/Makefile.am

    r112c79 r58603f  
    1818endif
    1919
    20 SOURCES = shared.cc lintree.cc bytebuf.cc thread.cc bytebuf.h lintree.h channel.h syncvar.h threadconf.h
     20SOURCES = shared.cc lintree.cc bytebuf.cc thread.cc bytebuf.h lintree.h channel.h syncvar.h threadconf.h thread.h singthreads.h
    2121systhreads_la_SOURCES   = $(SOURCES)
    2222systhreads_la_CPPFLAGS  = ${MYINCLUDES} ${P_PROCS_CPPFLAGS_COMMON}

  • Tst/Short/absfact.res.gz.uu

    r112c79 r58603f  
    11begin 640 absfact.res.gz
    2 M'XL("`O%'%T"`V%B<V9A8W0N<F5S`.V<6V_;.!;'W_,IB&*`D2PI(0]OX@3V
     2M'XL("#A40%T"`V%B<V9A8W0N<F5S`.V<6V_;.!;'W_,IB&*`D2PI(0]OX@3V
    33MP^QB%P,L]J'IVZ`IY-2>&IM-@EC=2/WT>VC9)%W1M]:=+K9*45OFS1+/__QX
    44M1)&^>?/7W_Y)"&$3\H_??B6OZF5]>;^8OKJ^N%GGP(1@XKO%PZ).TNL+^TXF
    … …  
    4949M#*-1J03F2$&I4D`1:O[7!P3'NA0H5L*P2@L4:C^MC:390$H+)4JI..#(I:AF
    5050M`O`K5L,CU4H8#8:656:KEH*6"D^.E>:$R`NDJD1A7WFAN?T>C."-890:.P*#
    51 MPC:5+!4#$XGKR\"5G_W]GUV&&;6'84/4=E[I&]C5U;%%^LR(4W:?HR;,BKD5
    52 MQDIKP&XGM?TD5.WQ<9ZP<1[/I+"2MDHN(2(SXW<",!/^%-!ZOBO_M)[QZ@GH
    53 M-4EV%4Y=BYOM/W,LU`H$>-'"*&%%`T%W[MH.-#\E#AQH;;O21&:^@)XT56LG
    54 M*UM1K(;;8[76'^R!LLB&;:`P4.JL)@<:@Q%0\3U,[ED"%%D2P).EUS]=_!>^
    55 '1X0-\DP`````
     51MPC:5+!4#$XGKR\"5G_W]GUV&&;6'84/4=E[I&]C5U;%%^LR(H\(YK35L!74V
     52M6++T%;PJ4'-KV/92VV@JZEB=$OR)@J,<,U$J.PEF5C>9)8O-/QB_18"9\#>"
     53MUA-A^:?U5%A/6:])LJMPZEK<[`N:8Z%6(-F+%D8)*QH(^GG7/J'Y*0'B@'';
     54ME28R)0;TI#E<.XO9BF(U#A^KMWX4`)1%=G(#A0%?9S4YT!BE@(KO87+/$J#(
     551DH"J++W^Z>*_BSX7&PM-````
    5656`
    5757end

  • Tst/Short/absfact.stat

    r112c79 r58603f  
    1 1 >> tst_memory_0 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:190632
    2 1 >> tst_memory_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2486272
    3 1 >> tst_memory_2 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2486272
    4 1 >> tst_timer :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:18
    5 2 >> tst_memory_0 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:195048
    6 2 >> tst_memory_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    7 2 >> tst_memory_2 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    8 2 >> tst_timer :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:50
    9 3 >> tst_memory_0 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:195048
    10 3 >> tst_memory_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    11 3 >> tst_memory_2 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    12 3 >> tst_timer :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:1
    13 4 >> tst_memory_0 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:196464
    14 4 >> tst_memory_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    15 4 >> tst_memory_2 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    16 4 >> tst_timer :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:9
    17 5 >> tst_memory_0 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:195968
    18 5 >> tst_memory_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    19 5 >> tst_memory_2 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2740224
    20 5 >> tst_timer_1 :: 1562166539:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:83
     11 >> tst_memory_0 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:192200
     21 >> tst_memory_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2490368
     31 >> tst_memory_2 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2490368
     41 >> tst_timer :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:18
     52 >> tst_memory_0 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:196928
     62 >> tst_memory_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     72 >> tst_memory_2 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     82 >> tst_timer :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:45
     93 >> tst_memory_0 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:196928
     103 >> tst_memory_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     113 >> tst_memory_2 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     123 >> tst_timer :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:1
     134 >> tst_memory_0 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:210376
     144 >> tst_memory_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     154 >> tst_memory_2 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2621440
     164 >> tst_timer :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:10
     175 >> tst_memory_0 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:209856
     185 >> tst_memory_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2715648
     195 >> tst_memory_2 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2715648
     205 >> tst_timer_1 :: 1564496952:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:80

  • Tst/Short/bug_zp.res.gz.uu

    r112c79 r58603f  
    11begin 640 bug_zp.res.gz
    2 M'XL("-%LPED"`V)U9U]Z<"YR97,`E5!-3\)`$+WW5[P0#B6I"[O(AQ#6Q'@A
    3 M,5[P1K19RK9.`FTSNX3*KW>KI.K1T[QY;U[FS6Q>'M?/`*3&T_H!/>^\.-"N
    4 MMXPV5T5I!#*EDGP\6$9MA=;8G8KT4HO2GH7SQG?C8XT.WPH,AS@;+A-0CIKI
    5 M:$$.OJIP,%Q8Y!4C-YFO^&/1V28"3&4!7DW&T_E,CD8R:9)]_9-IJJ\NNEC$
    6 MC;J1(5G@[Y&]&PZ*97*>LC_+XJ-I6D*]J;O!/Z>W\G716M(`5C+:JN]6=HEF
    7 7O\Z>BZ^/M6\YN3A$ZT>?0,1X.V<!````
     2M'XL("$.:.5T"`V)U9U]Z<"YR97,`55#+;L(P$+S[*T:(0Y""P>99$$:J>D&J
     3M>H$;:B.3.F$E2"+;B)2OK].BM)QV=V9'L[/;W<OF#8!0>-T\H^.=YR<Z=)9L
     4M>V>D0@`3*LA'O25K*I3"X9(GMXH7YLJ=U[Y='RFT_9AC,,!5VR(&9:@LG0W(
     5MP9<E3MKF!EEID>G4E_9KT<HF'):*''8U&4WG,S$<BKB./ZN_FZ;JKJ*;053+
     6MO@B7!7R-]*AM8(PEYRE],(O.NFX`^2&?>FPOWA>-)`G-2K"]_!U%ZS'[%V3.
     73?W[0!+VX*)AUV3<,'@((.0$`````
    88`
    99end

  • Tst/Short/bug_zp.stat

    r112c79 r58603f  
    1 1 >> tst_memory_0 :: 1505914065:4103, 64 bit:4.1.0:x86_64-Linux:nepomuck:83592
    2 1 >> tst_memory_1 :: 1505914065:4103, 64 bit:4.1.0:x86_64-Linux:nepomuck:2215936
    3 1 >> tst_memory_2 :: 1505914065:4103, 64 bit:4.1.0:x86_64-Linux:nepomuck:2215936
    4 1 >> tst_timer_1 :: 1505914065:4103, 64 bit:4.1.0:x86_64-Linux:nepomuck:0
     11 >> tst_memory_0 :: 1564056131:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:85680
     21 >> tst_memory_1 :: 1564056131:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2150192
     31 >> tst_memory_2 :: 1564056131:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:2191296
     41 >> tst_timer_1 :: 1564056131:4120, 64 bit:4.1.2:x86_64-Linux:nepomuck:5

  • doc/NEWS.texi

    r112c79 r58603f  
    2121
    2222@heading News for version @value{VERSION}
     23
     24Changes in the kernel/build system:
     25@itemize
     26@item improved gcd and multiplication via FLINT
     27@end itemize
     28
     29@heading News for version 4-1-2
    2330New libraries:
    2431@itemize

  • factory/FLINTconvert.cc

    r112c79 r58603f  
    163163  if (f.isImm ())
    164164  {
    165     fmpz_set_si (fmpq_numref (result), f.intval());
    166     fmpz_set_si (fmpq_denref (result), 1);
     165    fmpq_set_si (result, f.intval(), 1);
    167166  }
    168167  else if(f.inQ())
    … …  
    182181    fmpz_set_mpz (fmpq_numref (result), gmp_val);
    183182    mpz_clear (gmp_val);
    184     fmpz_set_si (fmpq_denref (result), 1);
     183    fmpz_one(fmpq_denref(result));
    185184  }
    186185  else
    … …  
    574573  else
    575574  {
    576     int c=f.intval();
    577     if (c<0) c+=getCharacteristic();
     575    int c=f.intval(); // with Off(SW_SYMMETRIC_FF): 0<=c<p
    578576    nmod_mpoly_push_term_ui_ui(result,c,exp,ctx);
    579577  }
    … …  
    608606  ulong * exp = (ulong*)Alloc(N*sizeof(ulong));
    609607  memset(exp,0,N*sizeof(ulong));
     608  bool save_sym_ff= isOn (SW_SYMMETRIC_FF);
     609  if (save_sym_ff) Off (SW_SYMMETRIC_FF);
    610610  convFlint_RecPP( f, exp, res, ctx, N );
     611  if (save_sym_ff) On(SW_SYMMETRIC_FF);
    611612  Free(exp,N*sizeof(ulong));
    612613}
    … …  
    692693  convFactoryPFlintMP(F,f,ctx,N);
    693694  convFactoryPFlintMP(G,g,ctx,N);
    694   nmod_mpoly_init3(res,lF+lG,bits+1,ctx);
     695  nmod_mpoly_init(res,ctx);
    695696  nmod_mpoly_mul(res,f,g,ctx);
    696697  nmod_mpoly_clear(g,ctx);
    … …  
    713714  convFactoryPFlintMP(F,f,ctx,N);
    714715  convFactoryPFlintMP(G,g,ctx,N);
    715   fmpq_mpoly_init3(res,lF+lG,bits+1,ctx);
     716  fmpq_mpoly_init(res,ctx);
    716717  fmpq_mpoly_mul(res,f,g,ctx);
    717718  fmpq_mpoly_clear(g,ctx);
    … …  
    737738  convFactoryPFlintMP(F,f,ctx,N);
    738739  convFactoryPFlintMP(G,g,ctx,N);
    739   nmod_mpoly_init3(res,lf,bits,ctx);
     740  nmod_mpoly_init(res,ctx);
    740741  int ok=nmod_mpoly_gcd(res,f,g,ctx);
    741742  nmod_mpoly_clear(g,ctx);
    … …  
    750751  return RES;
    751752}
     753
     754static CanonicalForm b_content ( const CanonicalForm & f )
     755{
     756    if ( f.inCoeffDomain() )
     757        return f;
     758    else
     759    {
     760        CanonicalForm result = 0;
     761        CFIterator i;
     762        for ( i = f; i.hasTerms() && (!result.isOne()); i++ )
     763            result=bgcd( b_content(i.coeff()) , result );
     764        return result;
     765    }
     766}
     767
    752768
    753769CanonicalForm gcdFlintMP_QQ(const CanonicalForm& F, const CanonicalForm& G)
    … …  
    779795    }
    780796    RES=convFlintMPFactoryP(res,ctx,N);
     797    // gcd(2x,4x) should be 2x, so RES should also have the gcd(lc(F),lc(G))
     798    RES*=bgcd(b_content(F),b_content(G));
    781799  }
    782800  fmpq_mpoly_clear(res,ctx);

  • factory/Makefile.am

    r112c79 r58603f  
    101101# factory header files
    102102factory_headers = \
     103                globaldefs.h \
    103104                cf_assert.h \
    104105                canonicalform.h \
    … …  
    173174libfactory_includedir = ${includedir}/factory
    174175
    175 nodist_libfactory_include_HEADERS = factory.h factoryconf.h globaldefs.h
     176nodist_libfactory_include_HEADERS = factory.h factoryconf.h
    176177
    177178

  • factory/cf_char.cc

    r112c79 r58603f  
    3333        theDegree = 1;
    3434        CFFactory::settype( FiniteFieldDomain );
     35        ff_big = c > cf_getSmallPrime( cf_getNumSmallPrimes()-1 );
     36        if (c!=theCharacteristic)
     37        {
     38          if (c > 536870909) factoryError("characteristic is too large(max is 2^29)");
     39          ff_setprime( c );
     40        }
    3541        theCharacteristic = c;
    36         ff_big = c > cf_getSmallPrime( cf_getNumSmallPrimes()-1 );
    37         if (c > 536870909) factoryError("characteristic is too large(max is 2^29)");
    38         ff_setprime( c );
    3942    }
    4043}

  • factory/cf_gcd.cc

    r112c79 r58603f  
    104104    if ( isOn( SW_USE_FL_GCD_P)
    105105    && (CFFactory::gettype() != GaloisFieldDomain)
    106     && (getCharacteristic()>500)
     106    && (getCharacteristic()>10)
    107107    &&(!hasAlgVar(fc)) && (!hasAlgVar(gc)))
    108108    {
    … …  
    329329                CanonicalForm cdF = bCommonDen( f );
    330330                CanonicalForm cdG = bCommonDen( g );
     331                CanonicalForm F = f * cdF, G = g * cdG;
    331332                Off( SW_RATIONAL );
    332                 CanonicalForm l = lcm( cdF, cdG );
    333                 On( SW_RATIONAL );
    334                 CanonicalForm F = f * l, G = g * l;
    335                 Off( SW_RATIONAL );
    336                 l = gcd_poly( F, G );
     333                CanonicalForm l = gcd_poly( F, G );
    337334                On( SW_RATIONAL );
    338335                return abs( l );

  • factory/cf_switches.cc

    r112c79 r58603f  
    3232#ifdef HAVE_NTL
    3333  On(SW_USE_CHINREM_GCD);
    34   //Off(SW_USE_NTL_SORT);
     34  On(SW_USE_NTL_SORT);
     35#endif
     36#ifdef HAVE_FLINT
    3537  On(SW_USE_FL_GCD_P);
    36   //On(SW_USE_FL_GCD_0);
     38  On(SW_USE_FL_GCD_0);
    3739#endif
    3840  On(SW_USE_EZGCD);

  • kernel/mod2.h

    r112c79 r58603f  
    7777#endif
    7878
    79 #define SINGULAR_PATCHLEVEL 0
     79#define SINGULAR_PATCHLEVEL 2
    8080#define SINGULAR_VERSION ((SINGULAR_MAJOR_VERSION*1000 + SINGULAR_MINOR_VERSION*100 + SINGULAR_SUB_VERSION*10)+SINGULAR_PATCHLEVEL)
    8181

  • kernel/preimage.cc

    r112c79 r58603f  
    8888  }
    8989
    90   assume(n_SetMap(theImageRing->cf, dst_r->cf) == ndCopyMap);
    91 
    9290  if (theImageRing->cf != dst_r->cf)
    9391  {

  • libpolys/coeffs/Makefile.am

    r112c79 r58603f  
    4444test_LDADD   = libcoeffs.la $(libcoeffs_la_LIBADD)
    4545
    46 EXTRA_DIST = rintegers2.cc rinteger3.cc
     46EXTRA_DIST = rintegers2.cc rintegers3.cc

  • libpolys/coeffs/rintegers2.cc

    r112c79 r58603f  
    515515}
    516516
     517static void nrzMPZ(mpz_t res, number &a, const coeffs)
     518{
     519  mpz_init_set(res, (mpz_ptr) a);
     520}
     521
    517522static number nrzFarey(number r, number N, const coeffs R)
    518523{
    … …  
    621626  r->cfInit = nrzInit;
    622627  r->cfInitMPZ = nrzInitMPZ;
     628  r->cfMPZ = nrzMPZ;
    623629  r->cfSize  = nrzSize;
    624630  r->cfInt  = nrzInt;

  • libpolys/polys/clapconv.cc

    r112c79 r58603f  
    8787  p1=p;
    8888  l=l/2;
    89   while(l>0) { p=pNext(p); l--; }
     89  while(l>1) { p=pNext(p); l--; }
    9090  p2=pNext(p);
    9191  pNext(p)=NULL;
    … …  
    101101
    102102#define MIN_CONV_LEN 7
    103 CanonicalForm convSingPFactoryP( poly p, const ring r )
     103static CanonicalForm convSingPFactoryP_intern( poly p, int l, BOOLEAN & setChar,const ring r )
    104104{
    105105  CanonicalForm result = 0;
    106106  int e, n = rVar(r);
    107   BOOLEAN setChar=TRUE;
    108 
    109   int l;
    110   if ((l=pLength(p))>MIN_CONV_LEN)
     107  assume(l==pLength(p));
     108
     109  if (l>MIN_CONV_LEN)
    111110  {
    112111    poly p1,p2;
    113112    convPhalf(p,l,p1,p2);
    114     CanonicalForm P=convSingPFactoryP(p1,r);
    115     P+=convSingPFactoryP(p2,r);
     113    CanonicalForm P=convSingPFactoryP_intern(p1,l/2,setChar,r);
     114    P+=convSingPFactoryP_intern(p2,l-l/2,setChar,r);
    116115    convPunhalf(p1,p2);
    117116    return P;
    118117  }
     118  BOOLEAN setChar_loc=setChar;
     119  setChar=FALSE;
    119120  while ( p!=NULL )
    120121  {
    121     CanonicalForm term=r->cf->convSingNFactoryN(pGetCoeff( p ),setChar, r->cf);
     122    CanonicalForm term=r->cf->convSingNFactoryN(pGetCoeff( p ),setChar_loc, r->cf);
    122123    if (errorreported) break;
    123     setChar=FALSE;
     124    setChar_loc=FALSE;
    124125    for ( int i = 1; i <=n; i++ )
    125126    {
    … …  
    131132 }
    132133 return result;
     134}
     135
     136CanonicalForm convSingPFactoryP( poly p, const ring r )
     137{
     138  BOOLEAN setChar=TRUE;
     139  return convSingPFactoryP_intern(p,pLength(p),setChar,r);
    133140}
    134141

  • libpolys/polys/clapsing.cc

    r112c79 r58603f  
    5959  #ifdef HAVE_FLINT
    6060  #if __FLINT_RELEASE >= 20503
    61   if (rField_is_Zp(r) && (r->cf->ch>500))
     61  if (rField_is_Zp(r) && (r->cf->ch>10))
    6262  {
    6363    nmod_mpoly_ctx_t ctx;

  • libpolys/polys/flint_mpoly.cc

    r112c79 r58603f  
    6565#if 1
    6666// memory allocation is not thread safe; singular polynomials must be constructed in serial
     67
     68/*
     69    We agree the that result of a singular -> fmpq_mpoly conversion is
     70    readonly. This restricts the usage of the result in flint functions to
     71    const arguments. However, the real readonly conversion is currently only
     72    implemented in the threaded conversion below since it requires a scan of
     73    all coefficients anyways. The _fmpq_mpoly_clear_readonly_sing needs to
     74    be provided for a consistent interface in the polynomial operations.
     75*/
     76static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx)
     77{
     78    fmpq_mpoly_clear(a, ctx);
     79}
    6780
    6881void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r)
    … …  
    186199}
    187200
     201
     202/*
     203    In order that flint may sometimes borrow the large integer coeffs of
     204    polynomials over QQ (borrow means: simply point to the same GMP structs
     205    that singular has already allocated), we define the result of a
     206    singular -> fmpq_mpoly conversion to be readonly. This means we agree
     207    that
     208        - it can only be used as an const argument to a flint function
     209        - singular must not mutate the original coeffs while the readonly object is in use
     210*/
     211
     212static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx)
     213{
     214    if (fmpq_is_one(a->content))
     215    {
     216        if (a->zpoly->alloc > 0)
     217        {
     218            flint_free(a->zpoly->coeffs);
     219            flint_free(a->zpoly->exps);
     220        }
     221
     222        fmpq_clear(a->content);
     223    }
     224    else
     225    {
     226        fmpq_mpoly_clear(a, ctx);
     227    }
     228}
     229
    188230/* singular -> fmpq_mpoly conversion */
    189231
    … …  
    197239    std::vector<poly> markers;
    198240    ring r;
    199 
    200     convert_sing_to_fmpq_mpoly_base(slong num_threads_, fmpq_mpoly_struct * res_,
    201                              const fmpq_mpoly_ctx_struct * ctx_, const ring r_, poly p)
    202       : num_threads(num_threads_),
     241    fmpq_t content;
     242
     243    convert_sing_to_fmpq_mpoly_base(fmpq_mpoly_struct * res_,
     244                     const fmpq_mpoly_ctx_struct * ctx_, const ring r_, poly p)
     245      : num_threads(0),
    203246        res(res_),
    204247        ctx(ctx_),
    205248        r(r_)
    206249    {
     250        fmpq_t c;
     251        fmpq_init(c);
     252        fmpq_init(content);
     253        fmpq_zero(content);
     254
    207255        length = 0;
    208256        while (1)
    209257        {
    210258            if ((length % 4096) == 0)
     259            {
     260                my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     261                fmpq_gcd(content, content, c);
    211262                markers.push_back(p);
     263            }
    212264            if (p == NULL)
    213265                return;
    … …  
    215267            pIter(p);
    216268        }
     269
     270        fmpq_clear(c);
     271    }
     272
     273    ~convert_sing_to_fmpq_mpoly_base()
     274    {
     275        fmpq_clear(content);
    217276    }
    218277};
    … …  
    248307
    249308    flint_bitcnt_t required_bits = MPOLY_MIN_BITS;
    250     fmpq_zero(arg->content);
     309    fmpq_set(arg->content, base->content);
    251310
    252311    while (idx < arg->end_idx)
    253312    {
    254         my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    255         fmpq_gcd(arg->content, arg->content, c);
     313        number n = number(pGetCoeff(p));
     314
     315        if (fmpq_is_one(arg->content) && (SR_HDL(n)&SR_INT || n->s >= 3))
     316        {
     317            /* content is 1 and n is an integer, nothing to do */
     318        }
     319        else
     320        {
     321            my_convSingNFlintN_QQ(c, n);
     322            fmpq_gcd(arg->content, arg->content, c);
     323        }
    256324
    257325        #if SIZEOF_LONG==8
    … …  
    293361    }
    294362
     363    slong N = mpoly_words_per_exp(base->res->zpoly->bits, base->ctx->zctx->minfo);
     364    fmpz * res_coeffs = base->res->zpoly->coeffs;
     365    ulong * res_exps = base->res->zpoly->exps;
     366    flint_bitcnt_t res_bits = base->res->zpoly->bits;
     367
    295368    while (idx < arg->end_idx)
    296369    {
    297         my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    298         FLINT_ASSERT(!fmpq_is_zero(base->res->content));
    299         fmpq_div(t, c, base->res->content);
    300         FLINT_ASSERT(fmpz_is_one(fmpq_denref(t)));
    301 
    302         slong N = mpoly_words_per_exp(base->res->zpoly->bits, base->ctx->zctx->minfo);
    303         fmpz_swap(base->res->zpoly->coeffs + idx, fmpq_numref(t));
     370        if (fmpq_is_one(base->res->content))
     371        {
     372            // borrowing singular integers
     373            // the entry res_coeffs[idx] is junk, we should just overwrite it
     374
     375            number n = number(pGetCoeff(p));
     376
     377            if (SR_HDL(n)&SR_INT)
     378            {
     379                // n is a singular-small integer
     380                res_coeffs[idx] = SR_TO_INT(n);
     381            }
     382            else if (n->s<3)
     383            {
     384                // n is an element of QQ \ ZZ, should not happen
     385                assume(false);
     386            }
     387            else
     388            {
     389                // n is a singular-large integer, n may be flint-small
     390                res_coeffs[idx] = PTR_TO_COEFF(n->z);
     391                if (fmpz_fits_si(res_coeffs + idx))
     392                {
     393                    slong val = fmpz_get_si(res_coeffs + idx);
     394                    if (val >= COEFF_MIN && val <= COEFF_MAX)
     395                        res_coeffs[idx] = val;
     396                }
     397            }
     398        }
     399        else
     400        {
     401            my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     402            FLINT_ASSERT(!fmpq_is_zero(base->res->content));
     403            fmpq_div(t, c, base->res->content);
     404            FLINT_ASSERT(fmpz_is_one(fmpq_denref(t)));
     405            fmpz_swap(res_coeffs + idx, fmpq_numref(t));
     406        }
     407
    304408        #if SIZEOF_LONG==8
    305409        p_GetExpVL(p, (int64*)exp, base->r);
    306         mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, exp, base->res->zpoly->bits, base->ctx->zctx->minfo);
     410        mpoly_set_monomial_ui(res_exps + N*idx, exp, res_bits, base->ctx->zctx->minfo);
    307411        #else
    308412        p_GetExpV(p, (int*)exp, base->r);
    309         mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, &(exp[1]), base->res->zpoly->bits, base->ctx->minfo);
     413        mpoly_set_monomial_ui(res_exps + N*idx, &(exp[1]), res_bits, base->ctx->minfo);
    310414        #endif
    311415
    … …  
    324428    thread_pool_handle * handles;
    325429    slong num_handles;
    326     slong thread_limit = 1000;
    327 
    328     /* get workers */
     430    slong thread_limit = 1000; // TODO: should be paramter to this function
     431
     432    /* the constructor works out the length of p and sets some markers */
     433    convert_sing_to_fmpq_mpoly_base base(res, ctx, r, p);
     434
     435    /* sensibly limit thread count and get workers */
     436    thread_limit = FLINT_MIN(thread_limit, base.length/1024);
    329437    handles = NULL;
    330438    num_handles = 0;
    … …  
    340448    }
    341449
    342     convert_sing_to_fmpq_mpoly_base base(num_handles + 1, res, ctx, r, p);
    343 
    344450    /* fill in thread division points */
     451    base.num_threads = 1 + num_handles;
    345452    convert_sing_to_fmpq_mpoly_arg * args = new convert_sing_to_fmpq_mpoly_arg[base.num_threads];
    346453    slong cur_idx = 0;
    … …  
    348455    {
    349456        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads : base.length;
     457        FLINT_ASSERT(cur_idx <= base.length);
    350458        next_idx = FLINT_MAX(next_idx, cur_idx);
    351459        next_idx = FLINT_MIN(next_idx, base.length);
    … …  
    367475        required_bits = FLINT_MAX(required_bits, args[i].required_bits);
    368476
    369     /* initialize res with optimal bits */
    370     fmpq_mpoly_init3(res, base.length, mpoly_fix_bits(required_bits, ctx->zctx->minfo), ctx);
    371 
    372477    /* sign of content should match sign of first coeff */
    373     fmpq_zero(base.res->content);
     478    fmpq_t content;
     479    fmpq_init(content);
     480    fmpq_zero(content);
    374481    for (slong i = 0; i < base.num_threads; i++)
    375         fmpq_gcd(base.res->content, base.res->content, args[i].content);
     482        fmpq_gcd(content, content, args[i].content);
    376483    if (p != NULL)
    377484    {
    … …  
    380487        my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    381488        if (fmpq_sgn(c) < 0)
    382             fmpq_neg(base.res->content, base.res->content);
     489            fmpq_neg(content, content);
    383490        fmpq_clear(c);
    384491    }
     492
     493    /* initialize res with optimal bits */
     494    required_bits = mpoly_fix_bits(required_bits, ctx->zctx->minfo);
     495    if (fmpq_is_one(content))
     496    {
     497        /* initialize borrowed coeffs */
     498        slong N = mpoly_words_per_exp(required_bits, ctx->zctx->minfo);
     499        slong alloc = base.length;
     500        if (alloc != 0)
     501        {
     502            res->zpoly->coeffs = (fmpz *) flint_malloc(alloc*sizeof(fmpz));
     503            res->zpoly->exps   = (ulong *) flint_malloc(alloc*N*sizeof(ulong));
     504        }
     505        else
     506        {
     507            res->zpoly->coeffs = NULL;
     508            res->zpoly->exps = NULL;
     509        }
     510        res->zpoly->alloc = alloc;
     511        res->zpoly->length = 0;
     512        res->zpoly->bits = required_bits;
     513
     514        fmpq_init(res->content);
     515        fmpq_one(res->content);
     516    }
     517    else
     518    {
     519        /* initialize coeffs that will be created and destroyed */
     520        fmpq_mpoly_init3(res, base.length, required_bits, ctx);
     521        fmpq_swap(res->content, content);
     522    }
     523
     524    fmpq_clear(content);
    385525
    386526    /* fill in res->zpoly */
    … …  
    468608    thread_pool_handle * handles;
    469609    slong num_handles;
    470     slong thread_limit = 1000;
    471 
    472     /* get workers */
     610    slong thread_limit = 1000;// TODO: should be paramter to this function
     611
     612    /* sensibly limit threads and get workers */
     613    thread_limit = FLINT_MIN(thread_limit, f->zpoly->length/1024);
    473614    handles = NULL;
    474615    num_handles = 0;
    … …  
    485626
    486627    convert_fmpq_mpoly_to_sing_base base(num_handles + 1, f, ctx, r);
    487 
    488628    convert_fmpq_mpoly_to_sing_arg * args = new convert_fmpq_mpoly_to_sing_arg[base.num_threads];
    489629    slong cur_idx = 0;
    … …  
    492632        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->zpoly->length/base.num_threads
    493633                                                  : base.f->zpoly->length;
     634        FLINT_ASSERT(cur_idx <= base.f->zpoly->length);
    494635        next_idx = FLINT_MAX(next_idx, cur_idx);
    495636        next_idx = FLINT_MIN(next_idx, base.f->zpoly->length);
    … …  
    542683    ring r;
    543684
    544     convert_sing_to_nmod_mpoly_base(slong num_threads_, nmod_mpoly_struct * res_,
    545                             const nmod_mpoly_ctx_struct * ctx_, const ring r_, poly p)
    546       : num_threads(num_threads_),
     685    convert_sing_to_nmod_mpoly_base(nmod_mpoly_struct * res_,
     686                     const nmod_mpoly_ctx_struct * ctx_, const ring r_, poly p)
     687      : num_threads(0),
    547688        res(res_),
    548689        ctx(ctx_),
    … …  
    624765    }
    625766
     767    slong N = mpoly_words_per_exp(base->res->bits, base->ctx->minfo);
     768    ulong * res_coeffs = base->res->coeffs;
     769    ulong * res_exps = base->res->exps;
     770    flint_bitcnt_t res_bits = base->res->bits;
     771
    626772    while (idx < arg->end_idx)
    627773    {
    628         slong N = mpoly_words_per_exp(base->res->bits, base->ctx->minfo);
    629774        #if SIZEOF_LONG==8
    630775        p_GetExpVL(p, (int64*)exp, base->r);
    631         mpoly_set_monomial_ui(base->res->exps + N*idx, exp, base->res->bits, base->ctx->minfo);
     776        mpoly_set_monomial_ui(res_exps + N*idx, exp, res_bits, base->ctx->minfo);
    632777        #else
    633778        p_GetExpV(p, (int*)exp, base->r);
    634         mpoly_set_monomial_ui(base->res->exps + N*idx, &(exp[1]), base->res->bits, base->ctx->minfo);
     779        mpoly_set_monomial_ui(res_exps + N*idx, &(exp[1]), res_bits, base->ctx->minfo);
    635780        #endif
    636781
    637         base->res->coeffs[idx] = (ulong)(number(pGetCoeff(p)));
     782        res_coeffs[idx] = (ulong)(number(pGetCoeff(p)));
    638783
    639784        pIter(p);
    … …  
    649794    thread_pool_handle * handles;
    650795    slong num_handles;
    651     slong thread_limit = 1000;
    652 
    653     /* get workers */
     796    slong thread_limit = 1000; // TODO: should be paramter to this function
     797
     798    /* the constructor works out the length of p and sets some markers */
     799    convert_sing_to_nmod_mpoly_base base(res, ctx, r, p);
     800
     801    /* sensibly limit thread count and get workers */
     802    thread_limit = FLINT_MIN(thread_limit, base.length/1024);
    654803    handles = NULL;
    655804    num_handles = 0;
    … …  
    665814    }
    666815
    667     convert_sing_to_nmod_mpoly_base base(num_handles + 1, res, ctx, r, p);
    668 
    669816    /* fill in thread division points */
     817    base.num_threads = 1 + num_handles;
    670818    convert_sing_to_nmod_mpoly_arg * args = new convert_sing_to_nmod_mpoly_arg[base.num_threads];
    671819    slong cur_idx = 0;
    … …  
    674822        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads
    675823                                                  : base.length;
     824        FLINT_ASSERT(cur_idx <= base.length);
    676825        next_idx = FLINT_MAX(next_idx, cur_idx);
    677826        next_idx = FLINT_MIN(next_idx, base.length);
    … …  
    777926    thread_pool_handle * handles;
    778927    slong num_handles;
    779     slong thread_limit = 1000;
    780 
    781     /* get workers */
     928    slong thread_limit = 1000; // TODO: should be paramter to this function
     929
     930    /* sensibly limit threads and get workers */
     931    thread_limit = FLINT_MIN(thread_limit, f->length/1024);
    782932    handles = NULL;
    783933    num_handles = 0;
    … …  
    794944
    795945    convert_nmod_mpoly_to_sing_base base(num_handles + 1, f, ctx, r);
    796 
    797946    convert_nmod_mpoly_to_sing_arg * args = new convert_nmod_mpoly_to_sing_arg[base.num_threads];
    798947    slong cur_idx = 0;
    … …  
    801950        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->length/base.num_threads
    802951                                                  : base.f->length;
     952        FLINT_ASSERT(cur_idx <= base.f->length);
    803953        next_idx = FLINT_MAX(next_idx, cur_idx);
    804954        next_idx = FLINT_MIN(next_idx, base.f->length);
    … …  
    845995{
    846996  fmpq_mpoly_t pp,qq,res;
    847   convSingPFlintMP(pp,ctx,p,lp,r);
    848   convSingPFlintMP(qq,ctx,q,lq,r);
     997  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     998  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    849999  fmpq_mpoly_init(res,ctx);
    8501000  fmpq_mpoly_mul(res,pp,qq,ctx);
    8511001  poly pres=convFlintMPSingP(res,ctx,r);
    8521002  fmpq_mpoly_clear(res,ctx);
    853   fmpq_mpoly_clear(pp,ctx);
    854   fmpq_mpoly_clear(qq,ctx);
     1003  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     1004  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    8551005  fmpq_mpoly_ctx_clear(ctx);
    8561006  p_Test(pres,r);
    … …  
    8781028{
    8791029  fmpq_mpoly_t pp,qq,res;
    880   convSingPFlintMP(pp,ctx,p,lp,r);
    881   convSingPFlintMP(qq,ctx,q,lq,r);
     1030  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     1031  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    8821032  fmpq_mpoly_init(res,ctx);
    8831033  fmpq_mpoly_divides(res,pp,qq,ctx);
    8841034  poly pres = convFlintMPSingP(res,ctx,r);
    8851035  fmpq_mpoly_clear(res,ctx);
    886   fmpq_mpoly_clear(pp,ctx);
    887   fmpq_mpoly_clear(qq,ctx);
     1036  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     1037  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    8881038  fmpq_mpoly_ctx_clear(ctx);
    8891039  p_Test(pres,r);
    … …  
    9341084{
    9351085  fmpq_mpoly_t pp,qq,res;
    936   convSingPFlintMP(pp,ctx,p,lp,r);
    937   convSingPFlintMP(qq,ctx,q,lq,r);
     1086  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     1087  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    9381088  fmpq_mpoly_init(res,ctx);
    9391089  int ok=fmpq_mpoly_gcd(res,pp,qq,ctx);
    … …  
    9591109  }
    9601110  fmpq_mpoly_clear(res,ctx);
    961   fmpq_mpoly_clear(pp,ctx);
    962   fmpq_mpoly_clear(qq,ctx);
     1111  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     1112  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    9631113  fmpq_mpoly_ctx_clear(ctx);
    9641114  return pres;

  • libpolys/polys/operations/p_Mult_q.cc

    r112c79 r58603f  
    6767}
    6868
     69static void pqLengthApprox(poly p, poly q, int &lp, int &lq, const int min)
     70{
     71  int l = 0;
     72
     73  do
     74  {
     75    if (p == NULL)
     76    {
     77      lp=l;
     78      lq=l+(q!=NULL);
     79      return;
     80    }
     81    if (q == NULL) /* && p!=NULL */
     82    {
     83      lp=l+1;
     84      lq=l;
     85      return;
     86    }
     87    if (l>min) /* && p,q!=NULL */
     88    {
     89      lp=l; lq=l;
     90      return;
     91    }
     92    pIter(p);
     93    pIter(q);
     94    l++;
     95  }
     96  while (1);
     97}
     98
    6999
    70100static poly _p_Mult_q_Bucket(poly p, const int lp,
    … …  
    290320  poly pt;
    291321
    292   pqLength(p, q, lp, lq, MIN_LENGTH_FACTORY);
     322  // MIN_LENGTH_FACTORY must be >= MIN_LENGTH_FACTORY_QQ, MIN_FLINT_QQ, MIN_FLINT_Zp 20
     323  pqLengthApprox(p, q, lp, lq, MIN_LENGTH_FACTORY);
    293324
    294325  if (lp < lq)
    … …  
    309340    if (pure_polys && rField_is_Q(r) && !convSingRFlintR(ctx,r))
    310341    {
    311       lp=pLength(p);
    312       //printf("mul in flint\n");
    313       poly res=Flint_Mult_MP(p,lp,q,lq,ctx,r);
     342      // lq is a lower bound for the length of p and  q
     343      poly res=Flint_Mult_MP(p,lq,q,lq,ctx,r);
    314344      if (!copy)
    315345      {
    … …  
    325355    if (pure_polys && rField_is_Zp(r) && !convSingRFlintR(ctx,r))
    326356    {
    327       lp=pLength(p);
    328       //printf("mul in flint\n");
    329       poly res=Flint_Mult_MP(p,lp,q,lq,ctx,r);
     357      // lq is a lower bound for the length of p and  q
     358      poly res=Flint_Mult_MP(p,lq,q,lq,ctx,r);
    330359      if (!copy)
    331360      {
    … …  
    357386  {
    358387    lp=pLength(p);
    359     assume(lq == pLength(q));
     388    lq=pLength(q);
    360389    return _p_Mult_q_Bucket(p, lp, q, lq, copy, r);
    361390  }
Note: See TracChangeset for help on using the changeset viewer.