Changeset 1873199 in git for Singular


Ignore:
Timestamp:
Jan 19, 2002, 7:03:38 PM (22 years ago)
Author:
Olaf Bachmann <obachman@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
Children:
d70379bc408f4461637440a4488dcfdb4e84a4ff
Parents:
2e2ba930ca93281e4952f68c62d045ea28a4ab8d
Message:
added OPT_PROT stuff


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

Legend:

Unmodified
Added
Removed

  • Singular/fast_eval.cc

    r2e2ba9 r1873199  
    33#include "p_polys.h"
    44#include "fast_maps.h"
     5#include "febase.h"
    56
    67// substitute p everywhere the monomial occours,
    … …  
    4445}
    4546
    46 void maPoly_Eval(mapoly root, ring src_r, ideal dest_id, ring dest_r)
     47void maPoly_Eval(mapoly root, ring src_r, ideal dest_id, ring dest_r, int total_cost)
    4748{
    4849  // invert the list rooted at root:
    … …  
    6263  }
    6364
     65  total_cost /= 10;
     66  int next_print_cost = total_cost;
     67 
    6468  // the evaluation -----------------------------------------
    6569  mapoly p=root;
     70  int cost = 0;
     71 
    6672  while (p!=NULL)
    6773  {
    68      // look at each mapoly: compute its value in ->dest
    69      assume (p->dest==NULL);
    70      {
    71        if ((p->f1!=NULL)&&(p->f2!=NULL))
    72        {
    73          printf("found prod:");
    74          p_wrp(p->src,src_r);printf("=");
    75          p_wrp(p->f1->src,src_r);printf(" * ");
    76          p_wrp(p->f2->src,src_r);printf("\n");
    77          poly f1=p->f1->dest;
    78          poly f2=p->f2->dest;
    79          if (p->f1->ref>0) f1=p_Copy(f1,dest_r);
    80          else
    81          {
    82            // we own p->f1->dest now (in f1)
    83            p->f1->dest=NULL;
    84          }
    85          if (p->f2->ref>0) f2=p_Copy(f2,dest_r);
    86          else
    87          {
    88            // we own p->f2->dest now (in f2)
    89            p->f2->dest=NULL;
    90          }
    91          maMonomial_Free(p->f1,src_r, dest_r);
    92          maMonomial_Free(p->f2,src_r, dest_r);
    93          p->dest=p_Mult_q(f1,f2,dest_r);
    94        } /* factors : 2 */
    95        else
    96        {
    97          //printf("compute "); p_wrp(p->src,src_r);printf("\n");
    98          assume((p->f1==NULL) && (p->f2==NULL));
    99          //if(p->f1!=NULL) { printf(" but f1="); p_wrp(p->f1->src,src_r);printf("\n"); }
    100          //if(p->f2!=NULL) { printf(" but f2="); p_wrp(p->f2->src,src_r);printf("\n"); }
    101          // no factorization provided, use the classical method:
    102          p->dest=maPoly_EvalMon(p->src,src_r,dest_id->m,dest_r);
    103          //p_wrp(p->dest, dest_r); printf(" is p->dest\n");
    104        }
    105      } /* p->dest==NULL */
    106      // substitute the monomial: go through macoeff
    107      p->ref -= maPoly_Substitute(p->coeff, p->dest, dest_r);
    108      //printf("subst done\n");
    109      mapoly pp=p;
    110      p=p->next;
    111      if (pp->ref<=0)
    112      {
    113        maMonomial_Destroy(pp, src_r, dest_r);
    114      }
    115    }
     74    // look at each mapoly: compute its value in ->dest
     75    assume (p->dest==NULL);
     76    {
     77      if ((p->f1!=NULL)&&(p->f2!=NULL))
     78      {
     79        printf("found prod:");
     80        p_wrp(p->src,src_r);printf("=");
     81        p_wrp(p->f1->src,src_r);printf(" * ");
     82        p_wrp(p->f2->src,src_r);printf("\n");
     83        poly f1=p->f1->dest;
     84        poly f2=p->f2->dest;
     85        if (p->f1->ref>0) f1=p_Copy(f1,dest_r);
     86        else
     87        {
     88          // we own p->f1->dest now (in f1)
     89          p->f1->dest=NULL;
     90        }
     91        if (p->f2->ref>0) f2=p_Copy(f2,dest_r);
     92        else
     93        {
     94          // we own p->f2->dest now (in f2)
     95          p->f2->dest=NULL;
     96        }
     97        maMonomial_Free(p->f1,src_r, dest_r);
     98        maMonomial_Free(p->f2,src_r, dest_r);
     99        p->dest=p_Mult_q(f1,f2,dest_r);
     100      } /* factors : 2 */
     101      else
     102      {
     103        //printf("compute "); p_wrp(p->src,src_r);printf("\n");
     104        assume((p->f1==NULL) && (p->f2==NULL));
     105        //if(p->f1!=NULL) { printf(" but f1="); p_wrp(p->f1->src,src_r);printf("\n"); }
     106        //if(p->f2!=NULL) { printf(" but f2="); p_wrp(p->f2->src,src_r);printf("\n"); }
     107        // no factorization provided, use the classical method:
     108        p->dest=maPoly_EvalMon(p->src,src_r,dest_id->m,dest_r);
     109        //p_wrp(p->dest, dest_r); printf(" is p->dest\n");
     110      }
     111    } /* p->dest==NULL */
     112    // substitute the monomial: go through macoeff
     113    p->ref -= maPoly_Substitute(p->coeff, p->dest, dest_r);
     114    //printf("subst done\n");
     115    if (total_cost)
     116    {
     117      assume(TEST_OPT_PROT);
     118      cost += pDeg(p->src, src_r);
     119      if (cost > next_print_cost)
     120      {
     121        Print("-");
     122        next_print_cost += total_cost;
     123      }
     124    }
     125
     126    mapoly pp=p;
     127    p=p->next;
     128    if (pp->ref<=0)
     129    {
     130      maMonomial_Destroy(pp, src_r, dest_r);
     131    }
     132  }
    116133}

  • Singular/fast_maps.cc

    r2e2ba9 r1873199  
    77 *  Author:  obachman (Olaf Bachmann)
    88 *  Created: 02/01
    9  *  Version: $Id: fast_maps.cc,v 1.18 2002-01-19 17:55:39 Singular Exp $
     9 *  Version: $Id: fast_maps.cc,v 1.19 2002-01-19 18:03:38 obachman Exp $
    1010 *******************************************************************/
    1111#include "mod2.h"
    … …  
    2727/*******************************************************************************
    2828**
    29 *F  maMaxExp  . . . . . . . . . . . .  returns maximal exponent of result of map
     29*F  maMaxExp  . . . . . . . .  returns bound on maximal exponent of result of map
    3030*/
    3131// return maximal monomial if max_map_monomials are substituted into pi_m
    … …  
    8888
    8989
    90 #if 0
    91 // paste into extra.cc
    92     if(strcmp(sys_cmd,"map")==0)
    93     {
    94       ring image_r = currRing;
    95       map theMap = (map)h->Data();
    96       ideal image_id = (ideal) theMap;
    97       ring map_r = IDRING(idroot->get(theMap->preimage, myynest));
    98       ideal map_id = IDIDEAL(map_r->idroot->get(h->Next()->Name(), myynest));
    99 
    100       ring src_r, dest_r;
    101       maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r);
    102       mapoly mp;
    103       maideal mideal;
    104          
    105       maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);
    106       maPoly_Out(mp, src_r);
    107       return FALSE;
    108     }
    109 
    110 // use as Singular source template
    111 ring map_r = 32003, (a, b, c), Dp;
    112 ideal map_id = a, ab, c3 + ab, a2b + ab;
    113 
    114 
    115 ring image_r;
    116 ideal image_id = x2y3, y, z+x6;
    117 
    118 map Phi = map_r, image_id;
    119 
    120 
    121 system("map", Phi, map_id);
    122 
    123 #endif
    124 
    12590/*******************************************************************************
    12691**
    12792*F  debugging stuff
    12893*/
     94#ifndef NDEBUG
    12995void maMonomial_Out(mapoly monomial, ring src_r, ring dest_r = NULL)
    13096{
    … …  
    150116  }
    151117}
    152 
     118#endif
    153119
    154120/*******************************************************************************
    … …  
    271237}
    272238
     239/*******************************************************************************
     240**
     241*F  maMap_Create . . . . . . . . . . . . . . . . . . . . create stuff
     242*/
     243
    273244void maMap_CreatePolyIdeal(ideal map_id, ring map_r, ring src_r, ring dest_r,
    274245                           mapoly &mp, maideal &mideal)
    … …  
    299270{
    300271#if HAVE_SRC_R > 0
     272  int* weights = (int*) omAlloc0(map_r->N);
     273  int i;
     274  int n = min(map_r->N, IDELEMS(theMap));
     275
     276  for (i=0; i<n; i++)
     277  {
     278    weights[i] = pLength(theMap->m[i]);
     279  }
     280  src_r = rModifyRing_Wp(map_r, weights);
    301281#else
    302282  src_r = map_r;
    … …  
    311291}
    312292
    313 void maMap_KillRings(ring map_r, ring image_r, ring src_r, ring dest_r)
     293static void maMap_KillRings(ring map_r, ring image_r, ring src_r, ring dest_r)
    314294{
    315295  if (map_r != src_r)
    … …  
    318298    rKillModifiedRing_Simple(dest_r);
    319299}
     300
     301/*******************************************************************************
     302**
     303*F  misc  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  misc  stuff
     304*/
    320305   
    321306ideal maIdeal_2_Ideal(maideal m_id, ring dest_r)
    … …  
    332317}
    333318
     319void maPoly_GetLength(mapoly mp, ring src_r, int length, int total_cost)
     320{
     321  length = 0;
     322  total_cost = 0;
     323  while (mp != NULL)
     324  {
     325    length++;
     326    if (mp->src != NULL)
     327      total_cost += pDeg(mp->src, src_r);
     328    mp = mp->next;
     329  }
     330}
     331
     332 
     333/*******************************************************************************
     334**
     335*F  fast_map  . . . . . . . . . . . . . . . . . . . . . . . . . .the real thing
     336*/
    334337
    335338ideal fast_map(ideal map_id, ring map_r, ideal image_id, ring image_r)
    … …  
    337340  ring src_r, dest_r;
    338341  ideal dest_id, res_id;
     342  int length, total_cost = 0;
    339343
    340344  // construct rings we work in:
    … …  
    354358  maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);
    355359 
     360  if (TEST_OPT_PROT)
     361  {
     362    maPoly_GetLength(mp, src_r, length, total_cost);
     363    Print("map[%d:%d]{%d:%d}", src_r->bitmask, src_r->ExpL_Size, length, total_cost);
     364  }
     365 
    356366  // do the optimization step
    357367#if HAVE_MAP_OPTIMIZE > 0
    358368  maPoly_Optimize(mp, src_r);
    359369#endif
    360 
     370  if (TEST_OPT_PROT)
     371  {
     372    maPoly_GetLength(mp, src_r, length, total_cost);
     373    Print("{%d:%d}", length, total_cost);
     374  }
     375     
    361376  // do the actual evaluation
    362   maPoly_Eval(mp, src_r, dest_id, dest_r);
     377  maPoly_Eval(mp, src_r, dest_id, dest_r, total_cost);
    363378
    364379  // collect the results back into an ideal
    … …  
    382397
    383398
    384 #if 0
    385 
    386 /*******************************************************************************
    387 **
    388 *F  maMaxExp  . . . . . . . . . . . .  returns maximal exponent of result of map
    389 */
    390 
    391 // return maximal monomial if max_map_monomials are substituted into pi_m
    392 static poly maGetMaxExpP(poly* max_map_monomials,
    393                          int n_max_map_monomials, ring map_r,
    394                          poly pi_m, ring pi_r)
    395 {
    396   int n = min(pi_r->N, n_max_map_monomials);
    397   int i, j;
    398   Exponent_t e_i, e_j;
    399   poly m_i, map_j = p_Init(map_r);
    400  
    401   for (i=1; i <= n; i++)
    402   {
    403     e_i = p_GetExp(pi_m, i, pi_r);
    404     m_i = max_map_monomials[i-1];
    405     if (e_i > 0 && m_i != NULL && ! p_IsConstantComp(m_i, map_r))
    406     {
    407       for (j = 1; j<= map_r->N; j++)
    408       {
    409         e_j = p_GetExp(m_i, j, map_r);
    410         if (e_j > 0)
    411         {
    412           p_AddExp(map_j, j, e_j*e_i, map_r);
    413         }
    414       }
    415     }
    416   }
    417   return map_j;
    418 }
    419 
    420 // returns maximal monomial if map_id is applied to pi_id
    421 static poly maGetMaxExpP(ideal map_id, ring map_r, ideal pi_id, ring pi_r)
    422 {
    423   poly* max_map_monomials = (poly*) omAlloc(IDELEMS(map_id)*sizeof(poly));
    424   poly max_pi_i, max_map_i;
    425   poly max_map = p_Init(map_r);
    426  
    427   int i, j;
    428   for (i=0; i<IDELEMS(map_id); i++)
    429   {
    430     max_map_monomials[i] = p_GetMaxExpP(map_id->m[i], map_r);
    431   }
    432  
    433   for (i=0; i<IDELEMS(pi_id); i++)
    434   {
    435     max_pi_i = p_GetMaxExpP(pi_id->m[i], pi_r);
    436     max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r,
    437                               max_pi_i, pi_r);
    438     // get maximum
    439     for (j = 1; j<= map_r->N; j++)
    440     {
    441       if (p_GetExp(max_map, j, map_r) < p_GetExp(max_map_i, j, map_r))
    442         p_SetExp(max_map, j, p_GetExp(max_map_i, j, map_r), map_r);
    443     }
    444     p_LmFree(max_pi_i, pi_r);
    445     p_LmFree(max_map_i, map_r);
    446   }
    447   return max_map;
    448 }
    449 
    450 // returns maximal exponent if map_id is applied to pi_id
    451 static Exponent_t maGetMaxExp(ideal map_id, ring map_r, ideal pi_id, ring pi_r)
    452 {
    453   poly p = maGetMaxExpP(map_id, map_r, pi_id, pi_r);
    454   Exponent_t max = p_GetMaxExp(p, map_r);
    455   p_LmFree(p, map_r);
    456   return max;
    457 }
    458 
    459 // construct ring/map ideal  in/with which we perform computations
    460 // return TRUE if ordering changed (not yet implemented), false, otherwise
    461 static BOOLEAN maGetCompIdealRing(ideal map_id, ring map_r, ideal pi_id,
    462                                    ring pi_r, ideal &comp_map_id, ring &comp_r)
    463 {
    464   Exponent_t max_exp = maGetMaxExp(map_id, map_r, pi_id, pi_r);
    465 
    466   comp_r = rModifyRing(map_r, FALSE, !idIsModule(pi_id, pi_r), max_exp);
    467   if (comp_r != map_r)
    468   {
    469     if (TEST_OPT_PROT)
    470       Print("[%d:%d]", (long) comp_r->bitmask, comp_r->ExpL_Size);
    471     comp_map_id = idrShallowCopyR_NoSort(map_id, map_r, comp_r);
    472   }
    473   else
    474     comp_map_id = map_id;
    475  
    476   return FALSE;
    477 }
    478 
    479 static void maDestroyCompIdealRing(ideal map_id, ring map_r,
    480                                     ideal comp_map_id, ring &comp_r,
    481                                     ideal &result, BOOLEAN sort=FALSE)
    482 {
    483   if (map_r != comp_r)
    484   {
    485     result = idrMoveR(result, comp_r, map_r);
    486     id_ShallowDelete(&comp_map_id, comp_r);
    487     rKillModifiedRing(comp_r);
    488   }
    489 }
    490 
    491 /*******************************************************************************
    492 **
    493 *F  maEggt  . . . . . . . . . . . . . . . . . . . . . . . .  returns extended ggt
    494 */
    495 // return NULL if deg(ggt(m1, m2)) < 2
    496 // else return m = ggT(m1, m2) and q1, q2 such that m1 = q1*m m2 = q2*m
    497 static poly maEggT(const poly m1, const poly m2, poly &q1, poly &q2,const ring r)
    498 {
    499   int i;
    500   int dg = 0;
    501   poly ggt = NULL;
    502   for (i=1; i<=r->N; i++)
    503   {
    504     Exponent_t e1 = p_GetExp(m1, i, r);
    505     Exponent_t e2 = p_GetExp(m2, i, r);
    506     if (e1 > 0 && e2 > 0)
    507     {
    508       Exponent_t em = (e1 > e2 ? e2 : e1);
    509       if (dg < 2)
    510       {
    511         ggt = p_Init(r);
    512         q1 = p_Init(r);
    513         q2 = p_Init(r);
    514       }
    515       dg += em;
    516       p_SetExp(ggt, i, em, r);
    517       p_SetExp(q1, i, e1 - em, r);
    518       p_SetExp(q2, i, e2 - em, r);
    519     }
    520   }
    521   if (ggt != NULL)
    522   {
    523     p_Setm(ggt, r);
    524     p_Setm(q1, r);
    525     p_Setm(q2, r);
    526   }
    527   return ggt;
    528 }
    529 
    530 /*******************************************************************************
    531 **
    532 *F  maGetMapRing  . . . . . . . . . . . . gets ring and ideal with which we work
    533 */
    534 static ring maGetWeightedRing(ideal theMap, ring map_r) 
    535 {
    536   // we work in a ring with ordering Deg,WeightedDegree,vars
    537   // First, construct weighted degrees
    538   int* weights = (int*) omAlloc0(map_r->N);
    539   int i;
    540   int n = min(map_r->N, IDELEMS(theMap));
    541 
    542   for (i=0; i<n; i++)
    543   {
    544     weights[i] = pLength(theMap->m[i]);
    545   }
    546   return rModifyRing_Wp(map_r, weights);
    547 }
    548 static void maDestroyWeightedRing(ring r)
    549 {
    550   rKillModified_Wp_Ring(r);
    551 }
    552 
    553 
    554 
    555 
    556 
    557 static void maMap_2_maPoly(ideal source_id, ring source_r, ring weight_r, ring comp_r,
    558                            mapoly &mp, maideal &mideal)
    559 {
    560   mideal = (maideal) omAlloc0(sizeof(maideal_s));
    561   mideal->n = IDELEMS(source_id);
    562   mideal->buckets = (sBucket_pt*) omAlloc0(mideal->n*sizeof(sBucket_pt));
    563   int i;
    564   mp = NULL;
    565  
    566   for (i=0; i<mideal->n; i++)
    567   {
    568     if (source_id->m[i] != NULL)
    569     {
    570       mideal->buckets[i] = sBucketCreate(comp_r);
    571       mapInsertPoly(mp,
    572                     prShallowCopyR_NoSort(source_id->m[i], source_r, weight_r),
    573                     weight_r,                     
    574                     mideal->buckets[i]);
    575     }
    576     maPoly_Out(mp, source_r);
    577   }
    578 }
    579 
    580 
    581 
    582 
    583 
    584 static mapoly maFindBestggT(mapoly mp, mapoly in, poly ggT, poly fp, poly fq)
    585 {
    586   int ggt_deg = 0;
    587   poly p = mp->src;
    588   mapoly mq = NULL;
    589  
    590   ggT = NULL;
    591   fp = NULL;
    592   fq = NULL;
    593   while (in != NULL && in->factors > 1 && pFDeg(in->src) > ggt_deg)
    594   {
    595     poly q1, q2, q;
    596 //    q = maEggT(p, in->src, q1, q2);
    597     if (q != NULL)
    598     {
    599       if (pFDeg(q) > fft_deg)
    600       {
    601         if (ggT != NULL)
    602         {
    603           p_LmFree(ggT);
    604           p_LmFree(fp);
    605           p_LmFree(fq);
    606         }
    607         ggT = q;
    608         fp = q1;
    609         fq = q2;
    610       }
    611       else
    612       {
    613         p_LmFree(q);
    614         p_LmFree(q1);
    615         p_LmFree(q2);
    616       }
    617     }
    618   }
    619 }
    620 
    621      
    622 static mapoly maPrepareEval(mapoly mp, ring r)
    623 {
    624   mapoly res = NULL;
    625   mapoly next = mp;
    626  
    627   while (mp != NULL)
    628   {
    629     next = mp->next;
    630     mp->next = res;
    631     res = mp;
    632     if (mp->factors > 1 && mp->f1 == NULL && mp->f2 == NULL)
    633     {
    634       poly fp, fq, ggT;
    635       mapoly mq = maFindBestggT(mp, next, ggT, fp, fq);
    636       if (mq != NULL)
    637       {
    638         assume(fp != NULL);
    639         mp->f1 = maInsertMonomial(next, fp, r, NULL);
    640        
    641         if (ggT != NULL)
    642         {
    643           assume(fq != NULL);
    644           mapoly mggT = maInsertMonomial(next, ggT, r, NULL);
    645           mq->f1 = maInsertMonomial(next, fq, r, NULL);
    646           mq->f2 = mggT;
    647           mp->f2 = mggT;
    648           mggT->ref++;
    649         }
    650         else
    651         {
    652           assume(fq == NULL);
    653           mp->f2 = mq;
    654           mq->ref++;
    655         }
    656       }
    657     }
    658     mp = next;
    659   }
    660   return res;
    661 }
    662 /*******************************************************************************
    663 **
    664 *F  ideal_2_maideal . . . . . . . . . . . . . . . . .  converts ideal to maideal
    665 */
    666 mapoly ideal_2_maideal(ideal id, ring r, maideal mid, ring mr, ring comp_ring,
    667                        nMapFunc nMap)
    668 {
    669   mid->n = IDELEMS(id);
    670   mid->buckets = (sBucket_pt*)omAlloc(mid->n*sizeof(sBucket_pt));
    671   mapoly mpoly = NULL;
    672  
    673   for (int i=0; i<mid->n; i++)
    674   {
    675     mid->buckets[i] = sBucketCreate(comp_ring);
    676     mapoly mpoly_i = poly_2_mapoly(id->m[i], mid, nMap, mid->buckets[i]);
    677     mpoly = mapAdd(mpoly, mpoly_i);
    678   }
    679   return mpoly;
    680 }
    681 
    682 ideal maideal_2_ideal(ideal orig_id,
    683                       maideal mideal, ring comp_ring, ring dest_ring)
    684 {
    685   ideal id = idInit(IDELEMS(orig_id), orig_id->rank);
    686   for (int i=0; i<mid->n; i++)
    687   {
    688     sBucketDestroyAdd(mid->buckets[i], &(id->m[i]), &dummy);
    689   }
    690  
    691   if (comp_ring != dest_ring)
    692     id = idrMoveR_NoSort(id, comp_ring);
    693  
    694   return id;
    695 }
    696 
    697 #endif
    698399
    699400

  • Singular/fast_maps.h

    r2e2ba9 r1873199  
    88 *           bricken (Michael Brickenstein)
    99 *  Created: 01/02
    10  *  Version: $Id: fast_maps.h,v 1.13 2002-01-19 15:48:46 obachman Exp $
     10 *  Version: $Id: fast_maps.h,v 1.14 2002-01-19 18:03:38 obachman Exp $
    1111 *******************************************************************/
    1212
    … …  
    8585
    8686// evaluates mpoly and destroys it, on the fly
    87 void maPoly_Eval(mapoly mpoly, ring src_r, ideal dest_id, ring dest_r);
     87void maPoly_Eval(mapoly mpoly, ring src_r, ideal dest_id, ring dest_r, int total_cost);
    8888
    8989// creates mpoly and  mideal
Note: See TracChangeset for help on using the changeset viewer.