# source:git/kernel/maps/fast_maps.cc@213596e

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