[8716be] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[341696] | 4 | /* $Id$ */ |
---|
[1c4e9a] | 5 | #include "mod2.h" |
---|
[210e07] | 6 | #include <polys/monomials/ring.h> |
---|
[599326] | 7 | #include <kernel/fast_mult.h> |
---|
[210e07] | 8 | #include <polys/kbuckets.h> |
---|
[599326] | 9 | #include <kernel/febase.h> |
---|
[8716be] | 10 | |
---|
[78030c] | 11 | typedef poly fastmultrec(poly f, poly g, ring r); |
---|
[446c2d] | 12 | static const int pass_option=1; |
---|
[e95342] | 13 | static int mults=0; |
---|
[fca87a] | 14 | int Mults() |
---|
| 15 | { |
---|
[e95342] | 16 | return mults; |
---|
| 17 | } |
---|
[fca87a] | 18 | static void degsplit(poly p,int n,poly &p1,poly&p2, int vn, ring r) |
---|
| 19 | { |
---|
[7a05d36] | 20 | poly erg1_i, erg2_i; |
---|
| 21 | erg1_i=NULL; |
---|
| 22 | erg2_i=NULL; |
---|
[fca87a] | 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; |
---|
[7a05d36] | 34 | } |
---|
| 35 | erg1_i=p; |
---|
[fca87a] | 36 | } |
---|
| 37 | else |
---|
| 38 | { |
---|
| 39 | if (p2==NULL) |
---|
| 40 | { |
---|
| 41 | p2=p; |
---|
| 42 | } |
---|
| 43 | else |
---|
| 44 | { |
---|
| 45 | pNext(erg2_i)=p; |
---|
[7a05d36] | 46 | } |
---|
| 47 | erg2_i=p; |
---|
| 48 | } |
---|
| 49 | p=pNext(p); |
---|
| 50 | } |
---|
[fca87a] | 51 | if(erg2_i) |
---|
| 52 | { |
---|
[7a05d36] | 53 | pNext(erg2_i)=NULL; |
---|
[9cb2646] | 54 | } |
---|
[fca87a] | 55 | if (erg1_i) |
---|
| 56 | { |
---|
[7a05d36] | 57 | pNext(erg1_i)=NULL; |
---|
| 58 | } |
---|
[fca87a] | 59 | |
---|
[7a05d36] | 60 | } |
---|
[fca87a] | 61 | static void div_by_x_power_n(poly p, int n, int vn, ring r) |
---|
| 62 | { |
---|
| 63 | while(p) |
---|
| 64 | { |
---|
[78030c] | 65 | assume(p_GetExp(p,vn,r)>=n); |
---|
| 66 | int e=p_GetExp(p,vn,r); |
---|
| 67 | p_SetExp(p,vn,e-n,r); |
---|
[7a05d36] | 68 | p=pNext(p); |
---|
| 69 | } |
---|
| 70 | } |
---|
[cd5189] | 71 | |
---|
[fca87a] | 72 | static poly do_unifastmult(poly f,int df,poly g,int dg, int vn, fastmultrec rec, ring r) |
---|
| 73 | { |
---|
[7a05d36] | 74 | int n=1; |
---|
| 75 | if ((f==NULL)||(g==NULL)) return NULL; |
---|
[c70e58] | 76 | //int df=pGetExp(f,vn);//pFDeg(f); |
---|
| 77 | // int dg=pGetExp(g,vn);//pFDeg(g); |
---|
[fca87a] | 78 | |
---|
[7a05d36] | 79 | int dm; |
---|
[fca87a] | 80 | if(df>dg) |
---|
| 81 | { |
---|
[7a05d36] | 82 | dm=df; |
---|
[fca87a] | 83 | } |
---|
| 84 | else |
---|
| 85 | { |
---|
[7a05d36] | 86 | dm=dg; |
---|
| 87 | } |
---|
| 88 | while(n<=dm) |
---|
| 89 | { |
---|
| 90 | n*=2; |
---|
| 91 | } |
---|
[fca87a] | 92 | if(n==1) |
---|
| 93 | { |
---|
[78030c] | 94 | return(pp_Mult_qq(f,g,r)); |
---|
[7a05d36] | 95 | } |
---|
[fca87a] | 96 | |
---|
[7a05d36] | 97 | int pot=n/2; |
---|
| 98 | assume(pot*2==n); |
---|
[446c2d] | 99 | |
---|
| 100 | |
---|
| 101 | //splitting |
---|
[7a05d36] | 102 | poly f1=NULL; |
---|
| 103 | poly f0=NULL;//f-(x^(pot)*F1); |
---|
[78030c] | 104 | degsplit(p_Copy(f,r),pot,f1,f0,vn,r); |
---|
| 105 | div_by_x_power_n(f1,pot,vn,r); |
---|
[446c2d] | 106 | |
---|
[7a05d36] | 107 | poly g1=NULL; |
---|
| 108 | poly g0=NULL; |
---|
[78030c] | 109 | degsplit(p_Copy(g,r),pot,g1,g0,vn,r); |
---|
| 110 | div_by_x_power_n(g1,pot,vn,r); |
---|
[fca87a] | 111 | |
---|
[446c2d] | 112 | //p00, p11 |
---|
[78030c] | 113 | poly p00=rec(f0,g0,r);//unifastmult(f0,g0); |
---|
| 114 | poly p11=rec(f1,g1,r); |
---|
[446c2d] | 115 | |
---|
| 116 | //construct erg, factor |
---|
| 117 | poly erg=NULL; |
---|
[e95342] | 118 | poly factor=p_ISet(1,r); |
---|
[446c2d] | 119 | |
---|
[78030c] | 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); |
---|
[446c2d] | 123 | |
---|
| 124 | |
---|
| 125 | |
---|
[fca87a] | 126 | |
---|
| 127 | |
---|
| 128 | if((f1!=NULL) &&(f0!=NULL) &&(g0!=NULL) && (g1!=NULL)) |
---|
| 129 | { |
---|
[9cb2646] | 130 | //if(true){ |
---|
| 131 | //eat up f0,f1,g0,g1 |
---|
[78030c] | 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); |
---|
[446c2d] | 137 | |
---|
[fca87a] | 138 | |
---|
[9cb2646] | 139 | //eat up pbig |
---|
| 140 | poly sum=pbig; |
---|
[78030c] | 141 | p_SetExp(factor,vn,pot,r); |
---|
[fca87a] | 142 | |
---|
[9cb2646] | 143 | //eat up p00 |
---|
[78030c] | 144 | sum=p_Add_q(sum,p_Neg(p00,r),r); |
---|
[fca87a] | 145 | |
---|
[9cb2646] | 146 | //eat up p11 |
---|
[78030c] | 147 | sum=p_Add_q(sum,p_Neg(p11,r),r); |
---|
[7a05d36] | 148 | |
---|
[fca87a] | 149 | |
---|
[78030c] | 150 | sum=p_Mult_mm(sum,factor,r); |
---|
[446c2d] | 151 | |
---|
| 152 | |
---|
[9cb2646] | 153 | //eat up sum |
---|
[78030c] | 154 | erg=p_Add_q(sum,erg,r); |
---|
[fca87a] | 155 | } |
---|
| 156 | else |
---|
| 157 | { |
---|
| 158 | //eat up f0,f1,g0,g1 |
---|
[78030c] | 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); |
---|
[e95342] | 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); |
---|
[78030c] | 169 | erg=p_Add_q(erg,h,r); |
---|
[446c2d] | 170 | } |
---|
[fca87a] | 171 | |
---|
| 172 | |
---|
[78030c] | 173 | p_Delete(&factor,r); |
---|
[fca87a] | 174 | |
---|
| 175 | |
---|
[446c2d] | 176 | |
---|
[7a05d36] | 177 | return(erg); |
---|
| 178 | } |
---|
[cd5189] | 179 | |
---|
[c70e58] | 180 | // poly do_unifastmult_buckets(poly f,poly g){ |
---|
[fca87a] | 181 | |
---|
| 182 | |
---|
[446c2d] | 183 | |
---|
[cd5189] | 184 | |
---|
[c70e58] | 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); |
---|
[fca87a] | 189 | |
---|
[c70e58] | 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); |
---|
[fca87a] | 221 | |
---|
[c70e58] | 222 | // //p00 |
---|
| 223 | // //p11 |
---|
| 224 | // poly p00=unifastmult(f0,g0); |
---|
| 225 | // poly p11=unifastmult(f1,g1); |
---|
[fca87a] | 226 | |
---|
[446c2d] | 227 | |
---|
| 228 | |
---|
| 229 | |
---|
[c70e58] | 230 | // //eat up f0,f1,g0,g1 |
---|
[fca87a] | 231 | // poly pbig=unifastmult(pAdd(f0,f1),pAdd(g0,g1)); |
---|
[c70e58] | 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); |
---|
[446c2d] | 239 | |
---|
[c70e58] | 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); |
---|
[cd5189] | 249 | |
---|
[fca87a] | 250 | |
---|
[c70e58] | 251 | // pseudo_len=0; |
---|
[446c2d] | 252 | |
---|
[fca87a] | 253 | |
---|
[c70e58] | 254 | // pDelete(&factor); |
---|
[446c2d] | 255 | |
---|
[c70e58] | 256 | // int len=0; |
---|
| 257 | // poly erg=NULL; |
---|
| 258 | // kBucketClear(erg_bucket,&erg,&len); |
---|
| 259 | // kBucketDestroy(&erg_bucket); |
---|
[446c2d] | 260 | |
---|
[c70e58] | 261 | // return erg; |
---|
| 262 | // } |
---|
[9cb2646] | 263 | |
---|
[fca87a] | 264 | static inline int max(int a, int b) |
---|
| 265 | { |
---|
[9cb2646] | 266 | return (a>b)? a:b; |
---|
| 267 | } |
---|
[fca87a] | 268 | static inline int min(int a, int b) |
---|
| 269 | { |
---|
[9cb2646] | 270 | return (a>b)? b:a; |
---|
| 271 | } |
---|
[fca87a] | 272 | poly unifastmult(poly f,poly g, ring r) |
---|
| 273 | { |
---|
[c70e58] | 274 | int vn=1; |
---|
[cd5189] | 275 | if((f==NULL)||(g==NULL)) return NULL; |
---|
[78030c] | 276 | int df=p_GetExp(f,vn,r); |
---|
| 277 | int dg=p_GetExp(g,vn,r); |
---|
[c70e58] | 278 | if ((df==0)||(dg==0)) |
---|
[78030c] | 279 | return pp_Mult_qq(f,g,r); |
---|
[cd5189] | 280 | if (df*dg<100) |
---|
[78030c] | 281 | return pp_Mult_qq(f,g,r); |
---|
[cd5189] | 282 | // if (df*dg>10000) |
---|
| 283 | // return |
---|
| 284 | // do_unifastmult_buckets(f,g); |
---|
| 285 | //else |
---|
[78030c] | 286 | return do_unifastmult(f,df,g,dg,vn,unifastmult,r); |
---|
[cd5189] | 287 | |
---|
| 288 | } |
---|
| 289 | |
---|
[fca87a] | 290 | poly multifastmult(poly f, poly g, ring r) |
---|
| 291 | { |
---|
[e95342] | 292 | mults++; |
---|
[9cb2646] | 293 | if((f==NULL)||(g==NULL)) return NULL; |
---|
| 294 | if (pLength(f)*pLength(g)<100) |
---|
[78030c] | 295 | return pp_Mult_qq(f,g,r); |
---|
[9cb2646] | 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; |
---|
[fca87a] | 303 | for(i=1;i<=rVar(r);i++) |
---|
| 304 | { |
---|
[9cb2646] | 305 | poly p; |
---|
| 306 | int df=0; |
---|
| 307 | int dg=0; |
---|
| 308 | //max min max Strategie |
---|
| 309 | p=f; |
---|
[fca87a] | 310 | while(p) |
---|
| 311 | { |
---|
[78030c] | 312 | df=max(df,p_GetExp(p,i,r)); |
---|
[9cb2646] | 313 | p=pNext(p); |
---|
| 314 | } |
---|
[fca87a] | 315 | if(df>can_crit) |
---|
| 316 | { |
---|
[9cb2646] | 317 | p=g; |
---|
[fca87a] | 318 | while(p) |
---|
| 319 | { |
---|
| 320 | dg=max(dg,p_GetExp(p,i,r)); |
---|
| 321 | p=pNext(p); |
---|
[9cb2646] | 322 | } |
---|
| 323 | int crit=min(df,dg); |
---|
[fca87a] | 324 | if (crit>can_crit) |
---|
| 325 | { |
---|
| 326 | can_crit=crit; |
---|
| 327 | can_i=i; |
---|
| 328 | can_df=df; |
---|
| 329 | can_dg=dg; |
---|
[9cb2646] | 330 | } |
---|
| 331 | } |
---|
| 332 | } |
---|
| 333 | if(can_crit==0) |
---|
[78030c] | 334 | return pp_Mult_qq(f,g,r); |
---|
[9cb2646] | 335 | else |
---|
[78030c] | 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 | } |
---|
[9cb2646] | 341 | } |
---|
[fca87a] | 342 | poly pFastPower(poly f, int n, ring r) |
---|
| 343 | { |
---|
[955737] | 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; |
---|
[fca87a] | 349 | while(i_max*2<=n) |
---|
| 350 | { |
---|
[955737] | 351 | i_max*=2; |
---|
| 352 | pot_max++; |
---|
| 353 | } |
---|
| 354 | int field_size=pot_max+1; |
---|
[fca87a] | 355 | int* int_pot_array=(int*) omAlloc(field_size*sizeof(int)); |
---|
| 356 | poly* pot_array=(poly*) omAlloc(field_size*sizeof(poly)); |
---|
[955737] | 357 | int i; |
---|
| 358 | int pot=1; |
---|
| 359 | //initializing int_pot |
---|
[fca87a] | 360 | for(i=0;i<field_size;i++) |
---|
| 361 | { |
---|
[955737] | 362 | int_pot_array[i]=pot; |
---|
| 363 | pot*=2; |
---|
| 364 | } |
---|
| 365 | //calculating pot_array |
---|
| 366 | pot_array[0]=f; //do not delete it |
---|
[fca87a] | 367 | for(i=1;i<field_size;i++) |
---|
| 368 | { |
---|
[955737] | 369 | poly p=pot_array[i-1]; |
---|
[c99285f] | 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); |
---|
[955737] | 374 | } |
---|
[fca87a] | 375 | |
---|
[955737] | 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){ |
---|
[fca87a] | 388 | // work_n-=int_pot_array[i]; |
---|
| 389 | // poly prod=multifastmult(erg,pot_array[i],r); |
---|
| 390 | // pDelete(&erg); |
---|
| 391 | // erg=prod; |
---|
[955737] | 392 | // } |
---|
| 393 | |
---|
| 394 | // if(i!=0) pDelete(&pot_array[i]); |
---|
| 395 | // } |
---|
[fca87a] | 396 | |
---|
| 397 | |
---|
| 398 | for(i=field_size-1;i>=0;i--) |
---|
| 399 | { |
---|
[955737] | 400 | |
---|
| 401 | assume(work_n<2*int_pot_array[i]); |
---|
[fca87a] | 402 | if(int_pot_array[i]<=work_n) |
---|
| 403 | { |
---|
| 404 | work_n-=int_pot_array[i]; |
---|
| 405 | int_pot_array[i]=1; |
---|
[955737] | 406 | } |
---|
| 407 | else int_pot_array[i]=0; |
---|
[fca87a] | 408 | |
---|
[955737] | 409 | } |
---|
[fca87a] | 410 | for(i=0;i<field_size;i++) |
---|
| 411 | { |
---|
| 412 | if(int_pot_array[i]==1) |
---|
| 413 | { |
---|
[c99285f] | 414 | poly prod; |
---|
| 415 | if(rVar(r)==1) |
---|
[fca87a] | 416 | prod=multifastmult(erg,pot_array[i],r); |
---|
[c99285f] | 417 | else |
---|
[fca87a] | 418 | prod=pp_Mult_qq(erg,pot_array[i],r); |
---|
[955737] | 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 | } |
---|
[cbea799] | 429 | static omBin lm_bin=NULL; |
---|
| 430 | /*3 |
---|
| 431 | * compute for monomials p*q |
---|
| 432 | * destroys p, keeps q |
---|
| 433 | */ |
---|
| 434 | static void p_MonMultMB(poly p, poly q,ring r) |
---|
| 435 | { |
---|
| 436 | number x, y; |
---|
| 437 | // int i; |
---|
| 438 | |
---|
| 439 | y = p_GetCoeff(p,r); |
---|
[1c4e9a] | 440 | x = n_Mult(y,pGetCoeff(q),r->cf); |
---|
| 441 | n_Delete(&y,r->cf); |
---|
[cbea799] | 442 | p_SetCoeff0(p,x,r); |
---|
[1f637e] | 443 | //for (i=(currRing->N); i!=0; i--) |
---|
[cbea799] | 444 | //{ |
---|
| 445 | // pAddExp(p,i, pGetExp(q,i)); |
---|
| 446 | //} |
---|
| 447 | //p->Order += q->Order; |
---|
| 448 | p_ExpVectorAdd(p,q,r); |
---|
| 449 | } |
---|
[c94036] | 450 | /*3 |
---|
| 451 | * compute for monomials p*q |
---|
| 452 | * keeps p, q |
---|
| 453 | * uses bin only available in MCPower |
---|
| 454 | */ |
---|
| 455 | static poly p_MonMultCMB(poly p, poly q, ring r) |
---|
| 456 | { |
---|
| 457 | number x; |
---|
| 458 | poly res = p_Init(r,lm_bin); |
---|
| 459 | |
---|
[1c4e9a] | 460 | x = n_Mult(p_GetCoeff(p,r),p_GetCoeff(q,r),r->cf); |
---|
[c94036] | 461 | p_SetCoeff0(res,x,r); |
---|
| 462 | p_ExpVectorSum(res,p, q,r); |
---|
| 463 | return res; |
---|
| 464 | } |
---|
[040ee80] | 465 | static poly p_MonPowerMB(poly p, int exp, ring r) |
---|
| 466 | { |
---|
| 467 | int i; |
---|
| 468 | |
---|
[1c4e9a] | 469 | if(!n_IsOne(p_GetCoeff(p,r),r->cf)) |
---|
[040ee80] | 470 | { |
---|
| 471 | number x, y; |
---|
| 472 | y = p_GetCoeff(p,r); |
---|
[1c4e9a] | 473 | n_Power(y,exp,&x,r->cf); |
---|
| 474 | n_Delete(&y,r->cf); |
---|
[040ee80] | 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 | } |
---|
[c94036] | 484 | static 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){ |
---|
[fca87a] | 485 | |
---|
[040ee80] | 486 | int i; |
---|
[cbea799] | 487 | // poly term=p_Init(r); |
---|
[d2d22f] | 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 | // } |
---|
[fca87a] | 496 | |
---|
[d2d22f] | 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); |
---|
[cbea799] | 501 | poly erg=p_Init(r,lm_bin); |
---|
| 502 | p_SetCoeff0(erg, coef,r); |
---|
| 503 | //p_NSet(n_Copy(coef,r),r); |
---|
[040ee80] | 504 | for(i=0;i<f_len;i++){ |
---|
| 505 | if(exp[i]!=0){ |
---|
[cbea799] | 506 | ///poly term=p_Copy(f_terms[i],r); |
---|
[c94036] | 507 | poly term=term_pot[i][exp[i]]; |
---|
[fca87a] | 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 | |
---|
[cbea799] | 512 | //term->next=NULL; |
---|
[fca87a] | 513 | |
---|
[c94036] | 514 | //p_MonPowerMB(term, exp[i],r); |
---|
[cbea799] | 515 | p_MonMultMB(erg,term,r); |
---|
| 516 | //p_Delete(&term,r); |
---|
[040ee80] | 517 | } |
---|
[fca87a] | 518 | |
---|
[040ee80] | 519 | } |
---|
[cbea799] | 520 | zw=erg; |
---|
[040ee80] | 521 | } |
---|
| 522 | |
---|
| 523 | |
---|
| 524 | |
---|
[c94036] | 525 | static 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){ |
---|
[040ee80] | 526 | int i; |
---|
[fca87a] | 527 | |
---|
| 528 | if (pos<f_len-1) |
---|
| 529 | { |
---|
[cbea799] | 530 | poly zw_l=NULL; |
---|
[9a92d97] | 531 | number new_coef; |
---|
[fca87a] | 532 | for(i=0;i<=n-sum;i++) |
---|
| 533 | { |
---|
[040ee80] | 534 | exp[pos]=i; |
---|
[fca87a] | 535 | if(i==0) |
---|
| 536 | { |
---|
[1c4e9a] | 537 | new_coef=n_Copy(coef,r->cf); |
---|
[9a92d97] | 538 | } |
---|
[fca87a] | 539 | else |
---|
| 540 | { |
---|
| 541 | number old=new_coef; |
---|
[1c4e9a] | 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); |
---|
[fca87a] | 547 | old=new_coef; |
---|
[1c4e9a] | 548 | new_coef=n_IntDiv(new_coef,i_number,r->cf); |
---|
| 549 | n_Delete(&old,r->cf); |
---|
| 550 | n_Delete(&i_number,r->cf); |
---|
[9a92d97] | 551 | } |
---|
[fca87a] | 552 | //new_coef is |
---|
[9a92d97] | 553 | //(n ) |
---|
| 554 | //(exp[0]..exp[pos] 0 0 0 rest) |
---|
[cbea799] | 555 | poly zw_real=NULL; |
---|
[c94036] | 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); |
---|
[fca87a] | 557 | if (pos==f_len-2) |
---|
| 558 | { |
---|
| 559 | //get first small polys |
---|
| 560 | |
---|
| 561 | zw_real->next=zw_l; |
---|
| 562 | zw_l=zw_real; |
---|
[cbea799] | 563 | } |
---|
[1c4e9a] | 564 | //n_Delete(& new_coef,r->cf); |
---|
[040ee80] | 565 | } |
---|
[1c4e9a] | 566 | n_Delete(&new_coef,r->cf); |
---|
[fca87a] | 567 | if (pos==f_len-2) |
---|
| 568 | { |
---|
[cbea799] | 569 | int len=n-sum+1; |
---|
| 570 | kBucket_Add_q(erg_bucket,zw_l,&len); |
---|
| 571 | } |
---|
[040ee80] | 572 | return; |
---|
| 573 | } |
---|
[fca87a] | 574 | if(pos==f_len-1) |
---|
| 575 | { |
---|
[040ee80] | 576 | i=n-sum; |
---|
| 577 | exp[pos]=i; |
---|
[1c4e9a] | 578 | number new_coef=n_Copy(coef,r->cf);//n_IntDiv(coef,facult[i],r); //really consumed??????? |
---|
[c94036] | 579 | buildTermAndAdd(n,facult,f_terms,exp,f_len,erg_bucket,r, new_coef,zw, tmp,term_pot); |
---|
[cbea799] | 580 | // n_Delete(& new_coef,r); |
---|
[040ee80] | 581 | } |
---|
| 582 | assume(pos<=f_len-1); |
---|
| 583 | } |
---|
[fca87a] | 584 | poly pFastPowerMC(poly f, int n, ring r) |
---|
| 585 | { |
---|
[040ee80] | 586 | //only char=0 |
---|
[fca87a] | 587 | |
---|
[040ee80] | 588 | if(rChar(r)!=0) |
---|
| 589 | Werror("Char not 0, pFastPowerMC not implemented for this case"); |
---|
[cbea799] | 590 | if (n<=1) |
---|
[1903b7] | 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"); |
---|
[040ee80] | 594 | // number null_number=n_Init(0,r); |
---|
[fca87a] | 595 | number* facult=(number*) omAlloc((n+1)*sizeof(number)); |
---|
[1c4e9a] | 596 | facult[0]=n_Init(1,r->cf); |
---|
[040ee80] | 597 | int i; |
---|
[fca87a] | 598 | for(i=1;i<=n;i++) |
---|
| 599 | { |
---|
[1c4e9a] | 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); |
---|
[040ee80] | 603 | } |
---|
[cbea799] | 604 | |
---|
| 605 | lm_bin=omGetSpecBin(POLYSIZE + (r->ExpL_Size)*sizeof(long)); |
---|
| 606 | |
---|
[040ee80] | 607 | kBucket_pt erg_bucket= kBucketCreate(currRing); |
---|
| 608 | kBucketInit(erg_bucket,NULL,0); |
---|
| 609 | const int f_len=pLength(f); |
---|
[fca87a] | 610 | int* exp=(int*)omAlloc(f_len*sizeof(int)); |
---|
[040ee80] | 611 | //poly f_terms[f_len]; |
---|
[fca87a] | 612 | poly* f_terms=(poly*)omAlloc(f_len*sizeof(poly)); |
---|
| 613 | poly** term_potences=(poly**) omAlloc(f_len*sizeof(poly*)); |
---|
| 614 | |
---|
[040ee80] | 615 | poly f_iter=f; |
---|
[fca87a] | 616 | for(i=0;i<f_len;i++) |
---|
| 617 | { |
---|
[040ee80] | 618 | f_terms[i]=f_iter; |
---|
| 619 | f_iter=pNext(f_iter); |
---|
| 620 | } |
---|
[fca87a] | 621 | for(i=0;i<f_len;i++) |
---|
| 622 | { |
---|
| 623 | term_potences[i]=(poly*)omAlloc((n+1)*sizeof(poly)); |
---|
[c94036] | 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 | } |
---|
[040ee80] | 630 | assume(f_iter==NULL); |
---|
[cbea799] | 631 | poly zw=NULL; |
---|
| 632 | poly tmp=p_Init(r); |
---|
[1c4e9a] | 633 | number one=n_Init(1,r->cf); |
---|
[9a92d97] | 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); |
---|
[040ee80] | 635 | |
---|
| 636 | |
---|
[1c4e9a] | 637 | n_Delete(&one,r->cf); |
---|
[040ee80] | 638 | |
---|
| 639 | |
---|
[fca87a] | 640 | |
---|
[040ee80] | 641 | //free res |
---|
[fca87a] | 642 | |
---|
[040ee80] | 643 | //free facult |
---|
[fca87a] | 644 | for(i=0;i<=n;i++) |
---|
| 645 | { |
---|
[1c4e9a] | 646 | n_Delete(&facult[i],r->cf); |
---|
[040ee80] | 647 | } |
---|
[c94036] | 648 | int i2; |
---|
[fca87a] | 649 | for (i=0;i<f_len;i++) |
---|
| 650 | { |
---|
| 651 | for(i2=0;i2<=n;i2++) |
---|
| 652 | { |
---|
[c94036] | 653 | p_Delete(&term_potences[i][i2],r); |
---|
| 654 | } |
---|
| 655 | omfree(term_potences[i]); |
---|
[fca87a] | 656 | |
---|
[c94036] | 657 | } |
---|
| 658 | omfree(term_potences); |
---|
[cbea799] | 659 | p_Delete(&tmp,r); |
---|
[040ee80] | 660 | omfree(exp); |
---|
| 661 | omfree(facult); |
---|
| 662 | omfree(f_terms); |
---|
| 663 | int len=0; |
---|
| 664 | poly erg; |
---|
| 665 | kBucketClear(erg_bucket,&erg, &len); |
---|
[d2d22f] | 666 | kBucketDestroy(&erg_bucket); |
---|
[cbea799] | 667 | omUnGetSpecBin(&lm_bin); |
---|
[040ee80] | 668 | return erg; |
---|
| 669 | } |
---|