Changeset 4b1d39 in git


Ignore:
Timestamp:
Jan 20, 2002, 12:44:48 PM (22 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'd1b01e9d51ade4b46b745d3bada5c5f3696be3a8')
Children:
d704faffc1735628a9173c79d65d0b7331bcc012
Parents:
b8490e74a539f406ed4238becdea28121774e338
Message:
*hannes/bricken: unified source, cleanup


git-svn-id: file:///usr/local/Singular/svn/trunk@5773 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
2 deleted
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/extra.cc

    rb8490e7 r4b1d39  
    22*  Computer Algebra System SINGULAR      *
    33*****************************************/
    4 /* $Id: extra.cc,v 1.173 2002-01-19 17:20:58 obachman Exp $ */
     4/* $Id: extra.cc,v 1.174 2002-01-20 11:44:47 Singular Exp $ */
    55/*
    66* ABSTRACT: general interface to internals of Singular ("system" command)
    … …  
    8888#include "silink.h"
    8989#include "walk.h"
     90
     91#include "fast_maps.h"
    9092
    9193/*
    … …  
    651653
    652654#include "mod_raw.h"
    653 // #include "fast_maps.cc"
    654655   
    655656static BOOLEAN jjEXTENDED_SYSTEM(leftv res, leftv h)
    … …  
    783784    }
    784785    else
    785 /*==================== ring debug ==================================*/
    786     if(strcmp(sys_cmd,"map")==0)
    787     {
    788       ring image_r = currRing;
    789       map theMap = (map)h->Data();
    790       ideal image_id = (ideal) theMap;
    791       ring map_r = IDRING(idroot->get(theMap->preimage, myynest));
    792       ideal map_id = IDIDEAL(map_r->idroot->get(h->Next()->Name(), myynest));
    793 #if 0
    794       ring src_r, dest_r;
    795       maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r);
    796       mapoly mp;
    797       maideal mideal;
    798          
    799       maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);
    800       maPoly_Out(mp, src_r);
    801 #endif
    802       return FALSE;
    803     }
    804     else
    805786#endif
    806787/*==================== mtrack ==================================*/

  • Singular/fast_maps.cc

    rb8490e7 r4b1d39  
    77 *  Author:  obachman (Olaf Bachmann)
    88 *  Created: 02/01
    9  *  Version: $Id: fast_maps.cc,v 1.27 2002-01-20 10:01:48 Singular Exp $
     9 *  Version: $Id: fast_maps.cc,v 1.28 2002-01-20 11:44:48 Singular Exp $
    1010 *******************************************************************/
    1111#include "mod2.h"
    … …  
    3030*/
    3131// return maximal monomial if max_map_monomials are substituted into pi_m
    32 static poly maGetMaxExpP(poly* max_map_monomials, 
    33                          int n_max_map_monomials, ring map_r, 
     32static poly maGetMaxExpP(poly* max_map_monomials,
     33                         int n_max_map_monomials, ring map_r,
    3434                         poly pi_m, ring pi_r)
    3535{
    … …  
    3838  Exponent_t e_i, e_j;
    3939  poly m_i, map_j = p_Init(map_r);
    40  
     40
    4141  for (i=1; i <= n; i++)
    4242  {
    … …  
    5959
    6060// returns maximal exponent if map_id is applied to pi_id
    61 static Exponent_t maGetMaxExp(ideal map_id, ring map_r, ideal pi_id, ring pi_r)
     61static Exponent_t maGetMaxExp(ideal pi_id, ring pi_r, ideal map_id, ring map_r)
    6262{
    6363  Exponent_t max=0;
    6464  poly* max_map_monomials = (poly*) omAlloc(IDELEMS(map_id)*sizeof(poly));
    6565  poly max_pi_i, max_map_i;
    66  
     66
    6767  int i, j;
    6868  for (i=0; i<IDELEMS(map_id); i++)
    … …  
    7070    max_map_monomials[i] = p_GetMaxExpP(map_id->m[i], map_r);
    7171  }
    72  
     72
    7373  for (i=0; i<IDELEMS(pi_id); i++)
    7474  {
    7575    max_pi_i = p_GetMaxExpP(pi_id->m[i], pi_r);
    76     max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r, 
     76    max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r,
    7777                              max_pi_i, pi_r);
    7878    Exponent_t temp = p_GetMaxExp(max_map_i, map_r);
    … …  
    102102    p_wrp(monomial->dest, dest_r);
    103103  }
    104   if (monomial->f1!=NULL) { printf(" f1:%x", (unsigned int)monomial->f1->src); 
     104  if (monomial->f1!=NULL) { printf(" f1:%x", (unsigned int)monomial->f1->src);
    105105                            // p_wrp(monomial->f1->src, src_r);
    106                           }
    107   if (monomial->f2!=NULL) { printf(" f2:%x",(unsigned int)monomial->f2->src); 
     106                          }
     107  if (monomial->f2!=NULL) { printf(" f2:%x",(unsigned int)monomial->f2->src);
    108108                            // p_wrp(monomial->f2->src, src_r);
    109                           }
     109                          }
    110110  printf("\n");
    111111  fflush(stdout);
    … …  
    161161      }
    162162      while (next != NULL);
    163       if (mp->dest != NULL) 
     163      if (mp->dest != NULL)
    164164      {
    165165        assume(dest_r != NULL);
    … …  
    182182    return what;
    183183  }
    184  
     184
    185185  mapoly iter = into;
    186186  mapoly prev = NULL;
    187  
     187
    188188  Top:
    189189  p_LmCmpAction(iter->src, what->src, src_r, goto Equal, goto Greater, goto Smaller);
    190  
     190
    191191  Greater:
    192192  if (iter->next == NULL)
    … …  
    198198  iter = iter->next;
    199199  goto Top;
    200  
     200
    201201  Smaller:
    202202  if (prev == NULL)
    … …  
    209209  what->next = iter;
    210210  return what;
    211  
     211
    212212  Equal:
    213213  iter->ref += what->ref;
    214214  macoeff coeff = what->coeff;
    215   if (coeff != NULL) 
     215  if (coeff != NULL)
    216216  {
    217217    while (coeff->next != NULL) coeff = coeff->next;
    … …  
    232232{
    233233  poly next;
    234  
     234
    235235  while (what != NULL)
    236236  {
    … …  
    254254  int i;
    255255  mp = NULL;
    256  
     256
    257257  for (i=0; i<mideal->n; i++)
    258258  {
    … …  
    260260    {
    261261      mideal->buckets[i] = sBucketCreate(dest_r);
    262       maPoly_InsertPoly(mp, 
     262      maPoly_InsertPoly(mp,
    263263#ifdef PDEBUG
    264264                        prShallowCopyR(map_id->m[i], map_r, src_r),
    … …  
    266266                        prShallowCopyR_NoSort(map_id->m[i], map_r, src_r),
    267267#endif
    268                         src_r,                     
     268                        src_r,
    269269                        mideal->buckets[i]);
    270270    }
    … …  
    273273
    274274
    275 void maMap_CreateRings(ideal map_id, ring map_r, 
    276                        ideal image_id, ring image_r, 
     275void maMap_CreateRings(ideal map_id, ring map_r,
     276                       ideal image_id, ring image_r,
    277277                       ring &src_r, ring &dest_r, BOOLEAN &simple)
    278278{
    … …  
    294294  Exponent_t maxExp = maGetMaxExp(map_id, map_r, image_id, image_r);
    295295  if (maxExp <=  1) maxExp = 2;
    296   else if (maxExp > (Exponent_t) image_r->bitmask) 
     296  else if (maxExp > (Exponent_t) image_r->bitmask)
    297297    maxExp = (Exponent_t) image_r->bitmask;
    298298  dest_r = rModifyRing_Simple(image_r, TRUE, TRUE, maxExp,  simple);
    … …  
    314314*F  misc  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  misc  stuff
    315315*/
    316    
     316
    317317ideal maIdeal_2_Ideal(maideal m_id, ring dest_r)
    318318{
    319319  ideal res = idInit(m_id->n, 1);
    320320  int l;
    321  
     321
    322322  for (int i= 0; i < m_id->n; i++)
    323323  {
    … …  
    338338}
    339339
    340  
     340
    341341/*******************************************************************************
    342342**
    … …  
    351351  BOOLEAN no_sort;
    352352
    353   // construct rings we work in: 
     353  // construct rings we work in:
    354354  // src_r: Wp with Weights set to length of poly in image_id
    355355  // dest_r: Simple ring without degree ordering and short exponents
    356356  maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r, no_sort);
    357  
     357
    358358  // construct dest_id
    359359  if (dest_r != image_r)
    … …  
    366366  maideal mideal;
    367367  maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);
    368  
     368
    369369  if (TEST_OPT_PROT)
    370370  {
    … …  
    372372    Print("map[%d:%d]{%d:", dest_r->bitmask, dest_r->ExpL_Size, length);
    373373  }
    374  
     374
    375375  // do the optimization step
    376376#if HAVE_MAP_OPTIMIZE > 0
    … …  
    382382    Print("%d}", length);
    383383  }
    384      
     384
    385385  // do the actual evaluation
    386386  maPoly_Eval(mp, src_r, dest_id, dest_r, length);
    387387  if (TEST_OPT_PROT) Print(".");
    388  
     388
    389389  // collect the results back into an ideal
    390390  ideal res_dest_id = maIdeal_2_Ideal(mideal, dest_r);
    … …  
    397397    if (no_sort)
    398398      res_image_id = idrShallowCopyR_NoSort(res_dest_id, dest_r, image_r);
    399     else 
     399    else
    400400      res_image_id = idrShallowCopyR(res_dest_id, dest_r, image_r);
    401401    id_ShallowDelete(&res_dest_id, dest_r);
    … …  
    403403  else
    404404    res_image_id = res_dest_id;
    405  
     405
    406406  if (TEST_OPT_PROT) Print(".");
    407407  // clean-up the rings
    … …  
    415415
    416416
    417 
    418 
    419    
    420    
    421                    
    422  
    423    
    424 
    425 
    426 
    427 
    428 
    429 
    430 
    431 
    432  
    433  
    434  
    435  
    436 
    437  
    438      
     417/**********************************************************************
     418* Evaluation stuff                                                    *
     419**********************************************************************/
     420
     421// substitute p everywhere the monomial occours,
     422// return the number of substitutions
     423static int maPoly_Substitute(macoeff c, poly p, ring dest_r)
     424{
     425  // substitute the monomial: go through macoeff
     426  int len=pLength(p);
     427  int done=0;
     428  while (c!=NULL)
     429  {
     430    done++;
     431    poly t=pp_Mult_nn(p,c->n,dest_r);
     432    sBucket_Add_p(c->bucket, t, len);
     433    c=c->next;
     434  }
     435  return done;
     436}
     437
     438static poly maPoly_EvalMon(poly src, ring src_r, poly* dest_id, ring dest_r)
     439{
     440  int i;
     441  int e;
     442  poly p=NULL;
     443  poly pp;
     444  for(i=1;i<=src_r->N;i++)
     445  {
     446    e=p_GetExp(src,i,src_r);
     447    if (e>0)
     448    {
     449      pp=dest_id[i-1];
     450      if (pp==NULL)
     451      {
     452        p_Delete(&p,dest_r);
     453        return NULL;
     454      }
     455      if ((p==NULL) /* && (e>0)*/)
     456      {
     457        p=p_Copy(pp /*dest_id[i-1]*/,dest_r);
     458        e--;
     459      }
     460      while (e>0)
     461      {
     462        p=p_Mult_q(p,p_Copy(pp /*dest_id[i-1]*/,dest_r),dest_r);
     463        e--;
     464      }
     465    }
     466  }
     467  if (p==NULL) p=p_ISet(1,dest_r);
     468  return p;
     469}
     470
     471void maPoly_Eval(mapoly root, ring src_r, ideal dest_id, ring dest_r, int total_cost)
     472{
     473  // invert the list rooted at root:
     474  if ((root!=NULL) && (root->next!=NULL))
     475  {
     476    mapoly q=root->next;
     477    mapoly qn;
     478    root->next=NULL;
     479    do
     480    {
     481      qn=q->next;
     482      q->next=root;
     483      root=q;
     484      q=qn;
     485    }
     486    while (qn !=NULL);
     487  }
     488
     489  total_cost /= 10;
     490  int next_print_cost = total_cost;
     491
     492  // the evaluation -----------------------------------------
     493  mapoly p=root;
     494  int cost = 0;
     495
     496  while (p!=NULL)
     497  {
     498    // look at each mapoly: compute its value in ->dest
     499    assume (p->dest==NULL);
     500    {
     501      if ((p->f1!=NULL)&&(p->f2!=NULL))
     502      {
     503#if 0
     504        printf("found prod:");
     505        p_wrp(p->src,src_r);printf("=");
     506        p_wrp(p->f1->src,src_r);printf(" * ");
     507        p_wrp(p->f2->src,src_r);printf("\n");
     508#endif
     509        poly f1=p->f1->dest;
     510        poly f2=p->f2->dest;
     511        if (p->f1->ref>0) f1=p_Copy(f1,dest_r);
     512        else
     513        {
     514          // we own p->f1->dest now (in f1)
     515          p->f1->dest=NULL;
     516        }
     517        if (p->f2->ref>0) f2=p_Copy(f2,dest_r);
     518        else
     519        {
     520          // we own p->f2->dest now (in f2)
     521          p->f2->dest=NULL;
     522        }
     523        maMonomial_Free(p->f1,src_r, dest_r);
     524        maMonomial_Free(p->f2,src_r, dest_r);
     525        p->dest=p_Mult_q(f1,f2,dest_r);
     526      } /* factors : 2 */
     527      else
     528      {
     529        //printf("compute "); p_wrp(p->src,src_r);printf("\n");
     530        assume((p->f1==NULL) && (p->f2==NULL));
     531        //if(p->f1!=NULL) { printf(" but f1="); p_wrp(p->f1->src,src_r);printf("\n"); }
     532        //if(p->f2!=NULL) { printf(" but f2="); p_wrp(p->f2->src,src_r);printf("\n"); }
     533        // no factorization provided, use the classical method:
     534        p->dest=maPoly_EvalMon(p->src,src_r,dest_id->m,dest_r);
     535        //p_wrp(p->dest, dest_r); printf(" is p->dest\n");
     536      }
     537    } /* p->dest==NULL */
     538    // substitute the monomial: go through macoeff
     539    p->ref -= maPoly_Substitute(p->coeff, p->dest, dest_r);
     540    //printf("subst done\n");
     541    if (total_cost)
     542    {
     543      assume(TEST_OPT_PROT);
     544      cost++;
     545      if (cost > next_print_cost)
     546      {
     547        Print("-");
     548        next_print_cost += total_cost;
     549      }
     550    }
     551
     552    mapoly pp=p;
     553    p=p->next;
     554    //p_wrp(pp->src, src_r);
     555    if (pp->ref<=0)
     556    {
     557      //printf(" (%x) killed\n",pp);
     558      maMonomial_Destroy(pp, src_r, dest_r);
     559    }
     560    //else
     561    //  printf(" (%x) not killed, ref=%d\n",pp,pp->ref);
     562  }
     563}
     564
     565
     566/*******************************************************************************
     567**
     568*F  maEggt  . . . . . . . . . . . . . . . . . . . . . . . .  returns extended ggt
     569*/
     570// return NULL if deg(ggt(m1, m2)) < 2
     571// else return m = ggT(m1, m2) and q1, q2 such that m1 = q1*m m2 = q2*m
     572static poly maEggT(const poly m1, const poly m2, poly &q1, poly &q2,const ring r)
     573{
     574
     575  int i;
     576  int dg = 0;
     577  poly ggt = NULL;
     578  q1 = p_Init(r);
     579  q2 = p_Init(r);
     580  ggt=p_Init(r);
     581
     582  for (i=1;i<=r->N;i++) {
     583    Exponent_t e1 = p_GetExp(m1, i, r);
     584    Exponent_t e2 = p_GetExp(m2, i, r);
     585    if (e1 > 0 && e2 > 0){
     586      Exponent_t em = (e1 > e2 ? e2 : e1);
     587      dg += em;
     588      p_SetExp(ggt, i, em, r);
     589      p_SetExp(q1, i, e1 - em, r);
     590      p_SetExp(q2, i, e2 - em, r);
     591    }
     592    else {
     593      p_SetExp(q1, i, e1, r);
     594      p_SetExp(q2, i, e2, r);
     595    }
     596  }
     597  if (dg>1)
     598  {
     599    p_Setm(ggt, r);
     600    p_Setm(q1, r);
     601    p_Setm(q2, r);
     602
     603
     604  }
     605  else {
     606    p_LmFree(ggt, r);
     607    p_LmFree(q1, r);
     608    p_LmFree(q2, r);
     609    ggt = NULL;
     610  }
     611  return ggt;
     612}
     613
     614/********************************************************************
     615 **                                                                 *
     616 * maFindBestggT                                                    *
     617 * finds ggT with the highest cost                                  *
     618 *******************************************************************/
     619
     620static mapoly maFindBestggT(mapoly mp, mapoly & choice, mapoly & fp, mapoly & fq,const ring r)
     621{
     622  int ggt_deg = 0;
     623  poly p = mp->src;
     624  mapoly iter = choice;
     625  poly ggT = NULL;
     626  fp = NULL;
     627  fq = NULL;
     628  poly fp_p=NULL;
     629  poly fq_p=NULL;
     630  choice=NULL;
     631  while ((iter != NULL) && (pDeg(iter->src, r) > ggt_deg))
     632  {
     633    //    maMonomial_Out(iter, r, NULL);
     634    poly q1, q2, q;
     635
     636    q = maEggT(p, iter->src, q1, q2,r);
     637    if (q != NULL)
     638    {
     639      assume((q1!=NULL)&&(q2!=NULL));
     640      if (pDeg(q,r) > ggt_deg)
     641      {
     642        choice=iter;
     643        if (ggT != NULL)
     644        {
     645          p_LmFree(ggT, r);
     646          p_LmFree(fp_p, r);
     647          p_LmFree(fq_p, r);
     648
     649        }
     650        ggt_deg = pDeg(q, r);
     651        ggT = q;
     652        fp_p = q1;
     653        fq_p = q2;
     654      }
     655      else
     656      {
     657        p_LmFree(q, r);
     658        p_LmFree(q1, r);
     659        p_LmFree(q2, r);
     660      }
     661    }
     662    iter=iter->next;
     663  }
     664  if(ggT!=NULL){
     665
     666    int dq =pTotaldegree(fq_p,r);
     667    if (dq!=0){
     668      fq=maPoly_InsertMonomial(mp, fq_p, r, NULL);
     669      fp=maPoly_InsertMonomial(mp, fp_p, r, NULL);
     670      return maPoly_InsertMonomial(mp, ggT, r, NULL);
     671    }
     672    else {
     673      fq=NULL;
     674      p_LmFree(fq_p, r);
     675      p_LmFree(ggT, r);
     676      fp=maPoly_InsertMonomial(mp, fp_p, r, NULL);
     677      choice->ref++;
     678      return choice;
     679    }
     680  }
     681  else {
     682    return NULL;
     683  }
     684
     685
     686}
     687
     688/********************************************************************
     689 **                                                                 *
     690 * maPoly_Optimize                                                  *
     691 * adds and integrates subexpressions                               *
     692 *******************************************************************/
     693
     694void maPoly_Optimize(mapoly mpoly, ring src_r){
     695  assume(mpoly!=NULL && mpoly->src!=NULL);
     696  mapoly iter = mpoly;
     697  mapoly choice;
     698  mapoly ggT=NULL;
     699  mapoly fp=NULL;
     700  mapoly fq=NULL;
     701  while (iter->next!=NULL){
     702    choice=iter->next;
     703    if ((iter->f1==NULL)){
     704      ggT=maFindBestggT(iter, choice, fp, fq,src_r);
     705      if (choice!=NULL){
     706        iter->f1=fp;
     707        iter->f2=ggT;
     708        if (fq!=NULL){
     709          ggT->ref++;
     710          choice->f1=fq;
     711          choice->f2=ggT;
     712        }
     713      }
     714
     715    }
     716    iter=iter->next;
     717  }
     718}
Note: See TracChangeset for help on using the changeset viewer.