source: git/kernel/fast_mult.cc @ fcb8022

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