source: git/kernel/fast_mult.cc @ 3a0e1a

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