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

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