![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to up.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to up.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/up.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3: #include <math.h>
! 4:
! 5: struct oEGT eg_chrem,eg_fore,eg_back;
! 6:
! 7: /*
! 8: int up_kara_mag=15;
! 9: int up_tkara_mag=15;
! 10: */
! 11: /*
! 12: int up_kara_mag=30;
! 13: int up_tkara_mag=20;
! 14: */
! 15: int up_kara_mag=20;
! 16: int up_tkara_mag=15;
! 17:
! 18: #if defined(sparc)
! 19: int up_fft_mag=50;
! 20: #else
! 21: int up_fft_mag=80;
! 22: #endif
! 23:
! 24: int up_rem_mag=1;
! 25: int debug_up;
! 26: int up_lazy;
! 27: extern int lm_lazy;
! 28: extern int current_ff;
! 29: extern int GC_dont_gc;
! 30:
! 31: void monicup(a,b)
! 32: UP a;
! 33: UP *b;
! 34: {
! 35: UP w;
! 36:
! 37: if ( !a )
! 38: *b = 0;
! 39: else {
! 40: w = W_UPALLOC(0); w->d = 0;
! 41: divnum(0,(Num)ONE,(Num)a->c[a->d],(Num *)&w->c[0]);
! 42: mulup(a,w,b);
! 43: }
! 44: }
! 45:
! 46: void simpup(a,b)
! 47: UP a;
! 48: UP *b;
! 49: {
! 50: int i,d;
! 51: UP c;
! 52:
! 53: if ( !a )
! 54: *b = 0;
! 55: else {
! 56: d = a->d;
! 57: c = UPALLOC(d);
! 58:
! 59: for ( i = 0; i <= d; i++ )
! 60: simpnum(a->c[i],&c->c[i]);
! 61: for ( i = d; i >= 0 && !c->c[i]; i-- );
! 62: if ( i < 0 )
! 63: *b = 0;
! 64: else {
! 65: c->d = i;
! 66: *b = c;
! 67: }
! 68: }
! 69: }
! 70:
! 71: void simpnum(a,b)
! 72: Num a;
! 73: Num *b;
! 74: {
! 75: LM lm;
! 76: GF2N gf;
! 77: GFPN gfpn;
! 78:
! 79: if ( !a )
! 80: *b = 0;
! 81: else
! 82: switch ( NID(a) ) {
! 83: case N_LM:
! 84: simplm((LM)a,&lm); *b = (Num)lm; break;
! 85: case N_GF2N:
! 86: simpgf2n((GF2N)a,&gf); *b = (Num)gf; break;
! 87: case N_GFPN:
! 88: simpgfpn((GFPN)a,&gfpn); *b = (Num)gfpn; break;
! 89: default:
! 90: *b = a; break;
! 91: }
! 92: }
! 93:
! 94: void uremp(p1,p2,rp)
! 95: P p1,p2;
! 96: P *rp;
! 97: {
! 98: UP n1,n2,r;
! 99:
! 100: if ( !p1 || NUM(p1) )
! 101: *rp = p1;
! 102: else {
! 103: ptoup(p1,&n1); ptoup(p2,&n2);
! 104: remup(n1,n2,&r);
! 105: uptop(r,rp);
! 106: }
! 107: }
! 108:
! 109: void ugcdp(p1,p2,rp)
! 110: P p1,p2;
! 111: P *rp;
! 112: {
! 113: UP n1,n2,r;
! 114:
! 115: ptoup(p1,&n1); ptoup(p2,&n2);
! 116: gcdup(n1,n2,&r);
! 117: uptop(r,rp);
! 118: }
! 119:
! 120: void reversep(p1,d,rp)
! 121: P p1;
! 122: Q d;
! 123: P *rp;
! 124: {
! 125: UP n1,r;
! 126:
! 127: if ( !p1 )
! 128: *rp = 0;
! 129: else {
! 130: ptoup(p1,&n1);
! 131: reverseup(n1,QTOS(d),&r);
! 132: uptop(r,rp);
! 133: }
! 134: }
! 135:
! 136: void invmodp(p1,d,rp)
! 137: P p1;
! 138: Q d;
! 139: P *rp;
! 140: {
! 141: UP n1,r;
! 142:
! 143: if ( !p1 )
! 144: error("invmodp : invalid input");
! 145: else {
! 146: ptoup(p1,&n1);
! 147: invmodup(n1,QTOS(d),&r);
! 148: uptop(r,rp);
! 149: }
! 150: }
! 151:
! 152: void addup(n1,n2,nr)
! 153: UP n1,n2;
! 154: UP *nr;
! 155: {
! 156: UP r,t;
! 157: int i,d1,d2;
! 158: Num *c,*c1,*c2;
! 159:
! 160: if ( !n1 ) {
! 161: *nr = n2;
! 162: } else if ( !n2 ) {
! 163: *nr = n1;
! 164: } else {
! 165: if ( n2->d > n1->d ) {
! 166: t = n1; n1 = n2; n2 = t;
! 167: }
! 168: d1 = n1->d; d2 = n2->d;
! 169: *nr = r = UPALLOC(d1);
! 170: c1 = n1->c; c2 = n2->c; c = r->c;
! 171: for ( i = 0; i <= d2; i++ )
! 172: addnum(0,*c1++,*c2++,c++);
! 173: for ( ; i <= d1; i++ )
! 174: *c++ = *c1++;
! 175: for ( i = d1, c = r->c; i >=0 && !c[i]; i-- );
! 176: if ( i < 0 )
! 177: *nr = 0;
! 178: else
! 179: r->d = i;
! 180: }
! 181: }
! 182:
! 183: void subup(n1,n2,nr)
! 184: UP n1,n2;
! 185: UP *nr;
! 186: {
! 187: UP r;
! 188: int i,d1,d2,d;
! 189: Num *c,*c1,*c2;
! 190:
! 191: if ( !n1 ) {
! 192: chsgnup(n2,nr);
! 193: } else if ( !n2 ) {
! 194: *nr = n1;
! 195: } else {
! 196: d1 = n1->d; d2 = n2->d; d = MAX(d1,d2);
! 197: *nr = r = UPALLOC(d);
! 198: c1 = n1->c; c2 = n2->c; c = r->c;
! 199: if ( d1 >= d2 ) {
! 200: for ( i = 0; i <= d2; i++ )
! 201: subnum(0,*c1++,*c2++,c++);
! 202: for ( ; i <= d1; i++ )
! 203: *c++ = *c1++;
! 204: } else {
! 205: for ( i = 0; i <= d1; i++ )
! 206: subnum(0,*c1++,*c2++,c++);
! 207: for ( ; i <= d2; i++ )
! 208: chsgnnum(*c2++,c++);
! 209: }
! 210: for ( i = d, c = r->c; i >=0 && !c[i]; i-- );
! 211: if ( i < 0 )
! 212: *nr = 0;
! 213: else
! 214: r->d = i;
! 215: }
! 216: }
! 217:
! 218: void chsgnup(n1,nr)
! 219: UP n1;
! 220: UP *nr;
! 221: {
! 222: UP r;
! 223: int d1,i;
! 224: Num *c1,*c;
! 225:
! 226: if ( !n1 ) {
! 227: *nr = 0;
! 228: } else {
! 229: d1 = n1->d;
! 230: *nr = r = UPALLOC(d1); r->d = d1;
! 231: c1 = n1->c; c = r->c;
! 232: for ( i = 0; i <= d1; i++ )
! 233: chsgnnum(*c1++,c++);
! 234: }
! 235: }
! 236:
! 237: void hybrid_mulup(ff,n1,n2,nr)
! 238: int ff;
! 239: UP n1,n2;
! 240: UP *nr;
! 241: {
! 242: if ( !n1 || !n2 )
! 243: *nr = 0;
! 244: else if ( MAX(n1->d,n2->d) < up_fft_mag )
! 245: kmulup(n1,n2,nr);
! 246: else
! 247: switch ( ff ) {
! 248: case FF_GFP:
! 249: fft_mulup_lm(n1,n2,nr); break;
! 250: case FF_GF2N:
! 251: kmulup(n1,n2,nr); break;
! 252: default:
! 253: fft_mulup(n1,n2,nr); break;
! 254: }
! 255: }
! 256:
! 257: void hybrid_squareup(ff,n1,nr)
! 258: int ff;
! 259: UP n1;
! 260: UP *nr;
! 261: {
! 262: if ( !n1 )
! 263: *nr = 0;
! 264: else if ( n1->d < up_fft_mag )
! 265: ksquareup(n1,nr);
! 266: else
! 267: switch ( ff ) {
! 268: case FF_GFP:
! 269: fft_squareup_lm(n1,nr); break;
! 270: case FF_GF2N:
! 271: ksquareup(n1,nr); break;
! 272: default:
! 273: fft_squareup(n1,nr); break;
! 274: }
! 275: }
! 276:
! 277: void hybrid_tmulup(ff,n1,n2,d,nr)
! 278: int ff;
! 279: UP n1,n2;
! 280: int d;
! 281: UP *nr;
! 282: {
! 283: if ( !n1 || !n2 )
! 284: *nr = 0;
! 285: else if ( MAX(n1->d,n2->d) < up_fft_mag )
! 286: tkmulup(n1,n2,d,nr);
! 287: else
! 288: switch ( ff ) {
! 289: case FF_GFP:
! 290: trunc_fft_mulup_lm(n1,n2,d,nr); break;
! 291: case FF_GF2N:
! 292: tkmulup(n1,n2,d,nr); break;
! 293: default:
! 294: trunc_fft_mulup(n1,n2,d,nr); break;
! 295: }
! 296: }
! 297:
! 298: void mulup(n1,n2,nr)
! 299: UP n1,n2;
! 300: UP *nr;
! 301: {
! 302: UP r;
! 303: Num *pc1,*pc,*c1,*c2,*c;
! 304: Num mul,t,s;
! 305: int i,j,d1,d2;
! 306:
! 307: if ( !n1 || !n2 )
! 308: *nr = 0;
! 309: else {
! 310: d1 = n1->d; d2 = n2->d;
! 311: *nr = r = UPALLOC(d1+d2); r->d = d1+d2;
! 312: c1 = n1->c; c2 = n2->c; c = r->c;
! 313: lm_lazy = 1;
! 314: for ( i = 0; i <= d2; i++, c++ )
! 315: if ( mul = *c2++ )
! 316: for ( j = 0, pc1 = c1, pc = c; j <= d1; j++ ) {
! 317: mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
! 318: }
! 319: lm_lazy = 0;
! 320: if ( !up_lazy )
! 321: for ( i = 0, c = r->c; i <= r->d; i++ ) {
! 322: simpnum(c[i],&t); c[i] = t;
! 323: }
! 324: }
! 325: }
! 326:
! 327: /* nr = c*n1 */
! 328:
! 329: void mulcup(c,n1,nr)
! 330: Num c;
! 331: UP n1;
! 332: UP *nr;
! 333: {
! 334: int d;
! 335: UP r;
! 336: Num t;
! 337: Num *c1,*cr;
! 338: int i;
! 339:
! 340: if ( !c || !n1 )
! 341: *nr = 0;
! 342: else {
! 343: d = n1->d;
! 344: *nr = r = UPALLOC(d); r->d = d;
! 345: c1 = n1->c; cr = r->c;
! 346: lm_lazy = 1;
! 347: for ( i = 0; i <= d; i++ )
! 348: mulnum(0,c1[i],c,&cr[i]);
! 349: lm_lazy = 0;
! 350: if ( !up_lazy )
! 351: for ( i = 0; i <= d; i++ ) {
! 352: simpnum(cr[i],&t); cr[i] = t;
! 353: }
! 354: }
! 355: }
! 356:
! 357: void tmulup(n1,n2,d,nr)
! 358: UP n1,n2;
! 359: int d;
! 360: UP *nr;
! 361: {
! 362: UP r;
! 363: Num *pc1,*pc,*c1,*c2,*c;
! 364: Num mul,t,s;
! 365: int i,j,d1,d2,dr;
! 366:
! 367: if ( !n1 || !n2 )
! 368: *nr = 0;
! 369: else {
! 370: d1 = n1->d; d2 = n2->d;
! 371: dr = MAX(d-1,d1+d2);
! 372: *nr = r = UPALLOC(dr);
! 373: c1 = n1->c; c2 = n2->c; c = r->c;
! 374: lm_lazy = 1;
! 375: for ( i = 0; i <= d2; i++, c++ )
! 376: if ( mul = *c2++ )
! 377: for ( j = 0, pc1 = c1, pc = c; (j <= d1) && (i+j < d); j++ ) {
! 378: mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
! 379: }
! 380: lm_lazy = 0;
! 381: c = r->c;
! 382: if ( !up_lazy )
! 383: for ( i = 0; i < d; i++ ) {
! 384: simpnum(c[i],&t); c[i] = t;
! 385: }
! 386: for ( i = d-1; i >= 0 && !c[i]; i-- );
! 387: if ( i < 0 )
! 388: *nr = 0;
! 389: else
! 390: r->d = i;
! 391: }
! 392: }
! 393:
! 394: void squareup(n1,nr)
! 395: UP n1;
! 396: UP *nr;
! 397: {
! 398: UP r;
! 399: Num *c1,*c;
! 400: Num t,s,u;
! 401: int i,j,h,d1,d;
! 402:
! 403: if ( !n1 )
! 404: *nr = 0;
! 405: else if ( !n1->d ) {
! 406: *nr = r = UPALLOC(0); r->d = 0;
! 407: mulnum(0,n1->c[0],n1->c[0],&r->c[0]);
! 408: } else {
! 409: d1 = n1->d;
! 410: d = 2*d1;
! 411: *nr = r = UPALLOC(2*d); r->d = d;
! 412: c1 = n1->c; c = r->c;
! 413: lm_lazy = 1;
! 414: for ( i = 0; i <= d1; i++ )
! 415: mulnum(0,c1[i],c1[i],&c[2*i]);
! 416: for ( j = 1; j < d; j++ ) {
! 417: h = (j-1)/2;
! 418: for ( i = MAX(j-d1,0), s = 0; i <= h; i++ ) {
! 419: mulnum(0,c1[i],c1[j-i],&t); addnum(0,t,s,&u); s = u;
! 420: }
! 421: addnum(0,s,s,&t); addnum(0,t,c[j],&u); c[j] = u;
! 422: }
! 423: lm_lazy = 0;
! 424: if ( !up_lazy )
! 425: for ( i = 0, c = r->c; i <= d; i++ ) {
! 426: simpnum(c[i],&t); c[i] = t;
! 427: }
! 428: }
! 429: }
! 430:
! 431: void remup(n1,n2,nr)
! 432: UP n1,n2;
! 433: UP *nr;
! 434: {
! 435: UP w,r;
! 436:
! 437: if ( !n2 )
! 438: error("remup : division by 0");
! 439: else if ( !n1 || !n2->d )
! 440: *nr = 0;
! 441: else if ( n1->d < n2->d )
! 442: *nr = n1;
! 443: else {
! 444: w = W_UPALLOC(n1->d); copyup(n1,w);
! 445: remup_destructive(w,n2);
! 446: if ( w->d < 0 )
! 447: *nr = 0;
! 448: else {
! 449: *nr = r = UPALLOC(w->d); copyup(w,r);
! 450: }
! 451: }
! 452: }
! 453:
! 454: void remup_destructive(n1,n2)
! 455: UP n1,n2;
! 456: {
! 457: Num *c1,*c2;
! 458: int d1,d2,i,j;
! 459: Num m,t,s,mhc;
! 460:
! 461: c1 = n1->c; c2 = n2->c;
! 462: d1 = n1->d; d2 = n2->d;
! 463: chsgnnum(c2[d2],&mhc);
! 464: for ( i = d1; i >= d2; i-- ) {
! 465: simpnum(c1[i],&t); c1[i] = t;
! 466: if ( !c1[i] )
! 467: continue;
! 468: divnum(0,c1[i],mhc,&m);
! 469: lm_lazy = 1;
! 470: for ( j = d2; j >= 0; j-- ) {
! 471: mulnum(0,c2[j],m,&t); addnum(0,c1[i+j-d2],t,&s);
! 472: c1[i+j-d2] = s;
! 473: }
! 474: lm_lazy = 0;
! 475: }
! 476: for ( i = 0; i < d2; i++ ) {
! 477: simpnum(c1[i],&t); c1[i] = t;
! 478: }
! 479: for ( i = d2-1; i >= 0 && !c1[i]; i-- );
! 480: n1->d = i;
! 481: }
! 482:
! 483: void qrup(n1,n2,nq,nr)
! 484: UP n1,n2;
! 485: UP *nq,*nr;
! 486: {
! 487: UP w,r,q;
! 488: struct oUP t;
! 489: Num m;
! 490: int d;
! 491:
! 492: if ( !n2 )
! 493: error("qrup : division by 0");
! 494: else if ( !n1 ) {
! 495: *nq = 0; *nr = 0;
! 496: } else if ( !n2->d )
! 497: if ( !n1->d ) {
! 498: divnum(0,n1->c[0],n2->c[0],(Num *)nq); *nr = 0;
! 499: } else {
! 500: divnum(0,(Num)ONE,n2->c[0],&m);
! 501: t.d = 0; t.c[0] = m;
! 502: mulup(n1,&t,nq); *nr = 0;
! 503: }
! 504: else if ( n1->d < n2->d ) {
! 505: *nq = 0; *nr = n1;
! 506: } else {
! 507: w = W_UPALLOC(n1->d); copyup(n1,w);
! 508: qrup_destructive(w,n2);
! 509: d = n1->d-n2->d;
! 510: *nq = q = UPALLOC(d); q->d = d;
! 511: bcopy(w->c+n2->d,q->c,(d+1)*sizeof(Num));
! 512: if ( w->d < 0 )
! 513: *nr = 0;
! 514: else {
! 515: *nr = r = UPALLOC(w->d); copyup(w,r);
! 516: }
! 517: }
! 518: }
! 519:
! 520: void qrup_destructive(n1,n2)
! 521: UP n1,n2;
! 522: {
! 523: Num *c1,*c2;
! 524: int d1,d2,i,j;
! 525: Num q,t,s,hc;
! 526:
! 527: c1 = n1->c; c2 = n2->c;
! 528: d1 = n1->d; d2 = n2->d;
! 529: hc = c2[d2];
! 530: for ( i = d1; i >= d2; i-- ) {
! 531: simpnum(c1[i],&t); c1[i] = t;
! 532: if ( c1[i] ) {
! 533: divnum(0,c1[i],hc,&q);
! 534: lm_lazy = 1;
! 535: for ( j = d2; j >= 0; j-- ) {
! 536: mulnum(0,c2[j],q,&t); subnum(0,c1[i+j-d2],t,&s);
! 537: c1[i+j-d2] = s;
! 538: }
! 539: lm_lazy = 0;
! 540: } else
! 541: q = 0;
! 542: c1[i] = q;
! 543: }
! 544: for ( i = 0; i < d2; i++ ) {
! 545: simpnum(c1[i],&t); c1[i] = t;
! 546: }
! 547: for ( i = d2-1; i >= 0 && !c1[i]; i-- );
! 548: n1->d = i;
! 549: }
! 550:
! 551: void gcdup(n1,n2,nr)
! 552: UP n1,n2;
! 553: UP *nr;
! 554: {
! 555: UP w1,w2,t,r;
! 556: int d1,d2;
! 557:
! 558: if ( !n1 )
! 559: *nr = n2;
! 560: else if ( !n2 )
! 561: *nr = n1;
! 562: else if ( !n1->d || !n2->d ) {
! 563: *nr = r = UPALLOC(0); r->d = 0;
! 564: divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]);
! 565: } else {
! 566: d1 = n1->d; d2 = n2->d;
! 567: if ( d2 > d1 ) {
! 568: w1 = W_UPALLOC(d2); copyup(n2,w1);
! 569: w2 = W_UPALLOC(d1); copyup(n1,w2);
! 570: } else {
! 571: w1 = W_UPALLOC(d1); copyup(n1,w1);
! 572: w2 = W_UPALLOC(d2); copyup(n2,w2);
! 573: }
! 574: d1 = w1->d; d2 = w2->d;
! 575: while ( 1 ) {
! 576: remup_destructive(w1,w2);
! 577: if ( w1->d < 0 ) {
! 578: *nr = r = UPALLOC(w2->d); copyup(w2,r); return;
! 579: } else if ( !w1->d ) {
! 580: *nr = r = UPALLOC(0); r->d = 0;
! 581: divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]); return;
! 582: } else {
! 583: t = w1; w1 = w2; w2 = t;
! 584: }
! 585: }
! 586: }
! 587: }
! 588:
! 589: /* compute r s.t. a*r = 1 mod m */
! 590:
! 591: void extended_gcdup(a,m,r)
! 592: UP a,m;
! 593: UP *r;
! 594: {
! 595: UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;
! 596: Num i;
! 597:
! 598: one = UPALLOC(0); one->d = 0; one->c[0] = (Num)ONELM;
! 599: g1 = m; g2 = a;
! 600: a1 = one; a2 = 0;
! 601: b1 = 0; b2 = one;
! 602: /* a2*m+b2*a = g2 always holds */
! 603: while ( g2->d ) {
! 604: qrup(g1,g2,&quo,&rem); g1 = g2; g2 = rem;
! 605: mulup(a2,quo,&t); subup(a1,t,&a3); a1 = a2; a2 = a3;
! 606: mulup(b2,quo,&t); subup(b1,t,&b3); b1 = b2; b2 = b3;
! 607: }
! 608: /* now a2*m+b2*a = g2 (const) */
! 609: divnum(0,(Num)ONE,g2->c[0],&i);
! 610: inv = UPALLOC(0); inv->d = 0; inv->c[0] = i;
! 611: mulup(b2,inv,r);
! 612: }
! 613:
! 614: void reverseup(n1,d,nr)
! 615: UP n1;
! 616: int d;
! 617: UP *nr;
! 618: {
! 619: Num *c1,*c;
! 620: int i,d1;
! 621: UP r;
! 622:
! 623: if ( !n1 )
! 624: *nr = 0;
! 625: else if ( d < (d1 = n1->d) )
! 626: error("reverseup : invalid input");
! 627: else {
! 628: *nr = r = UPALLOC(d);
! 629: c = r->c;
! 630: bzero((char *)c,(d+1)*sizeof(Num));
! 631: c1 = n1->c;
! 632: for ( i = 0; i <= d1; i++ )
! 633: c[d-i] = c1[i];
! 634: for ( i = d; i >= 0 && !c[i]; i-- );
! 635: r->d = i;
! 636: if ( i < 0 )
! 637: *nr = 0;
! 638: }
! 639: }
! 640:
! 641: void invmodup(n1,d,nr)
! 642: UP n1;
! 643: int d;
! 644: UP *nr;
! 645: {
! 646: UP r;
! 647: Num s,t,u,hinv;
! 648: int i,j,dmin;
! 649: Num *w,*c,*cr;
! 650:
! 651: if ( !n1 || !n1->c[0] )
! 652: error("invmodup : invalid input");
! 653: else if ( !n1->d ) {
! 654: *nr = r = UPALLOC(0); r->d = 0;
! 655: divnum(0,(Num)ONE,n1->c[0],&r->c[0]);
! 656: } else {
! 657: *nr = r = UPALLOC(d);
! 658: w = (Num *)ALLOCA((d+1)*sizeof(Q));
! 659: bzero((char *)w,(d+1)*sizeof(Q));
! 660: dmin = MIN(d,n1->d);
! 661: divnum(0,(Num)ONE,n1->c[0],&hinv);
! 662: for ( i = 0, c = n1->c; i <= dmin; i++ )
! 663: mulnum(0,c[i],hinv,&w[i]);
! 664: cr = r->c;
! 665: cr[0] = w[0];
! 666: for ( i = 1; i <= d; i++ ) {
! 667: for ( j = 1, s = w[i]; j < i; j++ ) {
! 668: mulnum(0,cr[j],w[i-j],&t);
! 669: addnum(0,s,t,&u);
! 670: s = u;
! 671: }
! 672: chsgnnum(s,&cr[i]);
! 673: }
! 674: for ( i = 0; i <= d; i++ ) {
! 675: mulnum(0,cr[i],hinv,&t); cr[i] = t;
! 676: }
! 677: for ( i = d; i >= 0 && !cr[i]; i-- );
! 678: r->d = i;
! 679: }
! 680: }
! 681:
! 682: void pwrup(n,e,nr)
! 683: UP n;
! 684: Q e;
! 685: UP *nr;
! 686: {
! 687: UP y,z,t;
! 688: N en,qn;
! 689: int r;
! 690:
! 691: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)ONE;
! 692: z = n;
! 693: if ( !e )
! 694: *nr = y;
! 695: else if ( UNIQ(e) )
! 696: *nr = n;
! 697: else {
! 698: en = NM(e);
! 699: while ( 1 ) {
! 700: r = divin(en,2,&qn); en = qn;
! 701: if ( r ) {
! 702: mulup(z,y,&t); y = t;
! 703: if ( !en ) {
! 704: *nr = y;
! 705: return;
! 706: }
! 707: }
! 708: mulup(z,z,&t); z = t;
! 709: }
! 710: }
! 711: }
! 712:
! 713: int compup(n1,n2)
! 714: UP n1,n2;
! 715: {
! 716: int i,r;
! 717:
! 718: if ( !n1 )
! 719: return n2 ? -1 : 0;
! 720: else if ( !n2 )
! 721: return 1;
! 722: else if ( n1->d > n2->d )
! 723: return 1;
! 724: else if ( n1->d < n2->d )
! 725: return -1;
! 726: else {
! 727: for ( i = n1->d; i >= 0; i-- ) {
! 728: r = compnum(CO,n1->c[i],n2->c[i]);
! 729: if ( r )
! 730: return r;
! 731: }
! 732: return 0;
! 733: }
! 734: }
! 735:
! 736: void kmulp(vl,n1,n2,nr)
! 737: VL vl;
! 738: P n1,n2;
! 739: P *nr;
! 740: {
! 741: UP b1,b2,br;
! 742:
! 743: if ( !n1 || !n2 )
! 744: *nr = 0;
! 745: else if ( OID(n1) == O_N || OID(n2) == O_N )
! 746: mulp(vl,n1,n2,nr);
! 747: else {
! 748: ptoup(n1,&b1); ptoup(n2,&b2);
! 749: kmulup(b1,b2,&br);
! 750: uptop(br,nr);
! 751: }
! 752: }
! 753:
! 754: void ksquarep(vl,n1,nr)
! 755: VL vl;
! 756: P n1;
! 757: P *nr;
! 758: {
! 759: UP b1,br;
! 760:
! 761: if ( !n1 )
! 762: *nr = 0;
! 763: else if ( OID(n1) == O_N )
! 764: mulp(vl,n1,n1,nr);
! 765: else {
! 766: ptoup(n1,&b1);
! 767: ksquareup(b1,&br);
! 768: uptop(br,nr);
! 769: }
! 770: }
! 771:
! 772: void kmulup(n1,n2,nr)
! 773: UP n1,n2,*nr;
! 774: {
! 775: UP n,t,s,m,carry;
! 776: int d,d1,d2,len,i,l;
! 777: pointer *r,*r0;
! 778:
! 779: if ( !n1 || !n2 ) {
! 780: *nr = 0; return;
! 781: }
! 782: d1 = DEG(n1)+1; d2 = DEG(n2)+1;
! 783: if ( (d1 < up_kara_mag) && (d2 < up_kara_mag) )
! 784: mulup(n1,n2,nr);
! 785: else
! 786: kmulupmain(n1,n2,nr);
! 787: }
! 788:
! 789: void ksquareup(n1,nr)
! 790: UP n1,*nr;
! 791: {
! 792: int d1;
! 793: extern Q TWO;
! 794:
! 795: if ( !n1 ) {
! 796: *nr = 0; return;
! 797: }
! 798: d1 = DEG(n1)+1;
! 799: if ( (d1 < up_kara_mag) )
! 800: squareup(n1,nr);
! 801: else
! 802: ksquareupmain(n1,nr);
! 803: }
! 804:
! 805: void copyup(n1,n2)
! 806: UP n1,n2;
! 807: {
! 808: n2->d = n1->d;
! 809: bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));
! 810: }
! 811:
! 812: void c_copyup(n,len,p)
! 813: UP n;
! 814: int len;
! 815: pointer *p;
! 816: {
! 817: if ( n )
! 818: bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));
! 819: }
! 820:
! 821: void kmulupmain(n1,n2,nr)
! 822: UP n1,n2,*nr;
! 823: {
! 824: int d1,d2,h,len;
! 825: UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
! 826:
! 827: d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;
! 828: decompup(n1,h,&n1lo,&n1hi);
! 829: decompup(n2,h,&n2lo,&n2hi);
! 830: kmulup(n1lo,n2lo,&lo);
! 831: kmulup(n1hi,n2hi,&hi);
! 832: shiftup(hi,2*h,&t2);
! 833: addup(lo,t2,&t1);
! 834: addup(hi,lo,&mid1);
! 835: subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2); kmulup(s1,s2,&mid2);
! 836: addup(mid1,mid2,&mid);
! 837: shiftup(mid,h,&t2);
! 838: addup(t1,t2,nr);
! 839: }
! 840:
! 841: void ksquareupmain(n1,nr)
! 842: UP n1,*nr;
! 843: {
! 844: int d1,h,len;
! 845: UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
! 846:
! 847: d1 = DEG(n1)+1; h = (d1+1)/2;
! 848: decompup(n1,h,&n1lo,&n1hi);
! 849: ksquareup(n1hi,&hi); ksquareup(n1lo,&lo);
! 850: shiftup(hi,2*h,&t2);
! 851: addup(lo,t2,&t1);
! 852: addup(hi,lo,&mid1);
! 853: subup(n1hi,n1lo,&s1); ksquareup(s1,&mid2);
! 854: subup(mid1,mid2,&mid);
! 855: shiftup(mid,h,&t2);
! 856: addup(t1,t2,nr);
! 857: }
! 858:
! 859: void rembymulup(n1,n2,nr)
! 860: UP n1,n2;
! 861: UP *nr;
! 862: {
! 863: int d1,d2,d;
! 864: UP r1,r2,inv2,t,s,q;
! 865:
! 866: if ( !n2 )
! 867: error("rembymulup : division by 0");
! 868: else if ( !n1 || !n2->d )
! 869: *nr = 0;
! 870: else if ( (d1 = n1->d) < (d2 = n2->d) )
! 871: *nr = n1;
! 872: else {
! 873: d = d1-d2;
! 874: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
! 875: reverseup(n2,d2,&r2);
! 876: invmodup(r2,d,&inv2);
! 877: kmulup(r1,inv2,&t);
! 878: truncup(t,d+1,&s);
! 879: reverseup(s,d,&q);
! 880: kmulup(q,n2,&t); subup(n1,t,nr);
! 881: }
! 882: }
! 883:
! 884: /*
! 885: deg n2 = d
! 886: deg n1 <= 2*d-1
! 887: inv2 = inverse of reversep(n2) mod x^d
! 888: */
! 889:
! 890: void hybrid_rembymulup_special(ff,n1,n2,inv2,nr)
! 891: int ff;
! 892: UP n1,n2,inv2;
! 893: UP *nr;
! 894: {
! 895: int d1,d2,d;
! 896: UP r1,t,s,q,u;
! 897:
! 898: if ( !n2 )
! 899: error("hybrid_rembymulup : division by 0");
! 900: else if ( !n1 || !n2->d )
! 901: *nr = 0;
! 902: else if ( (d1 = n1->d) < (d2 = n2->d) )
! 903: *nr = n1;
! 904: else {
! 905: d = d1-d2;
! 906: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
! 907: hybrid_tmulup(ff,r1,inv2,d+1,&t);
! 908: truncup(t,d+1,&s);
! 909: reverseup(s,d,&q);
! 910:
! 911: hybrid_tmulup(ff,q,n2,d2,&t);
! 912: truncup(n1,d2,&s);
! 913: subup(s,t,nr);
! 914: }
! 915: }
! 916:
! 917: void rembymulup_special(n1,n2,inv2,nr)
! 918: UP n1,n2,inv2;
! 919: UP *nr;
! 920: {
! 921: int d1,d2,d;
! 922: UP r1,t,s,q,u;
! 923:
! 924: if ( !n2 )
! 925: error("rembymulup : division by 0");
! 926: else if ( !n1 || !n2->d )
! 927: *nr = 0;
! 928: else if ( (d1 = n1->d) < (d2 = n2->d) )
! 929: *nr = n1;
! 930: else {
! 931: d = d1-d2;
! 932: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
! 933: tkmulup(r1,inv2,d+1,&t);
! 934: truncup(t,d+1,&s);
! 935: reverseup(s,d,&q);
! 936:
! 937: tkmulup(q,n2,d2,&t);
! 938: truncup(n1,d2,&s);
! 939: subup(s,t,nr);
! 940: }
! 941: }
! 942:
! 943: /* *nr = n1*n2 mod x^d */
! 944:
! 945: void tkmulup(n1,n2,d,nr)
! 946: UP n1,n2;
! 947: int d;
! 948: UP *nr;
! 949: {
! 950: int m;
! 951: UP n1l,n1h,n2l,n2h,l,h,t,s,u,afo;
! 952:
! 953: if ( d < 0 )
! 954: error("tkmulup : invalid argument");
! 955: else if ( d == 0 )
! 956: *nr = 0;
! 957: else {
! 958: truncup(n1,d,&t); n1 = t;
! 959: truncup(n2,d,&t); n2 = t;
! 960: if ( !n1 || !n2 )
! 961: *nr = 0;
! 962: else if ( n1->d < up_tkara_mag || n2->d < up_tkara_mag )
! 963: tmulup(n1,n2,d,nr);
! 964: else {
! 965: m = (d+1)/2;
! 966: decompup(n1,m,&n1l,&n1h);
! 967: decompup(n2,m,&n2l,&n2h);
! 968: kmulup(n1l,n2l,&l);
! 969: tkmulup(n1l,n2h,d-m,&t);
! 970: tkmulup(n2l,n1h,d-m,&s);
! 971: addup(t,s,&u);
! 972: shiftup(u,m,&h);
! 973: addup(l,h,&t);
! 974: truncup(t,d,nr);
! 975: }
! 976: }
! 977: }
! 978:
! 979: /* n->n*x^d */
! 980:
! 981: void shiftup(n,d,nr)
! 982: UP n;
! 983: int d;
! 984: UP *nr;
! 985: {
! 986: int dr;
! 987: UP r;
! 988:
! 989: if ( !n )
! 990: *nr = 0;
! 991: else {
! 992: dr = n->d+d;
! 993: *nr = r = UPALLOC(dr); r->d = dr;
! 994: bzero(r->c,(dr+1)*sizeof(Num));
! 995: bcopy(n->c,r->c+d,(n->d+1)*sizeof(Num));
! 996: }
! 997: }
! 998:
! 999: void fft_rembymulup_special(n1,n2,inv2,nr)
! 1000: UP n1,n2,inv2;
! 1001: UP *nr;
! 1002: {
! 1003: int d1,d2,d;
! 1004: UP r1,t,s,q,u;
! 1005:
! 1006: if ( !n2 )
! 1007: error("rembymulup : division by 0");
! 1008: else if ( !n1 || !n2->d )
! 1009: *nr = 0;
! 1010: else if ( (d1 = n1->d) < (d2 = n2->d) )
! 1011: *nr = n1;
! 1012: else {
! 1013: d = d1-d2;
! 1014: reverseup(n1,d1,&t); truncup(t,d+1,&r1);
! 1015: trunc_fft_mulup(r1,inv2,d+1,&t);
! 1016: truncup(t,d+1,&s);
! 1017: reverseup(s,d,&q);
! 1018: trunc_fft_mulup(q,n2,d2,&t); truncup(t,d2,&u);
! 1019: truncup(n1,d2,&s);
! 1020: subup(s,u,nr);
! 1021: }
! 1022: }
! 1023:
! 1024: void set_degreeup(n,d)
! 1025: UP n;
! 1026: int d;
! 1027: {
! 1028: int i;
! 1029:
! 1030: for ( i = d; i >= 0 && !n->c[i]; i-- );
! 1031: n->d = i;
! 1032: }
! 1033:
! 1034: /* n -> n0 + x^d n1 */
! 1035:
! 1036: void decompup(n,d,n0,n1)
! 1037: UP n;
! 1038: int d;
! 1039: UP *n0,*n1;
! 1040: {
! 1041: int dn;
! 1042: UP r0,r1;
! 1043:
! 1044: if ( !n || d > n->d ) {
! 1045: *n0 = n; *n1 = 0;
! 1046: } else if ( d < 0 )
! 1047: error("decompup : invalid argument");
! 1048: else if ( d == 0 ) {
! 1049: *n0 = 0; *n1 = n;
! 1050: } else {
! 1051: dn = n->d;
! 1052: *n0 = r0 = UPALLOC(d-1);
! 1053: *n1 = r1 = UPALLOC(dn-d);
! 1054: bcopy(n->c,r0->c,d*sizeof(Num));
! 1055: bcopy(n->c+d,r1->c,(dn-d+1)*sizeof(Num));
! 1056: set_degreeup(r0,d-1);
! 1057: if ( r0->d < 0 )
! 1058: *n0 = 0;
! 1059: set_degreeup(r1,dn-d);
! 1060: if ( r1->d < 0 )
! 1061: *n1 = 0;
! 1062: }
! 1063: }
! 1064:
! 1065: /* n -> n mod x^d */
! 1066:
! 1067: void truncup(n1,d,nr)
! 1068: UP n1;
! 1069: int d;
! 1070: UP *nr;
! 1071: {
! 1072: int i;
! 1073: UP r;
! 1074:
! 1075: if ( !n1 || d > n1->d )
! 1076: *nr = n1;
! 1077: else if ( d < 0 )
! 1078: error("truncup : invalid argument");
! 1079: else if ( d == 0 )
! 1080: *nr = 0;
! 1081: else {
! 1082: *nr = r = UPALLOC(d-1);
! 1083: bcopy(n1->c,r->c,d*sizeof(Num));
! 1084: for ( i = d-1; i >= 0 && !r->c[i]; i-- );
! 1085: if ( i < 0 )
! 1086: *nr = 0;
! 1087: else
! 1088: r->d = i;
! 1089: }
! 1090: }
! 1091:
! 1092: int int_bits(t)
! 1093: int t;
! 1094: {
! 1095: int k;
! 1096:
! 1097: for ( k = 0; t; t>>=1, k++);
! 1098: return k;
! 1099: }
! 1100:
! 1101: /* n is assumed to be LM or integer coefficient */
! 1102:
! 1103: int maxblenup(n)
! 1104: UP n;
! 1105: {
! 1106: int m,r,i,d;
! 1107: Num *c;
! 1108:
! 1109: if ( !n )
! 1110: return 0;
! 1111: d = n->d; c = (Num *)n->c;
! 1112: for ( r = 0, i = 0; i <= d; i++ )
! 1113: if ( c[i] ) {
! 1114: switch ( NID(c[i]) ) {
! 1115: case N_LM:
! 1116: m = n_bits(((LM)c[i])->body);
! 1117: break;
! 1118: case N_Q:
! 1119: m = n_bits(((Q)c[i])->nm);
! 1120: break;
! 1121: default:
! 1122: error("maxblenup : invalid coefficient");
! 1123: }
! 1124: r = MAX(m,r);
! 1125: }
! 1126: return r;
! 1127: }
! 1128:
! 1129: void uptofmarray(mod,n,f)
! 1130: int mod;
! 1131: UP n;
! 1132: ModNum *f;
! 1133: {
! 1134: int d,i;
! 1135: unsigned int r;
! 1136: Num *c;
! 1137:
! 1138: if ( !n )
! 1139: return;
! 1140: else {
! 1141: d = n->d; c = (Num *)n->c;
! 1142: for ( i = 0; i <= d; i++ ) {
! 1143: if ( c[i] ) {
! 1144: switch ( NID(c[i]) ) {
! 1145: case N_LM:
! 1146: f[i] = (ModNum)rem(((LM)c[i])->body,mod);
! 1147: break;
! 1148: case N_Q:
! 1149: r = rem(NM((Q)c[i]),mod);
! 1150: if ( r && (SGN((Q)c[i])<0) )
! 1151: r = mod-r;
! 1152: f[i] = (ModNum)r;
! 1153: break;
! 1154: default:
! 1155: error("uptofmarray : invalid coefficient");
! 1156: }
! 1157: } else
! 1158: f[i] = 0;
! 1159: }
! 1160: }
! 1161: }
! 1162:
! 1163: void fmarraytoup(f,d,nr)
! 1164: ModNum *f;
! 1165: int d;
! 1166: UP *nr;
! 1167: {
! 1168: int i;
! 1169: Q *c;
! 1170: UP r;
! 1171:
! 1172: if ( d < 0 ) {
! 1173: *nr = 0;
! 1174: } else {
! 1175: *nr = r = UPALLOC(d); c = (Q *)r->c;
! 1176: for ( i = 0; i <= d; i++ ) {
! 1177: UTOQ((unsigned int)f[i],c[i]);
! 1178: }
! 1179: for ( i = d; i >= 0 && !c[i]; i-- );
! 1180: if ( i < 0 )
! 1181: *nr = 0;
! 1182: else
! 1183: r->d = i;
! 1184: }
! 1185: }
! 1186:
! 1187: /* f[i]: an array of length n */
! 1188:
! 1189: void uiarraytoup(f,n,d,nr)
! 1190: unsigned int **f;
! 1191: int n,d;
! 1192: UP *nr;
! 1193: {
! 1194: int i,j;
! 1195: unsigned int *fi;
! 1196: N ci;
! 1197: Q *c;
! 1198: UP r;
! 1199:
! 1200: if ( d < 0 ) {
! 1201: *nr = 0;
! 1202: } else {
! 1203: *nr = r = UPALLOC(d); c = (Q *)r->c;
! 1204: for ( i = 0; i <= d; i++ ) {
! 1205: fi = f[i];
! 1206: for ( j = n-1; j >= 0 && !fi[j]; j-- );
! 1207: j++;
! 1208: if ( j ) {
! 1209: ci = NALLOC(j);
! 1210: PL(ci) = j;
! 1211: bcopy(fi,BD(ci),j*sizeof(unsigned int));
! 1212: NTOQ(ci,1,c[i]);
! 1213: } else
! 1214: c[i] = 0;
! 1215: }
! 1216: for ( i = d; i >= 0 && !c[i]; i-- );
! 1217: if ( i < 0 )
! 1218: *nr = 0;
! 1219: else
! 1220: r->d = i;
! 1221: }
! 1222: }
! 1223:
! 1224: void adj_coefup(n,m,m2,nr)
! 1225: UP n;
! 1226: N m,m2;
! 1227: UP *nr;
! 1228: {
! 1229: int d;
! 1230: Q *c,*cr;
! 1231: Q mq;
! 1232: int i;
! 1233: UP r;
! 1234:
! 1235: if ( !n )
! 1236: *nr = 0;
! 1237: else {
! 1238: d = n->d; c = (Q *)n->c;
! 1239: *nr = r = UPALLOC(d); cr = (Q *)r->c;
! 1240: NTOQ(m,1,mq);
! 1241: for ( i = 0; i <= d; i++ ) {
! 1242: if ( !c[i] )
! 1243: continue;
! 1244: if ( cmpn(NM(c[i]),m2) > 0 )
! 1245: subq(c[i],mq,&cr[i]);
! 1246: else
! 1247: cr[i] = c[i];
! 1248: }
! 1249: for ( i = d; i >= 0 && !cr[i]; i-- );
! 1250: if ( i < 0 )
! 1251: *nr = 0;
! 1252: else
! 1253: r->d = i;
! 1254: }
! 1255: }
! 1256:
! 1257: /* n is assumed to have positive integer coefficients. */
! 1258:
! 1259: void remcup(n,mod,nr)
! 1260: UP n;
! 1261: N mod;
! 1262: UP *nr;
! 1263: {
! 1264: int i,d;
! 1265: Q *c,*cr;
! 1266: UP r;
! 1267: N t;
! 1268:
! 1269: if ( !n )
! 1270: *nr = 0;
! 1271: else {
! 1272: d = n->d; c = (Q *)n->c;
! 1273: *nr = r = UPALLOC(d); cr = (Q *)r->c;
! 1274: for ( i = 0; i <= d; i++ )
! 1275: if ( c[i] ) {
! 1276: remn(NM(c[i]),mod,&t);
! 1277: if ( t )
! 1278: NTOQ(t,1,cr[i]);
! 1279: else
! 1280: cr[i] = 0;
! 1281: } else
! 1282: cr[i] = 0;
! 1283: for ( i = d; i >= 0 && !cr[i]; i-- );
! 1284: if ( i < 0 )
! 1285: *nr = 0;
! 1286: else
! 1287: r->d = i;
! 1288: }
! 1289: }
! 1290:
! 1291: void fft_mulup_main(UP,UP,int,UP *);
! 1292:
! 1293: void fft_mulup(n1,n2,nr)
! 1294: UP n1,n2;
! 1295: UP *nr;
! 1296: {
! 1297: int d1,d2,d,b1,b2,h;
! 1298: UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;
! 1299:
! 1300: if ( !n1 || !n2 )
! 1301: *nr = 0;
! 1302: else {
! 1303: d1 = n1->d; d2 = n2->d;
! 1304: if ( MAX(d1,d2) < up_fft_mag )
! 1305: kmulup(n1,n2,nr);
! 1306: else {
! 1307: b1 = maxblenup(n1); b2 = maxblenup(n2);
! 1308: if ( fft_available(d1,b1,d2,b2) )
! 1309: fft_mulup_main(n1,n2,0,nr);
! 1310: else {
! 1311: d = MAX(d1,d2)+1; h = (d+1)/2;
! 1312: decompup(n1,h,&n1lo,&n1hi);
! 1313: decompup(n2,h,&n2lo,&n2hi);
! 1314: fft_mulup(n1lo,n2lo,&lo);
! 1315: fft_mulup(n1hi,n2hi,&hi);
! 1316: shiftup(hi,2*h,&t2);
! 1317: addup(lo,t2,&t1);
! 1318: addup(hi,lo,&mid1);
! 1319: subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2);
! 1320: fft_mulup(s1,s2,&mid2);
! 1321: addup(mid1,mid2,&mid);
! 1322: shiftup(mid,h,&t2);
! 1323: addup(t1,t2,nr);
! 1324: }
! 1325: }
! 1326: }
! 1327: }
! 1328:
! 1329: void trunc_fft_mulup(n1,n2,dbd,nr)
! 1330: UP n1,n2;
! 1331: int dbd;
! 1332: UP *nr;
! 1333: {
! 1334: int d1,d2,b1,b2,m;
! 1335: UP n1l,n1h,n2l,n2h,l,h,t,s,u;
! 1336:
! 1337: if ( !n1 || !n2 )
! 1338: *nr = 0;
! 1339: else if ( dbd < 0 )
! 1340: error("trunc_fft_mulup : invalid argument");
! 1341: else if ( dbd == 0 )
! 1342: *nr = 0;
! 1343: else {
! 1344: truncup(n1,dbd,&t); n1 = t;
! 1345: truncup(n2,dbd,&t); n2 = t;
! 1346: d1 = n1->d; d2 = n2->d;
! 1347: if ( MAX(d1,d2) < up_fft_mag )
! 1348: tkmulup(n1,n2,dbd,nr);
! 1349: else {
! 1350: b1 = maxblenup(n1); b2 = maxblenup(n2);
! 1351: if ( fft_available(d1,b1,d2,b2) )
! 1352: fft_mulup_main(n1,n2,dbd,nr);
! 1353: else {
! 1354: m = (dbd+1)/2;
! 1355: decompup(n1,m,&n1l,&n1h);
! 1356: decompup(n2,m,&n2l,&n2h);
! 1357: fft_mulup(n1l,n2l,&l);
! 1358: trunc_fft_mulup(n1l,n2h,dbd-m,&t);
! 1359: trunc_fft_mulup(n2l,n1h,dbd-m,&s);
! 1360: addup(t,s,&u);
! 1361: shiftup(u,m,&h);
! 1362: addup(l,h,&t);
! 1363: truncup(t,dbd,nr);
! 1364: }
! 1365: }
! 1366: }
! 1367: }
! 1368:
! 1369: void fft_squareup(n1,nr)
! 1370: UP n1;
! 1371: UP *nr;
! 1372: {
! 1373: int d1,d,h,len,b1;
! 1374: UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
! 1375:
! 1376: if ( !n1 )
! 1377: *nr = 0;
! 1378: else {
! 1379: d1 = n1->d;
! 1380: if ( d1 < up_fft_mag )
! 1381: ksquareup(n1,nr);
! 1382: else {
! 1383: b1 = maxblenup(n1);
! 1384: if ( fft_available(d1,b1,d1,b1) )
! 1385: fft_mulup_main(n1,n1,0,nr);
! 1386: else {
! 1387: d = d1+1; h = (d1+1)/2;
! 1388: decompup(n1,h,&n1lo,&n1hi);
! 1389: fft_squareup(n1hi,&hi);
! 1390: fft_squareup(n1lo,&lo);
! 1391: shiftup(hi,2*h,&t2);
! 1392: addup(lo,t2,&t1);
! 1393: addup(hi,lo,&mid1);
! 1394: subup(n1hi,n1lo,&s1);
! 1395: fft_squareup(s1,&mid2);
! 1396: subup(mid1,mid2,&mid);
! 1397: shiftup(mid,h,&t2);
! 1398: addup(t1,t2,nr);
! 1399: }
! 1400: }
! 1401: }
! 1402: }
! 1403:
! 1404: /*
! 1405: * dbd == 0 => n1 * n2
! 1406: * dbd > 0 => n1 * n2 mod x^dbd
! 1407: * n1 == n2 => squaring
! 1408: */
! 1409:
! 1410: void fft_mulup_main(n1,n2,dbd,nr)
! 1411: UP n1,n2;
! 1412: UP *nr;
! 1413: {
! 1414: ModNum *f1,*f2,*w,*fr;
! 1415: ModNum **frarray,**fa;
! 1416: unsigned int *modarray,*ma;
! 1417: int modarray_len;
! 1418: int frarray_index = 0;
! 1419: N m,m1,m2;
! 1420: int d1,d2,dmin,i,mod,root,d,cond,bound;
! 1421: UP r;
! 1422:
! 1423: if ( !n1 || !n2 ) {
! 1424: *nr = 0; return;
! 1425: }
! 1426: d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
! 1427: if ( !d1 || !d2 ) {
! 1428: mulup(n1,n2,nr); return;
! 1429: }
! 1430: m = ONEN;
! 1431: bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
! 1432: f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
! 1433: if ( n1 == n2 )
! 1434: f2 = 0;
! 1435: else
! 1436: f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
! 1437: w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
! 1438: modarray_len = 1024;
! 1439: modarray = (unsigned int *)MALLOC_ATOMIC(modarray_len*sizeof(unsigned int));
! 1440: frarray = (ModNum **)MALLOC(modarray_len*sizeof(ModNum *));
! 1441: for ( i = 0; i < NPrimes; i++ ) {
! 1442: FFT_primes(i,&mod,&root,&d);
! 1443: if ( (1<<d) < d1+d2+1 )
! 1444: continue;
! 1445: if ( frarray_index == modarray_len ) {
! 1446: ma = (unsigned int *)MALLOC_ATOMIC(2*modarray_len*sizeof(unsigned int));
! 1447: bcopy(modarray,ma,modarray_len*sizeof(unsigned int));
! 1448: modarray = ma;
! 1449: fa = (ModNum **)MALLOC(2*modarray_len*sizeof(ModNum *));
! 1450: bcopy(frarray,fa,modarray_len*sizeof(ModNum *));
! 1451: frarray = fa;
! 1452: modarray_len *= 2;
! 1453: }
! 1454: modarray[frarray_index] = mod;
! 1455: frarray[frarray_index++] = fr
! 1456: = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
! 1457: uptofmarray(mod,n1,f1);
! 1458: if ( !f2 )
! 1459: cond = FFT_pol_square(d1,f1,fr,i,w);
! 1460: else {
! 1461: uptofmarray(mod,n2,f2);
! 1462: cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
! 1463: }
! 1464: if ( cond )
! 1465: error("fft_mulup : error in FFT_pol_product");
! 1466: STON(mod,m1); muln(m,m1,&m2); m = m2;
! 1467: if ( n_bits(m) > bound ) {
! 1468: if ( !dbd )
! 1469: dbd = d1+d2+1;
! 1470: crup(frarray,MIN(d1+d2,dbd-1),modarray,frarray_index,m,&r);
! 1471: divin(m,2,&m2);
! 1472: adj_coefup(r,m,m2,nr);
! 1473: return;
! 1474: }
! 1475: }
! 1476: error("fft_mulup : FFT_primes exhausted");
! 1477: }
! 1478: #if 0
! 1479: /* inefficient version */
! 1480:
! 1481: void crup(f,d,mod,index,m,r)
! 1482: ModNum **f;
! 1483: int d;
! 1484: int *mod;
! 1485: int index;
! 1486: N m;
! 1487: UP *r;
! 1488: {
! 1489: N *cof,*c;
! 1490: int *inv;
! 1491: struct oUP w;
! 1492: int i;
! 1493: UP s,s1,s2;
! 1494: N t;
! 1495: UP *g;
! 1496: Q q;
! 1497: struct oEGT eg0,eg1;
! 1498:
! 1499: get_eg(&eg0);
! 1500: w.d = 0;
! 1501: inv = (int *)ALLOCA(index*sizeof(int));
! 1502: cof = (N *)ALLOCA(index*sizeof(N));
! 1503: c = (N *)ALLOCA(index*sizeof(N));
! 1504: g = (UP *)ALLOCA(index*sizeof(UP));
! 1505: up_lazy = 1;
! 1506: for ( i = 0, s = 0; i < index; i++ ) {
! 1507: divin(m,mod[i],&cof[i]);
! 1508: inv[i] = invm(rem(cof[i],mod[i]),mod[i]);
! 1509: STON(inv[i],t);
! 1510: muln(cof[i],t,&c[i]);
! 1511: NTOQ(c[i],1,q); w.c[0] = (Num)q;
! 1512: fmarraytoup(f[i],d,&g[i]);
! 1513: mulup(g[i],&w,&s1);
! 1514: addup(s,s1,&s2);
! 1515: s = s2;
! 1516: }
! 1517: up_lazy = 0;
! 1518: remcup(s,m,r);
! 1519: get_eg(&eg1);
! 1520: add_eg(&eg_chrem,&eg0,&eg1);
! 1521: }
! 1522:
! 1523: #else
! 1524: /* improved version */
! 1525:
! 1526: void crup(f,d,mod,index,m,r)
! 1527: ModNum **f;
! 1528: int d;
! 1529: int *mod;
! 1530: int index;
! 1531: N m;
! 1532: UP *r;
! 1533: {
! 1534: N cof,c,t,w,w1;
! 1535: struct oN fc;
! 1536: N *s;
! 1537: UP u;
! 1538: Q q;
! 1539: int inv,i,j,rlen;
! 1540:
! 1541: rlen = PL(m)+10; /* XXX */
! 1542: PL(&fc) = 1;
! 1543: c = NALLOC(rlen);
! 1544: w = NALLOC(rlen);
! 1545: w1 = NALLOC(rlen);
! 1546: s = (N *)MALLOC((d+1)*sizeof(N));
! 1547: for ( i = 0; i <= d; i++ ) {
! 1548: s[i] = NALLOC(rlen);
! 1549: PL(s[i]) = 0;
! 1550: bzero(BD(s[i]),rlen*sizeof(unsigned int));
! 1551: }
! 1552: for ( i = 0; i < index; i++ ) {
! 1553: divin(m,mod[i],&cof);
! 1554: inv = invm(rem(cof,mod[i]),mod[i]); STON(inv,t);
! 1555: _muln(cof,t,c);
! 1556: /* s += c*f[i] */
! 1557: for ( j = 0; j <= d; j++ )
! 1558: if ( f[i][j] ) {
! 1559: BD(&fc)[0] = f[i][j];
! 1560: _muln(c,&fc,w);
! 1561: _addn(s[j],w,w1);
! 1562: dupn(w1,s[j]);
! 1563: }
! 1564: }
! 1565: for ( i = d; i >= 0; i-- )
! 1566: if ( s[i] )
! 1567: break;
! 1568: if ( i < 0 )
! 1569: *r = 0;
! 1570: else {
! 1571: u = UPALLOC(i);
! 1572: DEG(u) = i;
! 1573: for ( j = 0; j <= i; j++ ) {
! 1574: if ( PL(s[j]) )
! 1575: NTOQ(s[j],1,q);
! 1576: else
! 1577: q = 0;
! 1578: COEF(u)[j] = (Num)q;
! 1579: }
! 1580: remcup(u,m,r);
! 1581: }
! 1582: }
! 1583: #endif
! 1584: