source: git/kernel/fast_maps.cc @ abe5c8

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