source: git/kernel/fast_mult.cc @ 6b9532

spielwiese
Last change on this file since 6b9532 was ec8a6b6, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: om(re)alloc -> om(Re)Alloc MACROS git-svn-id: file:///usr/local/Singular/svn/trunk@8124 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 12.6 KB
Line 
1#include "fast_mult.h"
2#include "kbuckets.h"
3#include "febase.h"
4typedef poly fastmultrec(poly f, poly g, ring r);
5static const int pass_option=1;
6static int mults=0;
7int Mults(){
8  return mults;
9}
10static void degsplit(poly p,int n,poly &p1,poly&p2, int vn, ring r){
11  poly erg1_i, erg2_i;
12  erg1_i=NULL;
13  erg2_i=NULL;
14  while(p){
15    if(p_GetExp(p,vn,r)>=n){
16      if (p1==NULL){
17        p1=p;
18      } else{
19        pNext(erg1_i)=p;
20      }
21      erg1_i=p;
22    } else{
23      if (p2==NULL){
24        p2=p;
25      } else{
26        pNext(erg2_i)=p;
27      }
28      erg2_i=p;
29    }
30    p=pNext(p);
31  }
32  if(erg2_i){
33    pNext(erg2_i)=NULL;
34  }
35  if (erg1_i){
36    pNext(erg1_i)=NULL;
37  }
38 
39}
40static void div_by_x_power_n(poly p, int n, int vn, ring r){
41  while(p){
42    assume(p_GetExp(p,vn,r)>=n);
43    int e=p_GetExp(p,vn,r);
44    p_SetExp(p,vn,e-n,r);
45    p=pNext(p);
46  }
47}
48
49static poly do_unifastmult(poly f,int df,poly g,int dg, int vn, fastmultrec rec, ring r){
50  int n=1;
51  if ((f==NULL)||(g==NULL)) return NULL;
52  //int df=pGetExp(f,vn);//pFDeg(f);
53  //  int dg=pGetExp(g,vn);//pFDeg(g);
54       
55  int dm;
56  if(df>dg){
57    dm=df;
58  }else{
59    dm=dg;
60  }
61  while(n<=dm)
62    {
63      n*=2;
64    }
65  if(n==1){
66    return(pp_Mult_qq(f,g,r));
67  }
68 
69  int pot=n/2;
70  assume(pot*2==n);
71
72
73  //splitting
74  poly f1=NULL;
75  poly f0=NULL;//f-(x^(pot)*F1);
76  degsplit(p_Copy(f,r),pot,f1,f0,vn,r);
77  div_by_x_power_n(f1,pot,vn,r);
78
79  poly g1=NULL;
80  poly g0=NULL;
81  degsplit(p_Copy(g,r),pot,g1,g0,vn,r);
82  div_by_x_power_n(g1,pot,vn,r);
83 
84  //p00, p11
85  poly p00=rec(f0,g0,r);//unifastmult(f0,g0);
86  poly p11=rec(f1,g1,r);
87
88  //construct erg, factor
89  poly erg=NULL;
90  poly factor=p_ISet(1,r);
91
92  p_SetExp(factor,vn,n,r);
93  erg=pp_Mult_mm(p11,factor,r);
94  erg=p_Add_q(erg,p_Copy(p00,r),r);
95
96 
97
98
99 
100  if((f1!=NULL) &&(f0!=NULL) &&(g0!=NULL) && (g1!=NULL)){
101    //if(true){
102    //eat up f0,f1,g0,g1
103    poly s1=p_Add_q(f0,f1,r);
104    poly s2=p_Add_q(g0,g1,r);
105    poly pbig=rec(s1,s2,r);
106    p_Delete(&s1,r);
107    p_Delete(&s2,r);
108
109 
110    //eat up pbig
111    poly sum=pbig;
112    p_SetExp(factor,vn,pot,r);
113 
114    //eat up p00
115    sum=p_Add_q(sum,p_Neg(p00,r),r);
116 
117    //eat up p11
118    sum=p_Add_q(sum,p_Neg(p11,r),r);
119
120 
121    sum=p_Mult_mm(sum,factor,r);
122
123
124    //eat up sum
125    erg=p_Add_q(sum,erg,r);
126  } else {
127    //eat up f0,f1,g0,g1
128    poly s1=rec(f0,g1,r);
129    poly s2=rec(g0,f1,r);
130    p_SetExp(factor,vn,pot,r);
131    poly h=p_Mult_mm(((s1!=NULL)?s1:s2),factor,r);
132    p_Delete(&f1,r);
133    p_Delete(&f0,r);
134    p_Delete(&g0,r);
135    p_Delete(&g1,r);
136    p_Delete(&p00,r);
137    p_Delete(&p11,r);
138    erg=p_Add_q(erg,h,r);
139  }
140 
141 
142  p_Delete(&factor,r);
143                 
144                 
145
146  return(erg);
147}
148
149// poly do_unifastmult_buckets(poly f,poly g){
150   
151 
152
153
154//   int n=1;
155//   if ((f==NULL)||(g==NULL)) return NULL;
156//   int df=pGetExp(f,1);//pFDeg(f);
157//     int dg=pGetExp(g,1);//pFDeg(g);
158       
159//   int dm;
160//   if(df>dg){
161//     dm=df;
162//   }else{
163//     dm=dg;
164//   }
165//   while(n<=dm)
166//     {
167//       n*=2;
168//     }
169//   int pseudo_len=0;
170//   if(n==1){
171//     return ppMult_qq(f,g);
172
173//   }
174//   kBucket_pt erg_bucket= kBucketCreate(currRing);
175//   kBucketInit(erg_bucket,NULL,0 /*pLength(P.p)*/);
176
177//   int pot=n/2;
178//   assume(pot*2==n);
179
180
181//   //splitting
182//   poly f1=NULL;
183//   poly f0=NULL;//f-(x^(pot)*F1);
184//   degsplit(pCopy(f),pot,f1,f0);
185//   div_by_x_power_n(f1,pot);
186//   poly g1=NULL;
187//   poly g0=NULL;
188//   degsplit(pCopy(g),pot,g1,g0);
189//   div_by_x_power_n(g1,pot);
190 
191//   //p00
192//   //p11
193//   poly p00=unifastmult(f0,g0);
194//   poly p11=unifastmult(f1,g1);
195 
196
197
198
199//   //eat up f0,f1,g0,g1
200//   poly pbig=unifastmult(pAdd(f0,f1),pAdd(g0,g1)); 
201//   poly factor=pOne();//pCopy(erg_mult);
202//   pSetExp(factor,1,n);
203//   pseudo_len=0;
204//   kBucket_Add_q(erg_bucket,ppMult_mm(p11,factor),&pseudo_len);
205//   pseudo_len=0;
206//   kBucket_Add_q(erg_bucket,pCopy(p00),&pseudo_len);
207//   pSetExp(factor,1,pot);
208
209//   //eat up pbig
210//   pseudo_len=0;
211//   kBucket_Add_q(erg_bucket,pMult_mm(pbig,factor),&pseudo_len);
212//   pseudo_len=0;
213//   //eat up p00
214//   kBucket_Add_q(erg_bucket,pMult_mm(pNeg(p00),factor),&pseudo_len);
215//   pseudo_len=0;
216//   //eat up p11
217//   kBucket_Add_q(erg_bucket,pMult_mm(pNeg(p11),factor),&pseudo_len);
218
219 
220//   pseudo_len=0;
221
222 
223//   pDelete(&factor);
224
225//   int len=0;
226//   poly erg=NULL;
227//   kBucketClear(erg_bucket,&erg,&len);
228//   kBucketDestroy(&erg_bucket);
229
230//   return erg;
231// }
232
233static inline int max(int a, int b){
234  return (a>b)? a:b;
235}
236static inline int min(int a, int b){
237  return (a>b)? b:a;
238}
239poly unifastmult(poly f,poly g, ring r){
240  int vn=1;
241  if((f==NULL)||(g==NULL)) return NULL;
242  int df=p_GetExp(f,vn,r);
243  int dg=p_GetExp(g,vn,r);
244  if ((df==0)||(dg==0))
245    return pp_Mult_qq(f,g,r);
246  if (df*dg<100)
247    return pp_Mult_qq(f,g,r);
248  // if (df*dg>10000)
249  //  return
250  //    do_unifastmult_buckets(f,g);
251  //else
252  return do_unifastmult(f,df,g,dg,vn,unifastmult,r);
253
254}
255
256poly multifastmult(poly f, poly g, ring r){
257  mults++;
258  if((f==NULL)||(g==NULL)) return NULL;
259  if (pLength(f)*pLength(g)<100)
260    return pp_Mult_qq(f,g,r);
261  //find vn
262  //determine df,dg simultaneously
263  int i;
264  int can_i=-1;
265  int can_df=0;
266  int can_dg=0;
267  int can_crit=0;
268  for(i=1;i<=rVar(r);i++){
269    poly p;
270    int df=0;
271    int dg=0;
272    //max min max Strategie
273    p=f;
274    while(p){
275      df=max(df,p_GetExp(p,i,r));
276      p=pNext(p);
277    }
278    if(df>can_crit){
279      p=g;
280      while(p){
281        dg=max(dg,p_GetExp(p,i,r));
282        p=pNext(p);
283      }
284      int crit=min(df,dg);
285      if (crit>can_crit){
286        can_crit=crit;
287        can_i=i;
288        can_df=df;
289        can_dg=dg;
290      }
291    }
292  }
293  if(can_crit==0)
294    return pp_Mult_qq(f,g,r);
295  else
296    {
297      poly erg=do_unifastmult(f,can_df,g,can_dg,can_i,multifastmult,r);
298      p_Normalize(erg,r);
299      return(erg);
300    }
301}
302poly pFastPower(poly f, int n, ring r){
303  if (n==1) return f;
304  if (n==0) return p_ISet(1,r);
305  assume(n>=0);
306  int i_max=1;
307  int pot_max=0;
308  while(i_max*2<=n){
309    i_max*=2;
310    pot_max++;
311  }
312  int field_size=pot_max+1;
313  int* int_pot_array=(int*) omalloc(field_size*sizeof(int));
314  poly* pot_array=(poly*) omalloc(field_size*sizeof(poly));
315  int i;
316  int pot=1;
317  //initializing int_pot
318  for(i=0;i<field_size;i++){
319    int_pot_array[i]=pot;
320    pot*=2;
321  }
322  //calculating pot_array
323  pot_array[0]=f; //do not delete it
324  for(i=1;i<field_size;i++){
325    poly p=pot_array[i-1];
326    if(rVar(r)==1)
327      pot_array[i]=multifastmult(p,p,r);
328    else
329      pot_array[i]=pp_Mult_qq(p,p,r);
330  }
331 
332
333
334  int work_n=n;
335  assume(work_n>=int_pot_array[field_size-1]);
336  poly erg=p_ISet(1,r);
337
338
339  //forward maybe faster, later
340  // for(i=field_size-1;i>=0;i--){
341
342//       assume(work_n<2*int_pot_array[i]);
343//       if(int_pot_array[i]<=work_n){
344//      work_n-=int_pot_array[i];
345//      poly prod=multifastmult(erg,pot_array[i],r);
346//      pDelete(&erg);
347//      erg=prod;
348//       }
349
350//       if(i!=0) pDelete(&pot_array[i]);
351//   }
352 
353 
354  for(i=field_size-1;i>=0;i--){
355
356      assume(work_n<2*int_pot_array[i]);
357      if(int_pot_array[i]<=work_n){
358        work_n-=int_pot_array[i];
359        int_pot_array[i]=1;
360      }
361      else int_pot_array[i]=0;
362   
363  }
364  for(i=0;i<field_size;i++){
365    if(int_pot_array[i]==1){
366      poly prod;
367      if(rVar(r)==1)
368        prod=multifastmult(erg,pot_array[i],r);
369      else
370        prod=pp_Mult_qq(erg,pot_array[i],r);
371      pDelete(&erg);
372      erg=prod;
373    }
374    if(i!=0) pDelete(&pot_array[i]);
375  }
376  //free res
377  omfree(pot_array);
378  omfree(int_pot_array);
379  return erg;
380}
381static omBin lm_bin=NULL;
382/*3
383* compute for monomials p*q
384* destroys p, keeps q
385*/
386static void p_MonMultMB(poly p, poly q,ring r)
387{
388  number x, y;
389  // int i;
390
391  y = p_GetCoeff(p,r);
392  x = n_Mult(y,pGetCoeff(q),r);
393  n_Delete(&y,r);
394  p_SetCoeff0(p,x,r);
395  //for (i=pVariables; i!=0; i--)
396  //{
397  //  pAddExp(p,i, pGetExp(q,i));
398  //}
399  //p->Order += q->Order;
400  p_ExpVectorAdd(p,q,r);
401}
402/*3
403* compute for monomials p*q
404* keeps p, q
405* uses bin only available in MCPower
406*/
407static poly p_MonMultCMB(poly p, poly q, ring r)
408{
409  number x;
410  int i;
411  poly res = p_Init(r,lm_bin);
412
413  x = n_Mult(p_GetCoeff(p,r),p_GetCoeff(q,r),r);
414  p_SetCoeff0(res,x,r);
415  p_ExpVectorSum(res,p, q,r);
416  return res;
417}
418static poly p_MonPowerMB(poly p, int exp, ring r)
419{
420  int i;
421
422  if(!n_IsOne(p_GetCoeff(p,r),r))
423  {
424    number x, y;
425    y = p_GetCoeff(p,r);
426    n_Power(y,exp,&x,r);
427    n_Delete(&y,r);
428    p_SetCoeff0(p,x,r);
429  }
430  for (i=rVar(r); i!=0; i--)
431  {
432    p_MultExp(p,i, exp,r);
433  }
434  p_Setm(p,r);
435  return p;
436}
437static void buildTermAndAdd(int n,number* facult,poly* f_terms,int* exp,int f_len,kBucket_pt erg_bucket,ring r, number coef, poly & zw, poly tmp, poly** term_pot){
438 
439  int i;
440  //  poly term=p_Init(r);
441
442   //  number denom=n_Init(1,r);
443//   for(i=0;i<f_len;i++){
444//     if(exp[i]!=0){
445//       number trash=denom;
446//       denom=n_Mult(denom,facult[exp[i]],r);
447//       n_Delete(&trash,r);
448//     }
449   
450//   }
451//   number coef=n_IntDiv(facult[n],denom,r);   //right function here?
452//   n_Delete(&denom,r);
453//   poly erg=p_NSet(coef,r);
454  poly erg=p_Init(r,lm_bin);
455  p_SetCoeff0(erg, coef,r);
456  //p_NSet(n_Copy(coef,r),r);
457  for(i=0;i<f_len;i++){
458    if(exp[i]!=0){
459      ///poly term=p_Copy(f_terms[i],r);
460      poly term=term_pot[i][exp[i]];
461        //tmp;
462        //p_ExpVectorCopy(term,f_terms[i],r);
463        //p_SetCoeff(term, n_Copy(p_GetCoeff(f_terms[i],r),r),r);
464       
465      //term->next=NULL;
466     
467      //p_MonPowerMB(term, exp[i],r);
468      p_MonMultMB(erg,term,r);
469      //p_Delete(&term,r);
470    }
471   
472  }
473  int pseudo_len=1;
474  zw=erg;
475  //  kBucket_Add_q(erg_bucket,erg,&pseudo_len);
476}
477
478
479
480static void MC_iterate(poly f, int n, ring r, int f_len,number* facult, int* exp,poly* f_terms,kBucket_pt erg_bucket,int pos,int sum, number coef, poly & zw, poly tmp, poly** term_pot){
481  int i;
482 
483  if (pos<f_len-1){
484    poly zw_l=NULL;
485    number new_coef;
486    for(i=0;i<=n-sum;i++){
487      exp[pos]=i;
488      if(i==0){
489        new_coef=n_Copy(coef,r);
490      }
491      else {
492        number old=new_coef;
493        number old_rest=n_Init(n-sum-(i-1),r);
494        new_coef=n_Mult(new_coef,old_rest,r);
495        n_Delete(&old_rest,r);
496        n_Delete(&old,r);
497        number i_number=n_Init(i,r);
498        old=new_coef;
499        new_coef=n_IntDiv(new_coef,i_number,r);
500        n_Delete(&old,r);
501        n_Delete(&i_number,r);
502      }
503      //new_coef is
504      //(n                          )
505      //(exp[0]..exp[pos] 0 0 0 rest)
506      poly zw_real=NULL;
507      MC_iterate(f, n, r, f_len,facult, exp,f_terms,erg_bucket,pos+1,sum+i,new_coef,zw_real,tmp,term_pot);
508      if (pos==f_len-2){
509        //get first small polys
510       
511        zw_real->next=zw_l;
512        zw_l=zw_real;
513      }
514      //n_Delete(& new_coef,r);
515    }
516    n_Delete(&new_coef,r);
517    if (pos==f_len-2){
518      int len=n-sum+1;
519      kBucket_Add_q(erg_bucket,zw_l,&len);
520    }
521    return;
522  }
523  if(pos==f_len-1){
524    i=n-sum;
525    exp[pos]=i;
526    number new_coef=nCopy(coef);//n_IntDiv(coef,facult[i],r); //really consumed???????
527    buildTermAndAdd(n,facult,f_terms,exp,f_len,erg_bucket,r, new_coef,zw, tmp,term_pot);
528    // n_Delete(& new_coef,r);
529  }
530  assume(pos<=f_len-1);
531}
532poly pFastPowerMC(poly f, int n, ring r){
533  //only char=0
534 
535  if(rChar(r)!=0)
536    Werror("Char not 0, pFastPowerMC not implemented for this case");
537  if (n<=1)
538    Werror("not implemented for so small n, recursion fails");//should be length(f)
539   if (pLength(f)<=1)
540    Werror("not implemented for so small lenght of f, recursion fails");
541  //  number null_number=n_Init(0,r);
542  number* facult=(number*) omalloc((n+1)*sizeof(number));
543  facult[0]=n_Init(1,r);
544  int i;
545  for(i=1;i<=n;i++){
546    number this_n=n_Init(i,r);
547    facult[i]=n_Mult(this_n,facult[i-1],r);
548    n_Delete(&this_n,r);
549  }
550
551  lm_bin=omGetSpecBin(POLYSIZE + (r->ExpL_Size)*sizeof(long));
552
553  kBucket_pt erg_bucket= kBucketCreate(currRing);
554  kBucketInit(erg_bucket,NULL,0);
555  const int f_len=pLength(f);
556  int* exp=(int*)omalloc(f_len*sizeof(int));
557  //poly f_terms[f_len];
558  poly* f_terms=(poly*)omalloc(f_len*sizeof(poly));
559  poly** term_potences=(poly**) omalloc(f_len*sizeof(poly*));
560 
561  poly f_iter=f;
562  for(i=0;i<f_len;i++){
563    f_terms[i]=f_iter;
564    f_iter=pNext(f_iter);
565  }
566  for(i=0;i<f_len;i++){
567    term_potences[i]=(poly*)omalloc((n+1)*sizeof(poly));
568    term_potences[i][0]=p_ISet(1,r);
569    int j;
570    for(j=1;j<=n;j++){
571      term_potences[i][j]=p_MonMultCMB(f_terms[i],term_potences[i][j-1],r);
572    }
573  }
574  assume(f_iter==NULL);
575  poly zw=NULL;
576  poly tmp=p_Init(r);
577  number one=n_Init(1,r);
578  MC_iterate(f,n,r,f_len,&facult[0], &exp[0], &f_terms[0],erg_bucket,0,0,one/*facult[n]*/,zw,tmp, term_potences);
579
580
581  n_Delete(&one,r);
582
583
584 
585  //free res
586 
587  //free facult
588  for(i=0;i<=n;i++){
589    n_Delete(&facult[i],r);
590  }
591  int i2;
592  for (i=0;i<f_len;i++){
593    for(i2=0;i2<=n;i2++){
594      p_Delete(&term_potences[i][i2],r);
595    }
596    omfree(term_potences[i]);
597   
598  }
599  omfree(term_potences);
600  p_Delete(&tmp,r);
601  omfree(exp);
602  omfree(facult);
603  omfree(f_terms);
604  int len=0;
605  poly erg;
606  kBucketClear(erg_bucket,&erg, &len);
607  kBucketDestroy(&erg_bucket);
608  omUnGetSpecBin(&lm_bin);
609  return erg;
610}
Note: See TracBrowser for help on using the repository browser.