source: git/kernel/fast_mult.cc @ 78030c

spielwiese
Last change on this file since 78030c was 78030c, checked in by Michael Brickenstein <bricken@…>, 19 years ago
+ ring as parameter git-svn-id: file:///usr/local/Singular/svn/trunk@7727 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 5.5 KB
Line 
1#include "fast_mult.h"
2#include "kbuckets.h"
3typedef poly fastmultrec(poly f, poly g, ring r);
4static const int pass_option=1;
5static void degsplit(poly p,int n,poly &p1,poly&p2, int vn, ring r){
6  poly erg1_i, erg2_i;
7  erg1_i=NULL;
8  erg2_i=NULL;
9  while(p){
10    if(p_GetExp(p,vn,r)>=n){
11      if (p1==NULL){
12        p1=p;
13      } else{
14        pNext(erg1_i)=p;
15      }
16      erg1_i=p;
17    } else{
18      if (p2==NULL){
19        p2=p;
20      } else{
21        pNext(erg2_i)=p;
22      }
23      erg2_i=p;
24    }
25    p=pNext(p);
26  }
27  if(erg2_i){
28    pNext(erg2_i)=NULL;
29  }
30  if (erg1_i){
31    pNext(erg1_i)=NULL;
32  }
33 
34}
35static void div_by_x_power_n(poly p, int n, int vn, ring r){
36  while(p){
37    assume(p_GetExp(p,vn,r)>=n);
38    int e=p_GetExp(p,vn,r);
39    p_SetExp(p,vn,e-n,r);
40    p=pNext(p);
41  }
42}
43
44static poly do_unifastmult(poly f,int df,poly g,int dg, int vn, fastmultrec rec, ring r){
45  int n=1;
46  if ((f==NULL)||(g==NULL)) return NULL;
47  //int df=pGetExp(f,vn);//pFDeg(f);
48  //  int dg=pGetExp(g,vn);//pFDeg(g);
49       
50  int dm;
51  if(df>dg){
52    dm=df;
53  }else{
54    dm=dg;
55  }
56  while(n<=dm)
57    {
58      n*=2;
59    }
60  if(n==1){
61    return(pp_Mult_qq(f,g,r));
62  }
63 
64  int pot=n/2;
65  assume(pot*2==n);
66
67
68  //splitting
69  poly f1=NULL;
70  poly f0=NULL;//f-(x^(pot)*F1);
71  degsplit(p_Copy(f,r),pot,f1,f0,vn,r);
72  div_by_x_power_n(f1,pot,vn,r);
73
74  poly g1=NULL;
75  poly g0=NULL;
76  degsplit(p_Copy(g,r),pot,g1,g0,vn,r);
77  div_by_x_power_n(g1,pot,vn,r);
78 
79  //p00, p11
80  poly p00=rec(f0,g0,r);//unifastmult(f0,g0);
81  poly p11=rec(f1,g1,r);
82
83  //construct erg, factor
84  poly erg=NULL;
85  poly factor=pOne();
86
87  p_SetExp(factor,vn,n,r);
88  erg=pp_Mult_mm(p11,factor,r);
89  erg=p_Add_q(erg,p_Copy(p00,r),r);
90
91 
92
93
94 
95  if((f1!=NULL) &&(f0!=NULL) &&(g0!=NULL) && (g1!=NULL)){
96    //if(true){
97    //eat up f0,f1,g0,g1
98    poly s1=p_Add_q(f0,f1,r);
99    poly s2=p_Add_q(g0,g1,r);
100    poly pbig=rec(s1,s2,r);
101    p_Delete(&s1,r);
102    p_Delete(&s2,r);
103
104 
105    //eat up pbig
106    poly sum=pbig;
107    p_SetExp(factor,vn,pot,r);
108 
109    //eat up p00
110    sum=p_Add_q(sum,p_Neg(p00,r),r);
111 
112    //eat up p11
113    sum=p_Add_q(sum,p_Neg(p11,r),r);
114
115 
116    sum=p_Mult_mm(sum,factor,r);
117
118
119    //eat up sum
120    erg=p_Add_q(sum,erg,r);
121  } else {
122    //eat up f0,f1,g0,g1
123    poly s1=rec(f0,g1,r);
124    poly s2=rec(g0,f1,r);
125    p_SetExp(factor,vn,pot,r);
126    poly h=p_Mult_mm(((s1!=NULL)?s1:s2),factor,r);
127    pDelete(&f1);
128    pDelete(&f0);
129    pDelete(&g0);
130    pDelete(&g1);
131    erg=p_Add_q(erg,h,r);
132  }
133 
134 
135  p_Delete(&factor,r);
136                 
137                 
138
139  return(erg);
140}
141
142// poly do_unifastmult_buckets(poly f,poly g){
143   
144 
145
146
147//   int n=1;
148//   if ((f==NULL)||(g==NULL)) return NULL;
149//   int df=pGetExp(f,1);//pFDeg(f);
150//     int dg=pGetExp(g,1);//pFDeg(g);
151       
152//   int dm;
153//   if(df>dg){
154//     dm=df;
155//   }else{
156//     dm=dg;
157//   }
158//   while(n<=dm)
159//     {
160//       n*=2;
161//     }
162//   int pseudo_len=0;
163//   if(n==1){
164//     return ppMult_qq(f,g);
165
166//   }
167//   kBucket_pt erg_bucket= kBucketCreate(currRing);
168//   kBucketInit(erg_bucket,NULL,0 /*pLength(P.p)*/);
169
170//   int pot=n/2;
171//   assume(pot*2==n);
172
173
174//   //splitting
175//   poly f1=NULL;
176//   poly f0=NULL;//f-(x^(pot)*F1);
177//   degsplit(pCopy(f),pot,f1,f0);
178//   div_by_x_power_n(f1,pot);
179//   poly g1=NULL;
180//   poly g0=NULL;
181//   degsplit(pCopy(g),pot,g1,g0);
182//   div_by_x_power_n(g1,pot);
183 
184//   //p00
185//   //p11
186//   poly p00=unifastmult(f0,g0);
187//   poly p11=unifastmult(f1,g1);
188 
189
190
191
192//   //eat up f0,f1,g0,g1
193//   poly pbig=unifastmult(pAdd(f0,f1),pAdd(g0,g1)); 
194//   poly factor=pOne();//pCopy(erg_mult);
195//   pSetExp(factor,1,n);
196//   pseudo_len=0;
197//   kBucket_Add_q(erg_bucket,ppMult_mm(p11,factor),&pseudo_len);
198//   pseudo_len=0;
199//   kBucket_Add_q(erg_bucket,pCopy(p00),&pseudo_len);
200//   pSetExp(factor,1,pot);
201
202//   //eat up pbig
203//   pseudo_len=0;
204//   kBucket_Add_q(erg_bucket,pMult_mm(pbig,factor),&pseudo_len);
205//   pseudo_len=0;
206//   //eat up p00
207//   kBucket_Add_q(erg_bucket,pMult_mm(pNeg(p00),factor),&pseudo_len);
208//   pseudo_len=0;
209//   //eat up p11
210//   kBucket_Add_q(erg_bucket,pMult_mm(pNeg(p11),factor),&pseudo_len);
211
212 
213//   pseudo_len=0;
214
215 
216//   pDelete(&factor);
217
218//   int len=0;
219//   poly erg=NULL;
220//   kBucketClear(erg_bucket,&erg,&len);
221//   kBucketDestroy(&erg_bucket);
222
223//   return erg;
224// }
225
226static inline int max(int a, int b){
227  return (a>b)? a:b;
228}
229static inline int min(int a, int b){
230  return (a>b)? b:a;
231}
232poly unifastmult(poly f,poly g, ring r){
233  int vn=1;
234  if((f==NULL)||(g==NULL)) return NULL;
235  int df=p_GetExp(f,vn,r);
236  int dg=p_GetExp(g,vn,r);
237  if ((df==0)||(dg==0))
238    return pp_Mult_qq(f,g,r);
239  if (df*dg<100)
240    return pp_Mult_qq(f,g,r);
241  // if (df*dg>10000)
242  //  return
243  //    do_unifastmult_buckets(f,g);
244  //else
245  return do_unifastmult(f,df,g,dg,vn,unifastmult,r);
246
247}
248
249poly multifastmult(poly f, poly g, ring r){
250  if((f==NULL)||(g==NULL)) return NULL;
251  if (pLength(f)*pLength(g)<100)
252    return pp_Mult_qq(f,g,r);
253  //find vn
254  //determine df,dg simultaneously
255  int i;
256  int can_i=-1;
257  int can_df=0;
258  int can_dg=0;
259  int can_crit=0;
260  for(i=1;i<=rVar(r);i++){
261    poly p;
262    int df=0;
263    int dg=0;
264    //max min max Strategie
265    p=f;
266    while(p){
267      df=max(df,p_GetExp(p,i,r));
268      p=pNext(p);
269    }
270    if(df>can_crit){
271      p=g;
272      while(p){
273        dg=max(dg,p_GetExp(p,i,r));
274        p=pNext(p);
275      }
276      int crit=min(df,dg);
277      if (crit>can_crit){
278        can_crit=crit;
279        can_i=i;
280        can_df=df;
281        can_dg=dg;
282      }
283    }
284  }
285  if(can_crit==0)
286    return pp_Mult_qq(f,g,r);
287  else
288    {
289      poly erg=do_unifastmult(f,can_df,g,can_dg,can_i,multifastmult,r);
290      p_Normalize(erg,r);
291      return(erg);
292    }
293}
Note: See TracBrowser for help on using the repository browser.