source: git/kernel/fast_maps.cc @ fbc7cb

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