source: git/kernel/fast_maps.cc @ 788529d

spielwiese
Last change on this file since 788529d was 6aef6f2, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: fast_maps.* ->kernel git-svn-id: file:///usr/local/Singular/svn/trunk@7855 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    fast_maps.cc
6 *  Purpose: implementation of fast maps
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 02/01
9 *  Version: $Id: fast_maps.cc,v 1.2 2005-04-21 13:31:32 Singular Exp $
10 *******************************************************************/
11#include "mod2.h"
12#include <omalloc.h>
13#include "p_polys.h"
14#include "prCopy.h"
15#include "ideals.h"
16#include "ring.h"
17#include "febase.h"
18#include "fast_maps.h"
19
20// define if you want to use special dest_ring
21#define HAVE_DEST_R 1
22// define if you want to use special src_ring
23#define HAVE_SRC_R 1
24// define if you want to use optimization step
25#define HAVE_MAP_OPTIMIZE 1
26
27/*******************************************************************************
28**
29*F  maMaxExp  . . . . . . . .  returns bound on maximal exponent of result of map
30*/
31// return maximal monomial if max_map_monomials are substituted into pi_m
32static poly maGetMaxExpP(poly* max_map_monomials,
33                         int n_max_map_monomials, ring map_r,
34                         poly pi_m, ring pi_r)
35{
36  int n = si_min(pi_r->N, n_max_map_monomials);
37  int i, j;
38  Exponent_t e_i, e_j;
39  poly m_i, map_j = p_Init(map_r);
40
41  for (i=1; i <= n; i++)
42  {
43    e_i = p_GetExp(pi_m, i, pi_r);
44    m_i = max_map_monomials[i-1];
45    if (e_i > 0 && m_i != NULL && ! p_IsConstantComp(m_i, map_r))
46    {
47      for (j = 1; j<= map_r->N; j++)
48      {
49        e_j = p_GetExp(m_i, j, map_r);
50        if (e_j > 0)
51        {
52          p_AddExp(map_j, j, e_j*e_i, map_r);
53        }
54      }
55    }
56  }
57  return map_j;
58}
59
60// returns maximal exponent if map_id is applied to pi_id
61static Exponent_t maGetMaxExp(ideal pi_id, ring pi_r, ideal map_id, ring map_r)
62{
63  Exponent_t max=0;
64  poly* max_map_monomials = (poly*) omAlloc(IDELEMS(map_id)*sizeof(poly));
65  poly max_pi_i, max_map_i;
66
67  int i, j;
68  for (i=0; i<IDELEMS(map_id); i++)
69  {
70    max_map_monomials[i] = p_GetMaxExpP(map_id->m[i], map_r);
71  }
72
73  for (i=0; i<IDELEMS(pi_id); i++)
74  {
75    max_pi_i = p_GetMaxExpP(pi_id->m[i], pi_r);
76    max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r,
77                              max_pi_i, pi_r);
78    Exponent_t temp = p_GetMaxExp(max_map_i, map_r);
79    if (temp > max){
80      max=temp;
81    }
82
83    p_LmFree(max_pi_i, pi_r);
84    p_LmFree(max_map_i, map_r);
85  }
86  return max;
87}
88
89
90/*******************************************************************************
91**
92*F  debugging stuff
93*/
94#ifndef NDEBUG
95void maMonomial_Out(mapoly monomial, ring src_r, ring dest_r)
96{
97  p_wrp(monomial->src, src_r);
98  printf(" ref:%d", monomial->ref);
99  if (dest_r != NULL)
100  {
101    printf(" dest:");
102    p_wrp(monomial->dest, dest_r);
103  }
104  if (monomial->f1!=NULL) { printf(" f1:%x", (unsigned int)monomial->f1->src);
105                            // p_wrp(monomial->f1->src, src_r);
106                          }
107  if (monomial->f2!=NULL) { printf(" f2:%x",(unsigned int)monomial->f2->src);
108                            // p_wrp(monomial->f2->src, src_r);
109                          }
110  printf("\n");
111  fflush(stdout);
112}
113
114void maPoly_Out(mapoly mpoly, ring src_r, ring dest_r)
115{
116  while (mpoly != NULL)
117  {
118    maMonomial_Out(mpoly, src_r, dest_r);
119    mpoly = mpoly->next;
120  }
121}
122#endif
123
124/*******************************************************************************
125**
126*F  mapolyCreate  . . . . . . . . . . . . . . .  Creates mapoly
127*/
128static omBin mapolyBin = omGetSpecBin(sizeof(mapoly_s));
129static omBin macoeffBin = omGetSpecBin(sizeof(macoeff_s));
130
131mapoly maMonomial_Create(poly p, ring r_p, sBucket_pt bucket)
132{
133  mapoly mp = (mapoly) omAlloc0Bin(mapolyBin);
134  //p_wrp(p,r_p);printf(" (%x) created\n",mp);
135  mp->src = p;
136  p->next = NULL;
137
138  if (bucket != NULL)
139  {
140    mp->coeff = (macoeff) omAlloc0Bin(macoeffBin);
141    mp->coeff->bucket = bucket;
142    mp->coeff->n = pGetCoeff(p);
143  }
144  mp->ref = 1;
145  return mp;
146}
147
148void maMonomial_Destroy(mapoly mp, ring src_r, ring dest_r)
149{
150  if (mp != NULL)
151  {
152    p_LmFree(mp->src, src_r);
153    if (mp->coeff != NULL)
154    {
155      macoeff coeff, next = mp->coeff;
156      do
157      {
158        coeff = next;
159        next = coeff->next;
160        omFreeBin(coeff, macoeffBin);
161      }
162      while (next != NULL);
163      if (mp->dest != NULL)
164      {
165        assume(dest_r != NULL);
166        p_Delete(&(mp->dest), dest_r);
167      }
168    }
169  }
170  omFreeBin(mp, mapolyBin);
171}
172
173/*******************************************************************************
174**
175*F  maPoly_InsertMonomial . . . . . . . . .insertion of a monomial into mapoly
176*/
177mapoly maPoly_InsertMonomial(mapoly &into, mapoly what, ring src_r)
178{
179  if (into == NULL)
180  {
181    into = what;
182    return what;
183  }
184
185  mapoly iter = into;
186  mapoly prev = NULL;
187
188  Top:
189  p_LmCmpAction(iter->src, what->src, src_r, goto Equal, goto Greater, goto Smaller);
190
191  Greater:
192  if (iter->next == NULL)
193  {
194    iter->next = what;
195    return what;
196  }
197  prev = iter;
198  iter = iter->next;
199  goto Top;
200
201  Smaller:
202  if (prev == NULL)
203  {
204    into = what;
205    what->next = iter;
206    return what;
207  }
208  prev->next = what;
209  what->next = iter;
210  return what;
211
212  Equal:
213  iter->ref += what->ref;
214  macoeff coeff = what->coeff;
215  if (coeff != NULL)
216  {
217    while (coeff->next != NULL) coeff = coeff->next;
218    coeff->next = iter->coeff;
219    iter->coeff = what->coeff;
220    what->coeff = NULL;
221  }
222  maMonomial_Free(what, src_r);
223  return iter;
224}
225
226mapoly maPoly_InsertMonomial(mapoly &into, poly p, ring src_r, sBucket_pt bucket)
227{
228  return maPoly_InsertMonomial(into, maMonomial_Create(p, src_r, bucket), src_r);
229}
230
231static void maPoly_InsertPoly(mapoly &into, poly what, ring src_r, sBucket_pt bucket)
232{
233  poly next;
234
235  while (what != NULL)
236  {
237    next = what->next;
238    maPoly_InsertMonomial(into, what, src_r, bucket);
239    what = next;
240  }
241}
242
243/*******************************************************************************
244**
245*F  maMap_Create . . . . . . . . . . . . . . . . . . . . create stuff
246*/
247
248void maMap_CreatePolyIdeal(ideal map_id, ring map_r, ring src_r, ring dest_r,
249                           mapoly &mp, maideal &mideal)
250{
251  mideal = (maideal) omAlloc0(sizeof(maideal_s));
252  mideal->n = IDELEMS(map_id);
253  mideal->buckets = (sBucket_pt*) omAlloc0(mideal->n*sizeof(sBucket_pt));
254  int i;
255  mp = NULL;
256
257  for (i=0; i<mideal->n; i++)
258  {
259    if (map_id->m[i] != NULL)
260    {
261      mideal->buckets[i] = sBucketCreate(dest_r);
262      maPoly_InsertPoly(mp,
263#ifdef PDEBUG
264                        prShallowCopyR(map_id->m[i], map_r, src_r),
265#else
266                        prShallowCopyR_NoSort(map_id->m[i], map_r, src_r),
267#endif
268                        src_r,
269                        mideal->buckets[i]);
270    }
271  }
272}
273
274
275void maMap_CreateRings(ideal map_id, ring map_r,
276                       ideal image_id, ring image_r,
277                       ring &src_r, ring &dest_r, BOOLEAN &simple)
278{
279#if HAVE_SRC_R > 0
280  int* weights = (int*) omAlloc0(map_r->N*sizeof(int));
281  int i;
282  int n = si_min(map_r->N, IDELEMS(image_id));
283
284  for (i=0; i<n; i++)
285  {
286    weights[i] = pLength(image_id->m[i]);
287  }
288  src_r = rModifyRing_Wp(map_r, weights);
289#else
290  src_r = map_r;
291#endif
292
293#if HAVE_DEST_R > 0
294  Exponent_t maxExp = maGetMaxExp(map_id, map_r, image_id, image_r);
295  if (maxExp <=  1) maxExp = 2;
296  else if (maxExp > (Exponent_t) image_r->bitmask)
297    maxExp = (Exponent_t) image_r->bitmask;
298  dest_r = rModifyRing_Simple(image_r, TRUE, TRUE, maxExp,  simple);
299#else
300  dest_r = image_r;
301#endif
302}
303
304static void maMap_KillRings(ring map_r, ring image_r, ring src_r, ring dest_r)
305{
306  if (map_r != src_r)
307    rKillModified_Wp_Ring(src_r);
308  if (image_r != dest_r)
309    rKillModifiedRing_Simple(dest_r);
310}
311
312/*******************************************************************************
313**
314*F  misc  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  misc  stuff
315*/
316
317ideal maIdeal_2_Ideal(maideal m_id, ring dest_r)
318{
319  ideal res = idInit(m_id->n, 1);
320  int l;
321
322  for (int i= 0; i < m_id->n; i++)
323  {
324    if (m_id->buckets[i]!=NULL)
325      sBucketDestroyAdd(m_id->buckets[i], &(res->m[i]), &l);
326  }
327  omFree(m_id);
328  return res;
329}
330
331void maPoly_GetLength(mapoly mp, int &length)
332{
333  length = 0;
334  while (mp != NULL)
335  {
336    length++;
337    mp = mp->next;
338  }
339}
340
341
342/*******************************************************************************
343**
344*F  fast_map  . . . . . . . . . . . . . . . . . . . . . . . . . .the real thing
345*/
346
347ideal fast_map(ideal map_id, ring map_r, ideal image_id, ring image_r)
348{
349  ring src_r, dest_r;
350  ideal dest_id, res_id;
351  int length = 0;
352  BOOLEAN no_sort;
353
354  // construct rings we work in:
355  // src_r: Wp with Weights set to length of poly in image_id
356  // dest_r: Simple ring without degree ordering and short exponents
357  maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r, no_sort);
358
359  // construct dest_id
360  if (dest_r != image_r)
361    dest_id = idrShallowCopyR(image_id, image_r, dest_r);
362  else
363    dest_id = image_id;
364
365  // construct mpoly and mideal
366  mapoly mp;
367  maideal mideal;
368  maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);
369
370  if (TEST_OPT_PROT)
371  {
372    maPoly_GetLength(mp, length);
373    Print("map[%d:%d]{%d:", dest_r->bitmask, dest_r->ExpL_Size, length);
374  }
375
376  // do the optimization step
377#if HAVE_MAP_OPTIMIZE > 0
378  if (mp!=NULL) maPoly_Optimize(mp, src_r);
379#endif
380  if (TEST_OPT_PROT)
381  {
382    maPoly_GetLength(mp, length);
383    Print("%d}", length);
384  }
385
386  // do the actual evaluation
387  maPoly_Eval(mp, src_r, dest_id, dest_r, length);
388  if (TEST_OPT_PROT) Print(".");
389
390  // collect the results back into an ideal
391  ideal res_dest_id = maIdeal_2_Ideal(mideal, dest_r);
392  if (TEST_OPT_PROT) Print(".");
393
394  // convert result back to image_r
395  ideal res_image_id;
396  if (dest_r != image_r)
397  {
398    //if (no_sort) see Old/m134si.tst
399    //  res_image_id = idrShallowCopyR_NoSort(res_dest_id, dest_r, image_r);
400    //else
401      res_image_id = idrShallowCopyR(res_dest_id, dest_r, image_r);
402    id_ShallowDelete(&res_dest_id, dest_r);
403  }
404  else
405    res_image_id = res_dest_id;
406
407  if (TEST_OPT_PROT) Print(".");
408  // clean-up the rings
409  maMap_KillRings(map_r, image_r, src_r, dest_r);
410
411  if (TEST_OPT_PROT)
412    Print("\n");
413
414  idTest(res_image_id);
415  return res_image_id;
416}
417
418
419/**********************************************************************
420* Evaluation stuff                                                    *
421**********************************************************************/
422
423// substitute p everywhere the monomial occours,
424// return the number of substitutions
425static int maPoly_Substitute(macoeff c, poly p, ring dest_r)
426{
427  // substitute the monomial: go through macoeff
428  int len=pLength(p);
429  int done=0;
430  while (c!=NULL)
431  {
432    done++;
433    poly t=pp_Mult_nn(p,c->n,dest_r);
434    sBucket_Add_p(c->bucket, t, len);
435    c=c->next;
436  }
437  return done;
438}
439
440static poly maPoly_EvalMon(poly src, ring src_r, poly* dest_id, ring dest_r)
441{
442  int i;
443  int e;
444  poly p=NULL;
445  poly pp;
446  for(i=1;i<=src_r->N;i++)
447  {
448    e=p_GetExp(src,i,src_r);
449    if (e>0)
450    {
451      pp=dest_id[i-1];
452      if (pp==NULL)
453      {
454        p_Delete(&p,dest_r);
455        return NULL;
456      }
457      if ((p==NULL) /* && (e>0)*/)
458      {
459        p=p_Copy(pp /*dest_id[i-1]*/,dest_r);
460        e--;
461      }
462      while (e>0)
463      {
464        p=p_Mult_q(p,p_Copy(pp /*dest_id[i-1]*/,dest_r),dest_r);
465        e--;
466      }
467    }
468  }
469  if (p==NULL) p=p_ISet(1,dest_r);
470  return p;
471}
472
473void maPoly_Eval(mapoly root, ring src_r, ideal dest_id, ring dest_r, int total_cost)
474{
475  // invert the list rooted at root:
476  if ((root!=NULL) && (root->next!=NULL))
477  {
478    mapoly q=root->next;
479    mapoly qn;
480    root->next=NULL;
481    do
482    {
483      qn=q->next;
484      q->next=root;
485      root=q;
486      q=qn;
487    }
488    while (qn !=NULL);
489  }
490
491  total_cost /= 10;
492  int next_print_cost = total_cost;
493
494  // the evaluation -----------------------------------------
495  mapoly p=root;
496  int cost = 0;
497
498  while (p!=NULL)
499  {
500    // look at each mapoly: compute its value in ->dest
501    assume (p->dest==NULL);
502    {
503      if ((p->f1!=NULL)&&(p->f2!=NULL))
504      {
505#if 0
506        printf("found prod:");
507        p_wrp(p->src,src_r);printf("=");
508        p_wrp(p->f1->src,src_r);printf(" * ");
509        p_wrp(p->f2->src,src_r);printf("\n");
510#endif
511        poly f1=p->f1->dest;
512        poly f2=p->f2->dest;
513        if (p->f1->ref>0) f1=p_Copy(f1,dest_r);
514        else
515        {
516          // we own p->f1->dest now (in f1)
517          p->f1->dest=NULL;
518        }
519        if (p->f2->ref>0) f2=p_Copy(f2,dest_r);
520        else
521        {
522          // we own p->f2->dest now (in f2)
523          p->f2->dest=NULL;
524        }
525        maMonomial_Free(p->f1,src_r, dest_r);
526        maMonomial_Free(p->f2,src_r, dest_r);
527        p->dest=p_Mult_q(f1,f2,dest_r);
528      } /* factors : 2 */
529      else
530      {
531        //printf("compute "); p_wrp(p->src,src_r);printf("\n");
532        assume((p->f1==NULL) && (p->f2==NULL));
533        //if(p->f1!=NULL) { printf(" but f1="); p_wrp(p->f1->src,src_r);printf("\n"); }
534        //if(p->f2!=NULL) { printf(" but f2="); p_wrp(p->f2->src,src_r);printf("\n"); }
535        // no factorization provided, use the classical method:
536        p->dest=maPoly_EvalMon(p->src,src_r,dest_id->m,dest_r);
537        //p_wrp(p->dest, dest_r); printf(" is p->dest\n");
538      }
539    } /* p->dest==NULL */
540    // substitute the monomial: go through macoeff
541    p->ref -= maPoly_Substitute(p->coeff, p->dest, dest_r);
542    //printf("subst done\n");
543    if (total_cost)
544    {
545      assume(TEST_OPT_PROT);
546      cost++;
547      if (cost > next_print_cost)
548      {
549        Print("-");
550        next_print_cost += total_cost;
551      }
552    }
553
554    mapoly pp=p;
555    p=p->next;
556    //p_wrp(pp->src, src_r);
557    if (pp->ref<=0)
558    {
559      //printf(" (%x) killed\n",pp);
560      maMonomial_Destroy(pp, src_r, dest_r);
561    }
562    //else
563    //  printf(" (%x) not killed, ref=%d\n",pp,pp->ref);
564  }
565}
566
567
568/*******************************************************************************
569**
570*F  maEggt  . . . . . . . . . . . . . . . . . . . . . . . .  returns extended ggt
571*/
572// return NULL if deg(ggt(m1, m2)) < 2
573// else return m = ggT(m1, m2) and q1, q2 such that m1 = q1*m m2 = q2*m
574static poly maEggT(const poly m1, const poly m2, poly &q1, poly &q2,const ring r)
575{
576
577  int i;
578  int dg = 0;
579  poly ggt = NULL;
580  q1 = p_Init(r);
581  q2 = p_Init(r);
582  ggt=p_Init(r);
583
584  for (i=1;i<=r->N;i++) {
585    Exponent_t e1 = p_GetExp(m1, i, r);
586    Exponent_t e2 = p_GetExp(m2, i, r);
587    if (e1 > 0 && e2 > 0){
588      Exponent_t em = (e1 > e2 ? e2 : e1);
589      dg += em;
590      p_SetExp(ggt, i, em, r);
591      p_SetExp(q1, i, e1 - em, r);
592      p_SetExp(q2, i, e2 - em, r);
593    }
594    else {
595      p_SetExp(q1, i, e1, r);
596      p_SetExp(q2, i, e2, r);
597    }
598  }
599  if (dg>1)
600  {
601    p_Setm(ggt, r);
602    p_Setm(q1, r);
603    p_Setm(q2, r);
604
605
606  }
607  else {
608    p_LmFree(ggt, r);
609    p_LmFree(q1, r);
610    p_LmFree(q2, r);
611    ggt = NULL;
612  }
613  return ggt;
614}
615
616/********************************************************************
617 **                                                                 *
618 * maFindBestggT                                                    *
619 * finds ggT with the highest cost                                  *
620 *******************************************************************/
621
622static mapoly maFindBestggT(mapoly mp, mapoly & choice, mapoly & fp, mapoly & fq,const ring r)
623{
624  int ggt_deg = 0;
625  poly p = mp->src;
626  mapoly iter = choice;
627  poly ggT = NULL;
628  fp = NULL;
629  fq = NULL;
630  poly fp_p=NULL;
631  poly fq_p=NULL;
632  choice=NULL;
633  while ((iter != NULL) && (pDeg(iter->src, r) > ggt_deg))
634  {
635    //    maMonomial_Out(iter, r, NULL);
636    poly q1, q2, q;
637
638    q = maEggT(p, iter->src, q1, q2,r);
639    if (q != NULL)
640    {
641      assume((q1!=NULL)&&(q2!=NULL));
642      if (pDeg(q,r) > ggt_deg)
643      {
644        choice=iter;
645        if (ggT != NULL)
646        {
647          p_LmFree(ggT, r);
648          p_LmFree(fp_p, r);
649          p_LmFree(fq_p, r);
650
651        }
652        ggt_deg = pDeg(q, r);
653        ggT = q;
654        fp_p = q1;
655        fq_p = q2;
656      }
657      else
658      {
659        p_LmFree(q, r);
660        p_LmFree(q1, r);
661        p_LmFree(q2, r);
662      }
663    }
664    iter=iter->next;
665  }
666  if(ggT!=NULL){
667
668    int dq =pTotaldegree(fq_p,r);
669    if (dq!=0){
670      fq=maPoly_InsertMonomial(mp, fq_p, r, NULL);
671      fp=maPoly_InsertMonomial(mp, fp_p, r, NULL);
672      return maPoly_InsertMonomial(mp, ggT, r, NULL);
673    }
674    else {
675      fq=NULL;
676      p_LmFree(fq_p, r);
677      p_LmFree(ggT, r);
678      fp=maPoly_InsertMonomial(mp, fp_p, r, NULL);
679      choice->ref++;
680      return choice;
681    }
682  }
683  else {
684    return NULL;
685  }
686
687
688}
689
690/********************************************************************
691 **                                                                 *
692 * maPoly_Optimize                                                  *
693 * adds and integrates subexpressions                               *
694 *******************************************************************/
695
696void maPoly_Optimize(mapoly mpoly, ring src_r){
697  assume(mpoly!=NULL && mpoly->src!=NULL);
698  mapoly iter = mpoly;
699  mapoly choice;
700  mapoly ggT=NULL;
701  mapoly fp=NULL;
702  mapoly fq=NULL;
703  while (iter->next!=NULL){
704    choice=iter->next;
705    if ((iter->f1==NULL)){
706      ggT=maFindBestggT(iter, choice, fp, fq,src_r);
707      if (choice!=NULL){
708        iter->f1=fp;
709        iter->f2=ggT;
710        if (fq!=NULL){
711          ggT->ref++;
712          choice->f1=fq;
713          choice->f2=ggT;
714        }
715      }
716
717    }
718    iter=iter->next;
719  }
720}
Note: See TracBrowser for help on using the repository browser.