Changeset b75d13 in git


Ignore:
Timestamp:
May 12, 2005, 11:19:44 AM (19 years ago)
Author:
Michael Brickenstein <bricken@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
Children:
43bf9840cbfdc5c39cd127418a3d1d486fa793be
Parents:
15772766c4e2d15841019de318a47bafaffe477f
Message:
*bricken: externalized F4


git-svn-id: file:///usr/local/Singular/svn/trunk@8170 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • kernel/Makefile.in

    r1577276 rb75d13  
    5050CXXTEMPLFLAGS   = @CXXTEMPLFLAGS@
    5151CPPFLAGS        = -I${srcdir} @CPPFLAGS@
    52 DEFS            = -DNDEBUG -DOM_NDEBUG -D@SING_UNAME@ @DEFS@
     52DEFS            = -DNDEBUG -D@SING_UNAME@ @DEFS@
    5353LDFLAGS         = @LDFLAGS@
    5454LD_DYN_FLAGS    = @LD_DYN_FLAGS@
    … …  
    113113    pDebug.cc pInline2.cc pInline1.cc pInline0.cc \
    114114    pShallowCopyDelete.cc fast_mult.cc digitech.cc\
    115     tgb.cc tgbgauss.cc
     115    tgb.cc tgbgauss.cc F4.cc
    116116
    117117# normal C source files
    … …  
    452452CXXFLAGSP       = -pg -O3 ${PIPE}
    453453CXXTEMPLFLAGSP  = -fno-implicit-templates
    454 DEFSP           = -DNDEBUG -DOM_NDEBUG -DDO_PROFILE -D@SING_UNAME@ @DEFS@ -DDL_TAIL=$(DL_TAILP)
     454DEFSP           = -DNDEBUG -DDO_PROFILE -D@SING_UNAME@ @DEFS@ -DDL_TAIL=$(DL_TAILP)
    455455LDFLAGSP        = -static @LDFLAGS@
    456456
    … …  
    458458CXXFLAGSB       = -g -O3 ${PIPE}
    459459CXXTEMPLFLAGSB  = -fno-implicit-templates
    460 DEFSB           = -D@SING_UNAME@ -DOM_NDEBUG -DNDEBUG @DEFS@ -DDL_TAIL=$(DL_TAILB)
     460DEFSB           = -D@SING_UNAME@ -DNDEBUG @DEFS@ -DDL_TAIL=$(DL_TAILB)
    461461LDFLAGSB        = -static @LDFLAGS@
    462462

  • kernel/tgb.cc

    r1577276 rb75d13  
    55*  Computer Algebra System SINGULAR     *
    66****************************************/
    7 /* $Id: tgb.cc,v 1.21 2005-05-12 08:25:07 bricken Exp $ */
     7/* $Id: tgb.cc,v 1.22 2005-05-12 09:19:43 bricken Exp $ */
    88/*
    99* ABSTRACT: slimgb and F4 implementation
    … …  
    1313#include "tgb_internal.h"
    1414#include "tgbgauss.h"
     15#include "F4.h"
    1516static const int bundle_size=100;
    1617#if 1
    … …  
    210211  return(pLmCmp(((red_object*) ap)->p,((red_object*) bp)->p));
    211212}
    212 static int monom_poly_crit(const void* ap1, const void* ap2){
    213   monom_poly* p1;
    214   monom_poly* p2;
    215   p1=((monom_poly*) ap1);
    216   p2=((monom_poly*)ap2);
    217   if(((unsigned long) p1->f)>((unsigned long) p2->f)) return 1;
    218   if(((unsigned long) p1->f)<((unsigned long)p2->f)) return -1; 
    219 
    220   return pLmCmp(p1->m,p2->m);
    221  
    222 }
    223 static int pLmCmp_func(const void* ap1, const void* ap2){
    224     poly p1,p2;
    225   p1=*((poly*) ap1);
    226   p2=*((poly*)ap2);
    227 
    228   return pLmCmp(p1,p2);
    229 }
     213
     214
    230215static int pLmCmp_func_inverted(const void* ap1, const void* ap2){
    231216    poly p1,p2;
    … …  
    236221}
    237222
    238 static int pair_better_gen2(const void* ap,const void* bp){
    239   return(-pair_better_gen(ap,bp));
    240 }
    241 static int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj){
     223int tgb_pair_better_gen2(const void* ap,const void* bp){
     224  return(-tgb_pair_better_gen(ap,bp));
     225}
     226int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj){
    242227  int i;
    243228  long not_sev=~obj.sev;
    … …  
    249234  return -1;
    250235}
    251 static int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev){
     236int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev){
    252237  int i;
    253238  long not_sev=~sev;
    … …  
    284269    }
    285270}
    286 static int posInPolys (poly*  p, int pn, poly qe,calc_dat* c)
    287 {
    288   if(pn==0) return 0;
    289 
    290   int length=pn-1;
    291   int i;
    292   int an = 0;
    293   int en= length;
    294 
    295   if (pLmCmp(qe,p[en])==1)
    296     return length+1;
    297 
    298   while(1)
    299     {
    300       //if (an >= en-1)
    301       if(en-1<=an)
    302       {
    303         if (pLmCmp(p[an],qe)==1) return an;
    304         return en;
    305       }
    306       i=(an+en) / 2;
    307         if (pLmCmp(p[i],qe)==1)
    308           en=i;
    309       else an=i;
    310     }
    311 }
    312 static int posInMonomPolys (monom_poly*  p, int pn, monom_poly & qe,calc_dat* c)
    313 {
    314   if(pn==0) return 0;
    315 
    316   int length=pn-1;
    317   int i;
    318   int an = 0;
    319   int en= length;
    320 
    321   if (monom_poly_crit(&qe,&p[en])==1)
    322     return length+1;
    323 
    324   while(1)
    325     {
    326       //if (an >= en-1)
    327       if(en-1<=an)
    328       {
    329         if (monom_poly_crit(&p[an],&qe)==1) return an;
    330         return en;
    331       }
    332       i=(an+en) / 2;
    333         if (monom_poly_crit(&p[i],&qe)==1)
    334           en=i;
    335       else an=i;
    336     }
    337 }
     271
    338272static BOOLEAN  ascending(int* i,int top){
    339273  if(top<1) return TRUE;
    … …  
    342276}
    343277
    344 sorted_pair_node**  merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,calc_dat* c){
     278sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,calc_dat* c){
    345279  int i;
    346280  int* a= (int*) omalloc(qn*sizeof(int));
    … …  
    1021955}
    1022956
    1023 static sorted_pair_node** add_to_basis_ideal_quotient(poly h, int i_pos, int j_pos,calc_dat* c, int* ip)
     957sorted_pair_node** add_to_basis_ideal_quotient(poly h, int i_pos, int j_pos,calc_dat* c, int* ip)
    1024958{
    1025959
    … …  
    12241158  }
    12251159  if(!ip){
    1226     qsort(nodes_final,spc_final,sizeof(sorted_pair_node*),pair_better_gen2);
     1160    qsort(nodes_final,spc_final,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
    12271161 
    12281162   
    1229     c->apairs=merge(c->apairs,c->pair_top+1,nodes_final,spc_final,c);
     1163    c->apairs=spn_merge(c->apairs,c->pair_top+1,nodes_final,spc_final,c);
    12301164    c->pair_top+=spc_final;
    12311165    clean_top_of_pair_list(c);
    … …  
    17681702  }
    17691703
    1770   qsort(big_sbuf,sum,sizeof(sorted_pair_node*),pair_better_gen2);
    1771   c->apairs=merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
     1704  qsort(big_sbuf,sum,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
     1705  c->apairs=spn_merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
    17721706  c->pair_top+=sum;
    17731707  clean_top_of_pair_list(c);
    … …  
    19451879
    19461880//try to fill, return FALSE iff queue is empty
    1947 static void simplify(monom_poly& h, calc_dat* c){
    1948   mp_array_list* F=c->F;
    1949   poly_array_list* F_minus=c->F_minus;
    1950   while(F)
    1951   {
    1952     assume(F_minus!=NULL);
    1953     int i;
    1954     int posm=posInMonomPolys (F->mp, F->size, h,c);
    1955 #ifdef TGB_DEBUG
    1956 #endif
    1957      //for(i=0;i<F->size;i++)
    1958     for(i=si_min(posm,F->size-1);i>=0;i--)
    1959     {
    1960       //      if((i>=posm)&&(F->mp[i].f!=h.f)) break;
    1961       if((i<posm)&&(F->mp[i].f!=h.f)) break;
    1962       if ((h.f==F->mp[i].f) &&(p_LmDivisibleBy(F->mp[i].m,h.m,c->r)))
    1963       {
    1964        
    1965         //      Print("found");
    1966        
    1967           //according to the algorithm you should test (!(pIsConstant(F[i].m)))
    1968           //but I think this is only because of bad formulation
    1969         int j;
    1970        
    1971         poly lm=pLmInit(h.f);
    1972         pSetCoeff(lm,nInit(1));
    1973 
    1974         lm=pMult_mm(lm,F->mp[i].m);
    1975 
    1976         int pos;
    1977         j=-1;
    1978         pos=posInPolys (F_minus->p, F_minus->size, lm,c);
    1979       if((F_minus->size>pos)&&(pLmEqual(lm,F_minus->p[pos])))
    1980         j=pos;
    1981       if((pos>0) &&(pLmEqual(lm,F_minus->p[pos-1])))
    1982         j=pos-1;
    1983       assume(j!=-1);
    1984       //        if(j==-1) Print("\n jAltert \n");
    1985 //      for(j=0;j<F_minus->size;j++)
    1986 //      {
    1987 //        if (pLmEqual(F_minus->p[j],lm))
    1988 //          break;
    1989 //      }
    1990 //      assume(j<F_minus->size);
    1991         pDelete(&lm);
    1992         if(j>=0)
    1993         {
    1994           pExpVectorSub(h.m,F->mp[i].m);
    1995           h.f=F_minus->p[j];
    1996         }
    1997         break;
    1998       }
    1999     }
    2000     F=F->next;
    2001     F_minus=F_minus->next;
    2002   }
    2003   assume(F_minus==NULL);
    2004 }
    20051881
    20061882//transfers ownership of m to mat
    … …  
    20371913
    20381914
    2039 static tgb_matrix* build_matrix(poly* p,int p_index,poly* done, int done_index, calc_dat* c){
    2040   tgb_matrix* t=new tgb_matrix(p_index,done_index);
    2041   int i, pos;
    2042   //  Print("\n 0:%s\n",pString(done[done_index-1]));
    2043   //Print("\n 1:%s\n",pString(done[done_index-2]));
    2044   assume((!(pLmEqual(done[done_index-1],done[done_index-2]))));
    2045 #ifdef TGB_DEGUG
    2046   for(i=0;i<done_index;i++)
    2047   {
    2048     int j;
    2049     for(j=0;j<i;j++)
    2050     {
    2051       assume((!(pLmEqual(done[i],done[j]))));
    2052     }
    2053   }
    2054 #endif
    2055   for(i=0;i<p_index;i++)
    2056   {
    2057     // Print("%i ter Eintrag:%s\n",i,pString(p[i]));
    2058     poly p_i=p[i];
    2059     while(p_i)
    2060     {
    2061 
    2062       int v=-1;
    2063       pos=posInPolys (done, done_index, p_i,c);
    2064       if((done_index>pos)&&(pLmEqual(p_i,done[pos])))
    2065         v=pos;
    2066       if((pos>0) &&(pLmEqual(p_i,done[pos-1])))
    2067         v=pos-1;
    2068       assume(v!=-1);
    2069       //v is ascending ordered, we need descending order
    2070       v=done_index-1-v;
    2071       number nt=t->get(i,v);
    2072       nDelete(&nt);
    2073       t->set(i,v,nCopy(p_i->coef));
    2074       p_i=p_i->next;
    2075     }
    2076   }
    2077   return t;
    2078 }
    2079 static tgb_sparse_matrix* build_sparse_matrix(poly* p,int p_index,poly* done, int done_index, calc_dat* c){
    2080   tgb_sparse_matrix* t=new tgb_sparse_matrix(p_index,done_index,c->r);
    2081   int i, pos;
    2082   //  Print("\n 0:%s\n",pString(done[done_index-1]));
    2083   //Print("\n 1:%s\n",pString(done[done_index-2]));
    2084   //  assume((!(pLmEqual(done[done_index-1],done[done_index-2]))));
    2085 #ifdef TGB_DEGUG
    2086   for(i=0;i<done_index;i++)
    2087   {
    2088     int j;
    2089     for(j=0;j<i;j++)
    2090     {
    2091       assume((!(pLmEqual(done[i],done[j]))));
    2092     }
    2093   }
    2094 #endif
    2095   for(i=0;i<p_index;i++)
    2096   {
    2097     // Print("%i ter Eintrag:%s\n",i,pString(p[i]));
    2098     mac_poly m=NULL;
    2099     mac_poly* set_this=&m;
    2100     poly p_i=p[i];
    2101     while(p_i)
    2102     {
    2103 
    2104       int v=-1;
    2105       pos=posInPolys (done, done_index, p_i,c);
    2106       if((done_index>pos)&&(pLmEqual(p_i,done[pos])))
    2107         v=pos;
    2108       if((pos>0) &&(pLmEqual(p_i,done[pos-1])))
    2109         v=pos-1;
    2110       assume(v!=-1);
    2111       //v is ascending ordered, we need descending order
    2112       v=done_index-1-v;
    2113       (*set_this)=new mac_poly_r();
    2114       (*set_this)->exp=v;
    2115        
    2116       (*set_this)->coef=nCopy(p_i->coef);
    2117        set_this=&(*set_this)->next;
    2118        p_i=p_i->next;
    2119 
    2120     }
    2121     init_with_mac_poly(t,i,m);
    2122   }
    2123   return t;
    2124 }
    2125 
    2126 
    2127 static int retranslate(poly* m,tgb_sparse_matrix* mat,poly* done, calc_dat* c){
    2128   int i;
    2129   int m_index=0;
    2130   for(i=0;i<mat->get_rows();i++)
    2131   {
    2132     if(mat->zero_row(i))
    2133     {
    2134       mat->free_row(i);
    2135       continue;
    2136     }
    2137     m[m_index++]= free_row_to_poly(mat, i, done, mat->get_columns());
    2138 
    2139   }
    2140 
    2141   delete mat;
    2142   return m_index;
    2143 
    2144 }
    2145 
    2146 //returns m_index and destroys mat
    2147 static int retranslate(poly* m,tgb_matrix* mat,poly* done, calc_dat* c){
    2148   int i;
    2149   int m_index=0;
    2150   for(i=0;i<mat->get_rows();i++)
    2151   {
    2152     if(mat->zero_row(i))
    2153     {
    2154       mat->free_row(i);
    2155       continue;
    2156     }
    2157    
    2158     m[m_index]=pInit();
    2159     int v=mat->min_col_not_zero_in_row(i);
    2160     //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v
    2161     int pos=mat->get_columns()-1-v;
    2162     int* ev=(int*) omalloc((c->r->N+1)*sizeof(int));
    2163     pGetExpV(done[pos],ev);
    2164     pSetExpV(m[m_index],ev);
    2165     omfree(ev);
    2166    
    2167     poly p=m[m_index];
    2168     pSetCoeff(p,mat->get(i,v));
    2169     while((v=mat->next_col_not_zero(i,v))!=mat->get_columns())
    2170     {
    2171       poly pn=pInit();
    2172      
    2173       //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v
    2174       pos=mat->get_columns()-1-v;
    2175       ev=(int*) omalloc((c->r->N+1)*sizeof(int));
    2176       pGetExpV(done[pos],ev);
    2177       pSetExpV(pn,ev);
    2178       omfree(ev);
    2179       pSetCoeff(pn,mat->get(i,v));
    2180       p->next=pn;
    2181       p=pn;
    2182     }
    2183     p->next=NULL;
    2184     mat->free_row(i, FALSE);
    2185     m_index++;
    2186   }
    2187 
    2188   delete mat;
    2189   return m_index;
    2190 
    2191 }
    2192 static void go_on_F4 (calc_dat* c){
    2193   //set limit of 1000 for multireductions, at the moment for
    2194   //programming reasons
    2195   int max_par;
    2196   if (c->is_homog)
    2197     max_par=PAR_N_F4;
    2198   else
    2199     max_par=200;
    2200   int done_size=max_par*2;
    2201   poly* done=(poly*) omalloc(done_size*sizeof(poly));
    2202   int done_index=0; //done_index must always be smaller than done_size
    2203   int chosen_size=max_par*2;
    2204   monom_poly* chosen=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly));
    2205   int chosen_index=0;
    2206   //  monom_poly* vgl=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly));
    2207   int i=0;
    2208   c->average_length=0;
    2209   for(i=0;i<c->n;i++){
    2210     c->average_length+=c->lengths[i];
    2211   }
    2212   c->average_length=c->average_length/c->n;
    2213   i=0;
    2214   poly* p;
    2215   int nfs=0;
    2216   int curr_deg=-1;
    2217 
    2218   //choose pairs and preprocess symbolically
    2219   while(chosen_index<max_par)
    2220   {
    2221    
    2222     //    sorted_pair_node* s=c->apairs[c->pair_top];
    2223     sorted_pair_node* s;
    2224     if(c->pair_top>=0)
    2225       s=top_pair(c);//here is actually chain criterium done
    2226     else
    2227       s=NULL;
    2228     if (!s) break;
    2229     nfs++;
    2230     if(curr_deg>=0)
    2231     {
    2232       if (s->deg >curr_deg) break;
    2233     }
    2234 
    2235     else curr_deg=s->deg;
    2236     quick_pop_pair(c);
    2237     if(s->i>=0)
    2238     {
    2239       //replace_pair(s->i,s->j,c);
    2240       if(s->i==s->j)
    2241       {
    2242         free_sorted_pair_node(s,c->r);
    2243         continue;
    2244       }
    2245     }
    2246     assume(s->i>=0);
    2247     monom_poly h;
    2248     BOOLEAN only_free=FALSE;
    2249     if(s->i>=0)
    2250     {
    2251      
    2252       poly lcm=pOne();
    2253      
    2254       pLcm(c->S->m[s->i], c->S->m[s->j], lcm);
    2255       pSetm(lcm);
    2256       assume(lcm!=NULL);
    2257       poly factor1=pCopy(lcm);
    2258       pExpVectorSub(factor1,c->S->m[s->i]);
    2259       poly factor2=pCopy(lcm);
    2260       pExpVectorSub(factor2,c->S->m[s->j]);
    2261 
    2262       if(done_index>=done_size)
    2263       {
    2264         done_size+=max_par;
    2265         done=(poly*) omrealloc(done,done_size*sizeof(poly));
    2266       }
    2267       done[done_index++]=lcm;
    2268       if(chosen_index+1>=chosen_size)
    2269       {
    2270         //max_par must be greater equal 2
    2271         chosen_size+=si_max(max_par,2);
    2272         chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
    2273       }
    2274       h.f=c->S->m[s->i];
    2275       h.m=factor1;
    2276       chosen[chosen_index++]=h;
    2277       h.f=c->S->m[s->j];
    2278       h.m=factor2;
    2279       chosen[chosen_index++]=h;
    2280     }
    2281     else
    2282     {
    2283       if(chosen_index>=chosen_size)
    2284       {
    2285         chosen_size+=max_par;
    2286         chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
    2287       }
    2288       h.f=s->lcm_of_lm;
    2289       h.m=pOne();
    2290       //pSetm(h.m); done by pOne
    2291       chosen[chosen_index++]=h;
    2292       //must carefull remember to destroy such a h;
    2293       poly_list_node* next=c->to_destroy;
    2294      
    2295       c->to_destroy=(poly_list_node*) omalloc(sizeof(poly_list_node));
    2296       c->to_destroy->p=h.f;
    2297       c->to_destroy->next=next;
    2298       only_free=TRUE;
    2299     }
    2300 
    2301     if(s->i>=0)
    2302       now_t_rep(s->j,s->i,c);
    2303 
    2304     if(!(only_free))
    2305       free_sorted_pair_node(s,c->r);
    2306     else
    2307       omfree(s);
    2308    
    2309  
    2310 
    2311   }
    2312   c->normal_forms+=nfs;
    2313   if (TEST_OPT_PROT)
    2314       Print("M[%i, ",nfs);
    2315   //next Step, simplify all pairs
    2316   for(i=0;i<chosen_index;i++)
    2317   {
    2318     //    PrintS("simplifying ");
    2319     //Print("from %s",pString(chosen[i].m));
    2320     //Print(", %s", pString(chosen[i].f));
    2321      simplify(chosen[i],c);
    2322     //Print(" to %s",pString(chosen[i].m));
    2323     //Print(", %s \n", pString(chosen[i].f));
    2324   }
    2325  
    2326   //next Step remove duplicate entries
    2327   qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit);
    2328   int pos=0;
    2329   for(i=1;i<chosen_index;i++)
    2330   {
    2331     if(((chosen[i].f!=chosen[pos].f))||(!(pLmEqual(chosen[i].m,chosen[pos].m))))
    2332       chosen[++pos]=chosen[i];
    2333     else pDelete(&(chosen[i].m));
    2334   }
    2335   if(chosen_index>0)
    2336     chosen_index=pos+1;
    2337   //next step process polys
    2338   int p_size=2*chosen_index;
    2339   int p_index=0;
    2340   p=(poly*) omalloc(p_size*sizeof(poly));
    2341   for(p_index=0;p_index<chosen_index;p_index++)
    2342   {
    2343     p[p_index]=ppMult_mm(chosen[p_index].f,chosen[p_index].m);
    2344   }
    2345 
    2346   qsort(done, done_index,sizeof(poly),pLmCmp_func);
    2347   pos=0;
    2348   //Print("Done_index:%i",done_index);
    2349   if(done_index>0)
    2350   {
    2351     pTest(done[0]);
    2352   }
    2353   for(i=1;i<done_index;i++)
    2354   {
    2355     pTest(done[i]);
    2356     if((!(pLmEqual(done[i],done[pos]))))
    2357       done[++pos]=done[i];
    2358     else pDelete(&(done[i]));
    2359   }
    2360   if(done_index>0)
    2361     done_index=pos+1;
    2362 #ifdef TGB_DEBUG
    2363   for(i=0;i<done_index;i++)
    2364   {
    2365     int j;
    2366     for(j=0;j<i;j++)
    2367     {
    2368       assume((!pLmEqual(done[j],done[i])));
    2369     }
    2370   }
    2371 #endif
    2372 #ifdef TGB_DEBUG
    2373   for(i=0;i<done_index;i++)
    2374   {
    2375     pTest(done[i]);
    2376   }
    2377 #endif
    2378   poly* m;
    2379   int m_index=0;
    2380   int m_size=0;
    2381   for(i=0;i<p_index;i++)
    2382   {
    2383     m_size+=pLength(p[i]);
    2384   }
    2385   m=(poly*) omalloc(m_size*sizeof(poly));
    2386   //q=(poly*) omalloc(m_size*sizeof(poly));
    2387 
    2388  
    2389 
    2390   for(i=0;i<p_index;i++)
    2391   {
    2392 
    2393     poly p_i=p[i];
    2394 
    2395     pTest(p[i]);
    2396 
    2397     while(p_i)
    2398     {
    2399 
    2400       m[m_index]=pLmInit(p_i);
    2401 
    2402       pSetCoeff(m[m_index],nInit(1));
    2403 
    2404        p_i=p_i->next;
    2405 
    2406       m_index++;
    2407 
    2408     }
    2409   }
    2410 
    2411   int q_size=m_index;
    2412   poly* q=(poly*) omalloc(q_size*sizeof(poly));
    2413   int q_index=0;
    2414   //next Step additional reductors
    2415 
    2416   while(m_index>0)
    2417   {
    2418 #ifdef TGB_DEBUG
    2419      
    2420       for(i=0;i<done_index;i++)
    2421       {
    2422         pTest(done[i]);
    2423       }
    2424 #endif
    2425 
    2426     qsort(m, m_index,sizeof(poly),pLmCmp_func);
    2427 
    2428    
    2429     pos=0;
    2430     #ifdef TGB_DEBUG
    2431      
    2432       for(i=0;i<done_index;i++)
    2433       {
    2434        
    2435         pTest(done[i]);
    2436       }
    2437 #endif
    2438     for(i=1;i<m_index;i++)
    2439     {
    2440       pTest(m[i]);
    2441       pTest(m[pos]);
    2442       if((!(pLmEqual(m[i],m[pos]))))
    2443         m[++pos]=m[i];
    2444       else pDelete(&(m[i]));
    2445     }
    2446     if(m_index>1)
    2447       m_index=pos+1;
    2448     if(done_size<done_index+m_index)
    2449     {
    2450       done_size=done_index+2*m_index;
    2451       done=(poly*) omrealloc(done,done_size*sizeof(poly));
    2452     }
    2453     if(chosen_size<chosen_index+m_index)
    2454     {
    2455       chosen_size=chosen_index+2*m_index;
    2456       chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
    2457     }
    2458     q_index=0;
    2459     if(q_size<m_index)
    2460     {
    2461       q_size=2*m_index;
    2462       omfree(q);
    2463       q=(poly*) omalloc(q_size*sizeof(poly));
    2464     }
    2465 
    2466     for(i=0;i<m_index;i++)
    2467     {
    2468 
    2469 #ifdef TGB_DEBUG
    2470       {
    2471         int my_i;
    2472         for(my_i=0;my_i<done_index;my_i++)
    2473         {
    2474           int my_j;
    2475           for(my_j=0;my_j<my_i;my_j++)
    2476           {
    2477             assume((!pLmEqual(done[my_j],done[my_i])));
    2478           }
    2479         }
    2480       }
    2481 #endif
    2482       BOOLEAN in_done=FALSE;
    2483       pTest(m[i]);
    2484 
    2485       pos=posInPolys (done, done_index, m[i],c);
    2486 #ifdef TGB_DEBUG
    2487       {
    2488         int my_i;
    2489         for (my_i=0;my_i<done_index;my_i++)
    2490           pTest(done[my_i]);
    2491       }
    2492 #endif     
    2493       if(((done_index>pos)&&(pLmEqual(m[i],done[pos]))) ||(pos>0) &&(pLmEqual(m[i],done[pos-1])))
    2494         in_done=TRUE;
    2495       if (!(in_done))
    2496       {
    2497        
    2498         int S_pos=kFindDivisibleByInS_easy(c->strat,m[i], pGetShortExpVector(m[i]));
    2499         if(S_pos>=0)
    2500         {
    2501           monom_poly h;
    2502           h.f=c->strat->S[S_pos];
    2503           h.m=pOne();
    2504           int* ev=(int*) omalloc((c->r->N+1)*sizeof(int));
    2505           pGetExpV(m[i],ev);
    2506           pSetExpV(h.m,ev);
    2507           omfree(ev);
    2508           pExpVectorSub(h.m,c->strat->S[S_pos]);
    2509           simplify(h,c);
    2510           q[q_index]=ppMult_mm(h.f,h.m);
    2511           chosen[chosen_index++]=h;
    2512           q_index++;
    2513         }
    2514         pTest(m[i]);
    2515 #ifdef TGB_DEBUG
    2516         {
    2517           int my_i;
    2518           for (my_i=0;my_i<done_index;my_i++)
    2519             pTest(done[my_i]);
    2520         }
    2521 #endif     
    2522         memmove(&(done[pos+1]),&(done[pos]), (done_index-pos)*sizeof(poly));
    2523         done[pos]=m[i];
    2524         done_index++;
    2525 #ifdef TGB_DEBUG
    2526         {
    2527           int my_i;
    2528           for (my_i=0;my_i<done_index;my_i++)
    2529           {
    2530             //      Print("Position %i pos %i size %i\n",my_i,pos,done_index);
    2531             pTest(done[my_i]);
    2532           }
    2533         }
    2534 #endif     
    2535       }
    2536       else
    2537         pDelete(&m[i]);
    2538 #ifdef TGB_DEBUG
    2539       {
    2540         int my_i;
    2541         for(my_i=0;my_i<done_index;my_i++)
    2542         {
    2543           pTest(done[my_i]);
    2544           int my_j;
    2545           for(my_j=0;my_j<my_i;my_j++)
    2546           {
    2547             assume((!pLmEqual(done[my_j],done[my_i])));
    2548           }
    2549         }
    2550       }
    2551 #endif
    2552     }
    2553     if(p_size<p_index+q_index)
    2554     {
    2555       p_size=p_index+2*q_index;
    2556       p=(poly*) omrealloc(p,p_size*sizeof(poly));
    2557     }
    2558     for (i=0;i<q_index;i++)
    2559       p[p_index++]=q[i];
    2560     m_index=0;
    2561     int sum=0;
    2562     for(i=0;i<p_index;i++)
    2563     {
    2564       sum+=pLength(p[i])-1;
    2565     }
    2566     if (m_size<sum)
    2567     {
    2568       omfree(m);
    2569       m=(poly*) omalloc(sum*sizeof(poly));
    2570     }
    2571     m_size=sum;
    2572     for(i=0;i<q_index;i++)
    2573     {
    2574       poly p_i=q[i]->next;
    2575       while(p_i)
    2576       {
    2577         m[m_index]=pLmInit(p_i);
    2578         pSetCoeff(m[m_index],nInit(1));
    2579         p_i=p_i->next;
    2580         m_index++;
    2581       }
    2582     }
    2583     //qsort(done, done_index,sizeof(poly),pLmCmp_func);
    2584 #ifdef TGB_DEBUG
    2585     for(i=0;i<done_index;i++)
    2586     {
    2587       pTest(done[i]);
    2588     }
    2589 #endif
    2590   }
    2591 #ifdef TGB_DEBUG
    2592   for(i=0;i<done_index;i++)
    2593   {
    2594     int j;
    2595     for(j=0;j<i;j++)
    2596     {
    2597       assume((!pLmEqual(done[j],done[i])));
    2598     }
    2599   }
    2600 #endif
    2601   omfree(m);
    2602   omfree(q);
    2603   if (TEST_OPT_PROT)
    2604     Print("Mat[%i x %i], ",chosen_index, done_index);
    2605   //next Step build matrix
    2606   #ifdef TGB_DEBUG
    2607   for(i=0;i<done_index;i++)
    2608   {
    2609     int j;
    2610     for(j=0;j<i;j++)
    2611     {
    2612       assume((!pLmEqual(done[j],done[i])));
    2613     }
    2614   }
    2615 #endif
    2616   assume(p_index==chosen_index);
    2617  
    2618  //  tgb_matrix* mat=build_matrix(p,p_index,done, done_index,c);
    2619  
    2620 //   simple_gauss2(mat);
    2621   tgb_sparse_matrix* mat=build_sparse_matrix(p,p_index,done, done_index,c);
    2622   simple_gauss(mat,c);
    2623   m_size=mat->get_rows();
    2624   m=(poly*) omalloc(m_size*sizeof(poly));
    2625   m_index=retranslate(m,mat,done,c);
    2626  
    2627   mat=NULL;
    2628   for(i=0;i<done_index;i++)
    2629     pDelete(&done[i]);
    2630   omfree(done);
    2631   done=NULL;
    2632   //next Step addElements to basis
    2633   int F_plus_size=m_index;
    2634   poly* F_plus=(poly*)omalloc(F_plus_size*sizeof(poly));
    2635   int F_plus_index=0;
    2636   int F_minus_size=m_index;
    2637   poly* F_minus=(poly*) omalloc(F_minus_size*sizeof(poly));
    2638   int F_minus_index=0;
    2639 
    2640   //better algorithm replace p by its monoms, qsort,delete duplicates and binary search for testing if monom is contained in array
    2641   qsort(p, p_index,sizeof(poly),pLmCmp_func);
    2642   for(i=0;i<p_index;i++)
    2643     pDelete(&p[i]->next);
    2644   for(i=m_index-1;i>=0;--i)
    2645   {
    2646     int j;
    2647     int pos=posInPolys (p,p_index, m[i],c);
    2648     BOOLEAN minus=FALSE;
    2649     if(((p_index>pos)&&(pLmEqual(m[i],p[pos]))) ||(pos>0) &&(pLmEqual(m[i],p[pos-1])))
    2650       minus=TRUE;
    2651    
    2652     if(minus)
    2653     {
    2654       F_minus[F_minus_index++]=m[i];
    2655       m[i]=NULL;
    2656     }
    2657     else
    2658     {
    2659       F_plus[F_plus_index++]=m[i];
    2660       m[i]=NULL;
    2661     }
    2662   }
    2663 
    2664 //   for(i=m_index-1;i>=0;--i)
    2665 //   {
    2666 //     int j;
    2667 //     BOOLEAN minus=FALSE;
    2668 //     for(j=0;j<p_index;j++)
    2669 //       if (pLmEqual(p[j],m[i]))
    2670 //       {
    2671 //      minus=TRUE;
    2672 //      break;
    2673 //       }
    2674 //     if(minus)
    2675 //     {
    2676 //       F_minus[F_minus_index++]=m[i];
    2677 //       m[i]=NULL;
    2678 //     }
    2679 //     else
    2680 //     {
    2681 //       F_plus[F_plus_index++]=m[i];
    2682 //       m[i]=NULL;
    2683 //     }
    2684 //   }
    2685   //in this order F_minus will be automatically ascending sorted
    2686   //to make this sure for foreign gauss
    2687   //uncomment the following line
    2688   //  qsort(F_minus, F_minus_index,sizeof(poly),pLmCmp_func);
    2689   assume((F_plus_index+F_minus_index)==m_index);
    2690   if (TEST_OPT_PROT)
    2691     Print("%i]", F_plus_index);
    2692   for(i=0;i<p_index;i++)
    2693     pDelete(&p[i]);
    2694   omfree(p);
    2695   p=NULL;
    2696   omfree(m);
    2697   m=NULL;
    2698 
    2699   //the F_minus list must be cleared separately at the end
    2700   mp_array_list** F_i;
    2701   poly_array_list** F_m_i;
    2702   F_i=&(c->F);
    2703   F_m_i=&(c->F_minus);
    2704 
    2705   while((*F_i)!=NULL)
    2706   {
    2707     assume((*F_m_i)!=NULL);
    2708     F_i=(&((*F_i)->next));
    2709     F_m_i=(&((*F_m_i)->next));
    2710   }
    2711   assume((*F_m_i)==NULL);
    2712   //should resize the array to save memory
    2713   //F and F_minus
    2714   qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit);//important for simplify
    2715   (*F_m_i)=(poly_array_list*) omalloc(sizeof(poly_array_list));
    2716   (*F_m_i)->size=F_minus_index;
    2717   (*F_m_i)->p=F_minus;
    2718   (*F_m_i)->next=NULL;
    2719   (*F_i)=(mp_array_list*) omalloc(sizeof(poly_array_list));
    2720   (*F_i)->size=chosen_index;
    2721   (*F_i)->mp=chosen;
    2722   (*F_i)->next=NULL;
    2723  
    2724   if(F_plus_index>0)
    2725   {
    2726     int j;
    2727     int* ibuf=(int*) omalloc(F_plus_index*sizeof(int));
    2728     sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(F_plus_index*sizeof(sorted_pair_node**));
    2729  
    2730     for(j=0;j<F_plus_index;j++)
    2731     {
    2732       int len;
    2733       poly p=F_plus[j];
    2734    
    2735     // delete buf[j];
    2736     //remember to free res here
    2737     //    p=redTailShort(p, c->strat);
    2738       sbuf[j]=add_to_basis_ideal_quotient(p,-1,-1,c,ibuf+j);
    2739    
    2740     }
    2741     int sum=0;
    2742     for(j=0;j<F_plus_index;j++)
    2743     {
    2744       sum+=ibuf[j];
    2745     }
    2746     sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*));
    2747     int partsum=0;
    2748     for(j=0;j<F_plus_index;j++)
    2749     {
    2750       memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*));
    2751       omfree(sbuf[j]);
    2752       partsum+=ibuf[j];
    2753     }
    2754 
    2755     qsort(big_sbuf,sum,sizeof(sorted_pair_node*),pair_better_gen2);
    2756     c->apairs=merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
    2757     c->pair_top+=sum;
    2758     clean_top_of_pair_list(c);
    2759     omfree(big_sbuf);
    2760     omfree(sbuf);
    2761     omfree(ibuf);
    2762   }
    2763   omfree(F_plus);
    2764 #ifdef TGB_DEBUG
    2765   int z;
    2766   for(z=1;z<=c->pair_top;z++)
    2767   {
    2768     assume(pair_better(c->apairs[z],c->apairs[z-1],c));
    2769   }
    2770 #endif
    2771 
    2772   return;
    2773 }
     1915
     1916
     1917
     1918
     1919
     1920
    27741921
    27751922
    … …  
    30942241  return(I);
    30952242}
    3096 static void now_t_rep(const int & arg_i, const int & arg_j, calc_dat* c){
     2243void now_t_rep(const int & arg_i, const int & arg_j, calc_dat* c){
    30972244  int i,j;
    30982245  if (arg_i==arg_j){
    … …  
    32332380  else return (c->apairs[c->pair_top--]);
    32342381}
    3235 static sorted_pair_node* top_pair(calc_dat* c){
     2382sorted_pair_node* top_pair(calc_dat* c){
    32362383  super_clean_top_of_pair_list(c);//yeah, I know, it's odd that I use a different proc here
    32372384
    … …  
    32392386  else return (c->apairs[c->pair_top]);
    32402387}
    3241 static sorted_pair_node* quick_pop_pair(calc_dat* c){
     2388sorted_pair_node* quick_pop_pair(calc_dat* c){
    32422389  if(c->pair_top<0) return NULL;
    32432390  else return (c->apairs[c->pair_top--]);
    … …  
    32602407  }
    32612408}
    3262 static void clean_top_of_pair_list(calc_dat* c){
     2409void clean_top_of_pair_list(calc_dat* c){
    32632410  while((c->pair_top>=0) && (c->apairs[c->pair_top]->i>=0) && (!state_is(UNCALCULATED,c->apairs[c->pair_top]->j, c->apairs[c->pair_top]->i,c))){
    32642411
    … …  
    32832430  else return(c->states[arg_j][arg_i]==state);
    32842431}
    3285 static void free_sorted_pair_node(sorted_pair_node* s, ring r){
     2432void free_sorted_pair_node(sorted_pair_node* s, ring r){
    32862433  if (s->i>=0)
    32872434    p_Delete(&s->lcm_of_lm,r);
    … …  
    33042451}
    33052452
    3306 static int pair_better_gen(const void* ap,const void* bp){
     2453static int tgb_pair_better_gen(const void* ap,const void* bp){
    33072454
    33082455  sorted_pair_node* a=*((sorted_pair_node**)ap);

  • kernel/tgb_internal.h

    r1577276 rb75d13  
    55*  Computer Algebra System SINGULAR     *
    66****************************************/
    7 /* $Id: tgb_internal.h,v 1.9 2005-05-12 08:25:08 bricken Exp $ */
     7/* $Id: tgb_internal.h,v 1.10 2005-05-12 09:19:43 bricken Exp $ */
    88/*
    99 * ABSTRACT: tgb internal .h file
    … …  
    154154static poly redNFTail (poly h,const int sl,kStrategy strat, int len);
    155155static poly redNF2 (poly h,calc_dat* c , int &len, number&  m,int n=0);
    156 static void free_sorted_pair_node(sorted_pair_node* s, ring r);
     156void free_sorted_pair_node(sorted_pair_node* s, ring r);
    157157static void shorten_tails(calc_dat* c, poly monom);
    158158static void replace_pair(int & i, int & j, calc_dat* c);
    … …  
    163163static int* make_connections(int from, poly bound, calc_dat* c);
    164164static int* make_connections(int from, int to, poly bound, calc_dat* c);
    165 static void now_t_rep(const int & arg_i, const int & arg_j, calc_dat* c);
     165void now_t_rep(const int & arg_i, const int & arg_j, calc_dat* c);
    166166static void soon_t_rep(const int & arg_i, const int & arg_j, calc_dat* c);
    167167static int pLcmDeg(poly a, poly b);
    … …  
    171171static sorted_pair_node* pop_pair(calc_dat* c);
    172172static BOOLEAN no_pairs(calc_dat* c);
    173 static void clean_top_of_pair_list(calc_dat* c);
     173void clean_top_of_pair_list(calc_dat* c);
    174174static void super_clean_top_of_pair_list(calc_dat* c);
    175175static BOOLEAN state_is(calc_state state, const int & i, const int & j, calc_dat* c);
    176176static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, calc_dat* c);
    177 static int pair_better_gen(const void* ap,const void* bp);
     177static int tgb_pair_better_gen(const void* ap,const void* bp);
    178178static poly redTailShort(poly h, kStrategy strat);
    179179static poly gcd_of_terms(poly p, ring r);
    … …  
    181181static poly kBucketGcd(kBucket* b, ring r);
    182182static void multi_reduction(red_object* los, int & losl, calc_dat* c);
    183 static sorted_pair_node* quick_pop_pair(calc_dat* c);
    184 static sorted_pair_node* top_pair(calc_dat* c);
     183sorted_pair_node* quick_pop_pair(calc_dat* c);
     184sorted_pair_node* top_pair(calc_dat* c);
     185sorted_pair_node** add_to_basis_ideal_quotient(poly h, int i_pos, int j_pos,calc_dat* c, int* ip);
     186sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,calc_dat* c);
     187int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj);
     188int tgb_pair_better_gen2(const void* ap,const void* bp);
     189int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev);
    185190//static int quality(poly p, int len, calc_dat* c);
    186191/**
Note: See TracChangeset for help on using the changeset viewer.