source: git/kernel/fast_mult.cc @ abe5c8

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