source: git/kernel/fast_mult.cc @ 8d1432e

spielwiese
Last change on this file since 8d1432e was 7b9b8e5, checked in by Hans Schoenemann <hannes@…>, 8 years ago
fix: Werror -> WerrorS if possible
  • Property mode set to 100644
File size: 13.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5#include <kernel/mod2.h>
6
7#include <polys/monomials/ring.h>
8#include <kernel/fast_mult.h>
9#include <polys/kbuckets.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}
465// unused
466#if 0
467static poly p_MonPowerMB(poly p, int exp, ring r)
468{
469  int i;
470
471  if(!n_IsOne(p_GetCoeff(p,r),r->cf))
472  {
473    number x, y;
474    y = p_GetCoeff(p,r);
475    n_Power(y,exp,&x,r->cf);
476    n_Delete(&y,r->cf);
477    p_SetCoeff0(p,x,r);
478  }
479  for (i=rVar(r); i!=0; i--)
480  {
481    p_MultExp(p,i, exp,r);
482  }
483  p_Setm(p,r);
484  return p;
485}
486#endif
487static 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){
488
489  int i;
490  //  poly term=p_Init(r);
491
492   //  number denom=n_Init(1,r);
493//   for(i=0;i<f_len;i++){
494//     if(exp[i]!=0){
495//       number trash=denom;
496//       denom=n_Mult(denom,facult[exp[i]],r);
497//       n_Delete(&trash,r);
498//     }
499
500//   }
501//   number coef=n_IntDiv(facult[n],denom,r);   //right function here?
502//   n_Delete(&denom,r);
503//   poly erg=p_NSet(coef,r);
504  poly erg=p_Init(r,lm_bin);
505  p_SetCoeff0(erg, coef,r);
506  //p_NSet(n_Copy(coef,r),r);
507  for(i=0;i<f_len;i++){
508    if(exp[i]!=0){
509      ///poly term=p_Copy(f_terms[i],r);
510      poly term=term_pot[i][exp[i]];
511        //tmp;
512        //p_ExpVectorCopy(term,f_terms[i],r);
513        //p_SetCoeff(term, n_Copy(p_GetCoeff(f_terms[i],r),r),r);
514
515      //term->next=NULL;
516
517      //p_MonPowerMB(term, exp[i],r);
518      p_MonMultMB(erg,term,r);
519      //p_Delete(&term,r);
520    }
521
522  }
523  zw=erg;
524}
525
526
527
528static 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){
529  int i;
530
531  if (pos<f_len-1)
532  {
533    poly zw_l=NULL;
534    number new_coef;
535    for(i=0;i<=n-sum;i++)
536    {
537      exp[pos]=i;
538      if(i==0)
539      {
540        new_coef=n_Copy(coef,r->cf);
541      }
542      else
543      {
544        number old=new_coef;
545        number old_rest=n_Init(n-sum-(i-1),r->cf);
546        new_coef=n_Mult(new_coef,old_rest,r->cf);
547        n_Delete(&old_rest,r->cf);
548        n_Delete(&old,r->cf);
549        number i_number=n_Init(i,r->cf);
550        old=new_coef;
551        new_coef=n_Div(new_coef,i_number,r->cf);
552        n_Normalize(new_coef,r->cf);
553        n_Delete(&old,r->cf);
554        n_Delete(&i_number,r->cf);
555      }
556      //new_coef is
557      //(n                          )
558      //(exp[0]..exp[pos] 0 0 0 rest)
559      poly zw_real=NULL;
560      MC_iterate(f, n, r, f_len,facult, exp,f_terms,erg_bucket,pos+1,sum+i,new_coef,zw_real,tmp,term_pot);
561      if (pos==f_len-2)
562      {
563        //get first small polys
564
565        zw_real->next=zw_l;
566        zw_l=zw_real;
567      }
568      //n_Delete(& new_coef,r->cf);
569    }
570    n_Delete(&new_coef,r->cf);
571    if (pos==f_len-2)
572    {
573      int len=n-sum+1;
574      kBucket_Add_q(erg_bucket,zw_l,&len);
575    }
576    return;
577  }
578  if(pos==f_len-1)
579  {
580    i=n-sum;
581    exp[pos]=i;
582    number new_coef=n_Copy(coef,r->cf);//n_IntDiv(coef,facult[i],r); //really consumed???????
583    buildTermAndAdd(n,facult,f_terms,exp,f_len,erg_bucket,r, new_coef,zw, tmp,term_pot);
584    // n_Delete(& new_coef,r);
585  }
586  assume(pos<=f_len-1);
587}
588poly pFastPowerMC(poly f, int n, ring r)
589{
590  //only char=0
591
592  if(rChar(r)!=0)
593    WerrorS("Char not 0, pFastPowerMC not implemented for this case");
594  if (n<=1)
595    WerrorS("not implemented for so small n, recursion fails");//should be length(f)
596   if (pLength(f)<=1)
597    WerrorS("not implemented for so small length of f, recursion fails");
598  //  number null_number=n_Init(0,r);
599  number* facult=(number*) omAlloc((n+1)*sizeof(number));
600  facult[0]=n_Init(1,r->cf);
601  int i;
602  for(i=1;i<=n;i++)
603  {
604    number this_n=n_Init(i,r->cf);
605    facult[i]=n_Mult(this_n,facult[i-1],r->cf);
606    n_Delete(&this_n,r->cf);
607  }
608
609  lm_bin=omGetSpecBin(POLYSIZE + (r->ExpL_Size)*sizeof(long));
610
611  kBucket_pt erg_bucket= kBucketCreate(currRing);
612  kBucketInit(erg_bucket,NULL,0);
613  const int f_len=pLength(f);
614  int* exp=(int*)omAlloc(f_len*sizeof(int));
615  //poly f_terms[f_len];
616  poly* f_terms=(poly*)omAlloc(f_len*sizeof(poly));
617  poly** term_potences=(poly**) omAlloc(f_len*sizeof(poly*));
618
619  poly f_iter=f;
620  for(i=0;i<f_len;i++)
621  {
622    f_terms[i]=f_iter;
623    f_iter=pNext(f_iter);
624  }
625  for(i=0;i<f_len;i++)
626  {
627    term_potences[i]=(poly*)omAlloc((n+1)*sizeof(poly));
628    term_potences[i][0]=p_ISet(1,r);
629    int j;
630    for(j=1;j<=n;j++){
631      term_potences[i][j]=p_MonMultCMB(f_terms[i],term_potences[i][j-1],r);
632    }
633  }
634  assume(f_iter==NULL);
635  poly zw=NULL;
636  poly tmp=p_Init(r);
637  number one=n_Init(1,r->cf);
638  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);
639
640
641  n_Delete(&one,r->cf);
642
643
644
645  //free res
646
647  //free facult
648  for(i=0;i<=n;i++)
649  {
650    n_Delete(&facult[i],r->cf);
651  }
652  int i2;
653  for (i=0;i<f_len;i++)
654  {
655    for(i2=0;i2<=n;i2++)
656    {
657      p_Delete(&term_potences[i][i2],r);
658    }
659    omfree(term_potences[i]);
660
661  }
662  omfree(term_potences);
663  p_Delete(&tmp,r);
664  omfree(exp);
665  omfree(facult);
666  omfree(f_terms);
667  int len=0;
668  poly erg;
669  kBucketClear(erg_bucket,&erg, &len);
670  kBucketDestroy(&erg_bucket);
671  omUnGetSpecBin(&lm_bin);
672  return erg;
673}
Note: See TracBrowser for help on using the repository browser.