Annotation of OpenXM_contrib2/asir2000/engine/Q.c, Revision 1.10
1.4 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.5 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.4 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.10 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.9 2002/08/14 04:49:52 noro Exp $
1.4 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52: #include "inline.h"
53:
1.7 noro 54: void addq(Q n1,Q n2,Q *nr)
1.1 noro 55: {
1.10 ! noro 56: N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
! 57: int sgn;
! 58: unsigned t,t1;
! 59:
! 60: if ( !n1 )
! 61: *nr = n2;
! 62: else if ( !n2 )
! 63: *nr = n1;
! 64: else if ( INT(n1) )
! 65: if ( INT(n2) ) {
! 66: nm1 = NM(n1); nm2 = NM(n2);
! 67: if ( SGN(n1) == SGN(n2) ) {
! 68: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 69: t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
! 70: if ( t < t1 ) {
! 71: m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
! 72: } else
! 73: UTON(t,m);
! 74: } else
! 75: addn(NM(n1),NM(n2),&m);
! 76: NTOQ(m,SGN(n1),*nr);
! 77: } else {
! 78: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 79: t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
! 80: if ( !t )
! 81: sgn = 0;
! 82: else if ( t > t1 ) {
! 83: sgn = -1; t = -((int)t); UTON(t,m);
! 84: } else {
! 85: sgn = 1; UTON(t,m);
! 86: }
! 87: } else
! 88: sgn = subn(NM(n1),NM(n2),&m);
! 89: if ( sgn ) {
! 90: if ( SGN(n1) == sgn )
! 91: NTOQ(m,1,*nr);
! 92: else
! 93: NTOQ(m,-1,*nr);
! 94: } else
! 95: *nr = 0;
! 96: }
! 97: } else {
! 98: kmuln(NM(n1),DN(n2),&m);
! 99: if ( SGN(n1) == SGN(n2) ) {
! 100: sgn = SGN(n1); addn(m,NM(n2),&nm);
! 101: } else
! 102: sgn = SGN(n1)*subn(m,NM(n2),&nm);
! 103: if ( sgn ) {
! 104: dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
! 105: } else
! 106: *nr = 0;
! 107: }
! 108: else if ( INT(n2) ) {
! 109: kmuln(NM(n2),DN(n1),&m);
! 110: if ( SGN(n1) == SGN(n2) ) {
! 111: sgn = SGN(n1); addn(m,NM(n1),&nm);
! 112: } else
! 113: sgn = SGN(n1)*subn(NM(n1),m,&nm);
! 114: if ( sgn ) {
! 115: dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
! 116: } else
! 117: *nr = 0;
! 118: } else {
! 119: gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
! 120: kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
! 121: if ( SGN(n1) == SGN(n2) ) {
! 122: sgn = SGN(n1); addn(nm1,nm2,&nm3);
! 123: } else
! 124: sgn = SGN(n1)*subn(nm1,nm2,&nm3);
! 125: if ( sgn ) {
! 126: gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
! 127: kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
! 128: if ( UNIN(dn) )
! 129: NTOQ(nm,sgn,*nr);
! 130: else
! 131: NDTOQ(nm,dn,sgn,*nr);
! 132: } else
! 133: *nr = 0;
! 134: }
1.1 noro 135: }
136:
1.7 noro 137: void subq(Q n1,Q n2,Q *nr)
1.1 noro 138: {
1.10 ! noro 139: Q m;
1.1 noro 140:
1.10 ! noro 141: if ( !n1 )
! 142: if ( !n2 )
! 143: *nr = 0;
! 144: else {
! 145: DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
! 146: }
! 147: else if ( !n2 )
! 148: *nr = n1;
! 149: else if ( n1 == n2 )
! 150: *nr = 0;
! 151: else {
! 152: DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
! 153: }
1.1 noro 154: }
155:
1.7 noro 156: void mulq(Q n1,Q n2,Q *nr)
1.1 noro 157: {
1.10 ! noro 158: N nm,nm1,nm2,dn,dn1,dn2,g;
! 159: int sgn;
! 160: unsigned u,l;
! 161:
! 162: if ( !n1 || !n2 ) *nr = 0;
! 163: else if ( INT(n1) )
! 164: if ( INT(n2) ) {
! 165: nm1 = NM(n1); nm2 = NM(n2);
! 166: if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
! 167: DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
! 168: if ( u ) {
! 169: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
! 170: } else {
! 171: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
! 172: }
! 173: } else
! 174: kmuln(nm1,nm2,&nm);
! 175: sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
! 176: } else {
! 177: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
! 178: kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
! 179: if ( UNIN(dn) )
! 180: NTOQ(nm,sgn,*nr);
! 181: else
! 182: NDTOQ(nm,dn,sgn,*nr);
! 183: }
! 184: else if ( INT(n2) ) {
! 185: gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
! 186: kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
! 187: if ( UNIN(dn) )
! 188: NTOQ(nm,sgn,*nr);
! 189: else
! 190: NDTOQ(nm,dn,sgn,*nr);
! 191: } else {
! 192: gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
! 193: gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
! 194: kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
! 195: if ( UNIN(dn) )
! 196: NTOQ(nm,sgn,*nr);
! 197: else
! 198: NDTOQ(nm,dn,sgn,*nr);
! 199: }
1.1 noro 200: }
201:
1.7 noro 202: void divq(Q n1,Q n2,Q *nq)
1.1 noro 203: {
1.10 ! noro 204: Q m;
1.1 noro 205:
1.10 ! noro 206: if ( !n2 ) {
! 207: error("division by 0");
! 208: *nq = 0;
! 209: return;
! 210: } else if ( !n1 )
! 211: *nq = 0;
! 212: else if ( n1 == n2 )
! 213: *nq = ONE;
! 214: else {
! 215: invq(n2,&m); mulq(n1,m,nq);
! 216: }
1.6 noro 217: }
218:
1.7 noro 219: void divsq(Q n1,Q n2,Q *nq)
1.6 noro 220: {
1.10 ! noro 221: int sgn;
! 222: N tn;
1.6 noro 223:
1.10 ! noro 224: if ( !n2 ) {
! 225: error("division by 0");
! 226: *nq = 0;
! 227: return;
! 228: } else if ( !n1 )
! 229: *nq = 0;
! 230: else if ( n1 == n2 )
! 231: *nq = ONE;
! 232: else {
! 233: divsn(NM(n1),NM(n2),&tn);
! 234: sgn = SGN(n1)*SGN(n2);
! 235: NTOQ(tn,sgn,*nq);
! 236: }
1.1 noro 237: }
238:
1.7 noro 239: void invq(Q n,Q *nr)
1.1 noro 240: {
1.10 ! noro 241: N nm,dn;
1.1 noro 242:
1.10 ! noro 243: if ( INT(n) )
! 244: if ( UNIN(NM(n)) )
! 245: if ( SGN(n) > 0 )
! 246: *nr = ONE;
! 247: else {
! 248: nm = ONEN; NTOQ(nm,SGN(n),*nr);
! 249: }
! 250: else {
! 251: nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
! 252: }
! 253: else if ( UNIN(NM(n)) ) {
! 254: nm = DN(n); NTOQ(nm,SGN(n),*nr);
! 255: } else {
! 256: nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
! 257: }
1.1 noro 258: }
259:
1.7 noro 260: void chsgnq(Q n,Q *nr)
1.1 noro 261: {
1.10 ! noro 262: Q t;
1.1 noro 263:
1.10 ! noro 264: if ( !n )
! 265: *nr = 0;
! 266: else {
! 267: DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
! 268: }
1.1 noro 269: }
270:
1.7 noro 271: void pwrq(Q n1,Q n,Q *nr)
1.1 noro 272: {
1.10 ! noro 273: N nm,dn;
! 274: int sgn;
1.1 noro 275:
1.10 ! noro 276: if ( !n || UNIQ(n1) )
! 277: *nr = ONE;
! 278: else if ( !n1 )
! 279: *nr = 0;
! 280: else if ( !INT(n) ) {
! 281: error("can't calculate fractional power."); *nr = 0;
! 282: } else if ( MUNIQ(n1) )
! 283: *nr = BD(NM(n))[0]%2 ? n1 : ONE;
! 284: else if ( PL(NM(n)) > 1 ) {
! 285: error("exponent too big."); *nr = 0;
! 286: } else {
! 287: sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
! 288: pwrn(NM(n1),BD(NM(n))[0],&nm);
! 289: if ( INT(n1) ) {
! 290: if ( UNIN(nm) )
! 291: if ( sgn == 1 )
! 292: *nr = ONE;
! 293: else
! 294: NTOQ(ONEN,sgn,*nr);
! 295: else if ( SGN(n) > 0 )
! 296: NTOQ(nm,sgn,*nr);
! 297: else
! 298: NDTOQ(ONEN,nm,sgn,*nr);
! 299: } else {
! 300: pwrn(DN(n1),BD(NM(n))[0],&dn);
! 301: if ( SGN(n) > 0 )
! 302: NDTOQ(nm,dn,sgn,*nr);
! 303: else if ( UNIN(nm) )
! 304: NTOQ(dn,sgn,*nr);
! 305: else
! 306: NDTOQ(dn,nm,sgn,*nr);
! 307: }
! 308: }
1.1 noro 309: }
310:
1.7 noro 311: int cmpq(Q q1,Q q2)
1.1 noro 312: {
1.10 ! noro 313: int sgn;
! 314: N t,s;
1.1 noro 315:
1.10 ! noro 316: if ( !q1 )
! 317: if ( !q2 )
! 318: return ( 0 );
! 319: else
! 320: return ( -SGN(q2) );
! 321: else if ( !q2 )
! 322: return ( SGN(q1) );
! 323: else if ( SGN(q1) != SGN(q2) )
! 324: return ( SGN(q1) );
! 325: else if ( INT(q1) )
! 326: if ( INT(q2) ) {
! 327: sgn = cmpn(NM(q1),NM(q2));
! 328: if ( !sgn )
! 329: return ( 0 );
! 330: else
! 331: return ( SGN(q1)==sgn?1:-1 );
! 332: } else {
! 333: kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
! 334: return ( SGN(q1) * sgn );
! 335: }
! 336: else if ( INT(q2) ) {
! 337: kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
! 338: return ( SGN(q1) * sgn );
! 339: } else {
! 340: kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
! 341: return ( SGN(q1) * sgn );
! 342: }
1.1 noro 343: }
344:
1.7 noro 345: void remq(Q n,Q m,Q *nr)
1.1 noro 346: {
1.10 ! noro 347: N q,r,s;
1.1 noro 348:
1.10 ! noro 349: if ( !n ) {
! 350: *nr = 0;
! 351: return;
! 352: }
! 353: divn(NM(n),NM(m),&q,&r);
! 354: if ( !r )
! 355: *nr = 0;
! 356: else if ( SGN(n) > 0 )
! 357: NTOQ(r,1,*nr);
! 358: else {
! 359: subn(NM(m),r,&s); NTOQ(s,1,*nr);
! 360: }
1.1 noro 361: }
362:
1.2 noro 363: /* t = [nC0 nC1 ... nCn] */
364:
1.7 noro 365: void mkbc(int n,Q *t)
1.1 noro 366: {
1.10 ! noro 367: int i;
! 368: N c,d;
1.1 noro 369:
1.10 ! noro 370: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
! 371: c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
! 372: kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
! 373: }
! 374: for ( ; i <= n; i++ )
! 375: t[i] = t[n-i];
1.2 noro 376: }
377:
378: /*
379: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
380: *
381: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
382: * where W(k,l,i) = i! * kCi * lCi
383: */
384:
1.7 noro 385: void mkwc(int k,int l,Q *t)
1.2 noro 386: {
1.10 ! noro 387: int i,n,up,low;
! 388: N nm,d,c;
1.2 noro 389:
1.10 ! noro 390: n = MIN(k,l);
! 391: for ( t[0] = ONE, i = 1; i <= n; i++ ) {
! 392: DM(k-i+1,l-i+1,up,low);
! 393: if ( up ) {
! 394: nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
! 395: } else {
! 396: nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
! 397: }
! 398: kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
! 399: }
1.3 noro 400: }
401:
402: /* mod m table */
403: /* XXX : should be optimized */
404:
1.7 noro 405: void mkwcm(int k,int l,int m,int *t)
1.3 noro 406: {
1.10 ! noro 407: int i,n;
! 408: Q *s;
1.3 noro 409:
1.10 ! noro 410: n = MIN(k,l);
! 411: s = (Q *)ALLOCA((n+1)*sizeof(Q));
! 412: mkwc(k,l,s);
! 413: for ( i = 0; i <= n; i++ ) {
! 414: t[i] = rem(NM(s[i]),m);
! 415: }
1.1 noro 416: }
417:
1.8 noro 418: #if 0
1.7 noro 419: void factorial(int n,Q *r)
1.1 noro 420: {
1.10 ! noro 421: Q t,iq,s;
! 422: unsigned int i,m,m0;
! 423: unsigned int c;
! 424:
! 425: if ( !n )
! 426: *r = ONE;
! 427: else if ( n < 0 )
! 428: *r = 0;
! 429: else {
! 430: for ( i = 1, t = ONE; (int)i <= n; ) {
! 431: for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
! 432: m0 = m; DM(m0,i,c,m)
! 433: }
! 434: if ( c ) {
! 435: m = m0; i--;
! 436: }
! 437: UTOQ(m,iq); mulq(t,iq,&s); t = s;
! 438: }
! 439: *r = t;
! 440: }
1.8 noro 441: }
442: #endif
443:
444: void partial_factorial(int s,int e,N *r);
445:
446: void factorial(int n,Q *r)
447: {
1.10 ! noro 448: N nm;
1.8 noro 449:
1.10 ! noro 450: if ( !n )
! 451: *r = ONE;
! 452: else if ( n < 0 )
! 453: *r = 0;
! 454: else {
! 455: partial_factorial(1,n,&nm);
! 456: NTOQ(nm,1,*r);
! 457: }
1.8 noro 458: }
459:
460: /* s*(s+1)*...*e */
461: void partial_factorial(int s,int e,N *r)
462: {
1.10 ! noro 463: unsigned int len,i,m,m0,c;
! 464: unsigned int *p,*p1,*p2;
! 465: N nm,r1,r2;
! 466:
! 467: /* XXX */
! 468: if ( e-s > 2000 ) {
! 469: m = (e+s)/2;
! 470: partial_factorial(s,m,&r1);
! 471: partial_factorial(m+1,e,&r2);
! 472: kmuln(r1,r2,r);
! 473: return;
! 474: }
! 475: /* find i s.t. 2^(i-1) < m <= 2^i */
! 476: for ( i = 0, m = 1; m < e; m <<=1, i++ );
! 477: /* XXX estimate of word length of the result */
! 478: len = (i*(e-s+1)+1+31)/32;
! 479: p = ALLOCA((len+1)*sizeof(int));
! 480: p1 = ALLOCA((len+1)*sizeof(int));
! 481: p[0] = s;
! 482: len = 1;
! 483: for ( i = s+1; (int)i <= e; ) {
! 484: for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) {
! 485: m0 = m; DM(m0,i,c,m)
! 486: }
! 487: if ( c ) {
! 488: m = m0; i--;
! 489: }
! 490: bzero(p1,(len+1)*sizeof(int));
! 491: muln_1(p,len,m,p1);
! 492: if ( p1[len] )
! 493: len++;
! 494: p2 = p; p = p1; p1 = p2;
! 495: }
! 496: *r = nm = NALLOC(len);
! 497: bcopy(p,BD(nm),sizeof(int)*len);
! 498: PL(nm) = len;
1.1 noro 499: }
500:
1.7 noro 501: void invl(Q a,Q mod,Q *ar)
1.1 noro 502: {
1.10 ! noro 503: Q f1,f2,a1,a2,q,m,s;
! 504: N qn,rn;
1.1 noro 505:
1.10 ! noro 506: a1 = ONE; a2 = 0;
! 507: DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
1.1 noro 508:
1.10 ! noro 509: while ( 1 ) {
! 510: divn(NM(f1),NM(f2),&qn,&rn);
! 511: if ( !qn )
! 512: q = 0;
! 513: else
! 514: NTOQ(qn,1,q);
! 515:
! 516: f1 = f2;
! 517: if ( !rn )
! 518: break;
! 519: else
! 520: NTOQ(rn,1,f2);
! 521:
! 522: mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
! 523: if ( !s )
! 524: a2 = 0;
! 525: else
! 526: remq(s,mod,&a2);
! 527: }
! 528: if ( SGN(a) < 0 )
! 529: chsgnq(a2,&m);
! 530: else
! 531: m = a2;
! 532:
! 533: if ( SGN(m) >= 0 )
! 534: *ar = m;
! 535: else
! 536: addq(m,mod,ar);
1.1 noro 537: }
538:
539: int kara_mag=100;
540:
1.7 noro 541: void kmuln(N n1,N n2,N *nr)
1.1 noro 542: {
1.10 ! noro 543: N n,t,s,m,carry;
! 544: int d,d1,d2,len,i,l;
! 545: int *r,*r0;
! 546:
! 547: if ( !n1 || !n2 ) {
! 548: *nr = 0; return;
! 549: }
! 550: d1 = PL(n1); d2 = PL(n2);
! 551: if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
! 552: muln(n1,n2,nr); return;
! 553: }
! 554: if ( d1 < d2 ) {
! 555: n = n1; n1 = n2; n2 = n;
! 556: d = d1; d1 = d2; d2 = d;
! 557: }
! 558: if ( d2 > (d1+1)/2 ) {
! 559: kmulnmain(n1,n2,nr); return;
! 560: }
! 561: d = (d1/d2)+((d1%d2)!=0);
! 562: len = (d+1)*d2;
! 563: r0 = (int *)ALLOCA(len*sizeof(int));
! 564: bzero((char *)r0,len*sizeof(int));
! 565: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
! 566: extractn(n1,i*d2,d2,&m);
! 567: if ( m ) {
! 568: kmuln(m,n2,&t);
! 569: addn(t,carry,&s);
! 570: copyn(s,d2,r);
! 571: extractn(s,d2,d2,&carry);
! 572: } else {
! 573: copyn(carry,d2,r);
! 574: carry = 0;
! 575: }
! 576: }
! 577: copyn(carry,d2,r);
! 578: for ( l = len - 1; !r0[l]; l-- );
! 579: l++;
! 580: *nr = t = NALLOC(l); PL(t) = l;
! 581: bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
1.1 noro 582: }
583:
1.7 noro 584: void extractn(N n,int index,int len,N *nr)
1.1 noro 585: {
1.10 ! noro 586: unsigned int *m;
! 587: int l;
! 588: N t;
! 589:
! 590: if ( !n ) {
! 591: *nr = 0; return;
! 592: }
! 593: m = BD(n)+index;
! 594: if ( (l = PL(n)-index) >= len ) {
! 595: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
! 596: l++;
! 597: }
! 598: if ( l <= 0 ) {
! 599: *nr = 0; return;
! 600: } else {
! 601: *nr = t = NALLOC(l); PL(t) = l;
! 602: bcopy((char *)m,(char *)BD(t),l*sizeof(int));
! 603: }
1.1 noro 604: }
605:
1.7 noro 606: void copyn(N n,int len,int *p)
1.1 noro 607: {
1.10 ! noro 608: if ( n )
! 609: bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
1.1 noro 610: }
611:
1.7 noro 612: void dupn(N n,N p)
1.1 noro 613: {
1.10 ! noro 614: if ( n )
! 615: bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
1.1 noro 616: }
617:
1.7 noro 618: void kmulnmain(N n1,N n2,N *nr)
1.1 noro 619: {
1.10 ! noro 620: int d1,d2,h,sgn,sgn1,sgn2,len;
! 621: N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
1.1 noro 622:
1.10 ! noro 623: d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
! 624: extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
! 625: extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
! 626: kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
! 627: len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
! 628: bzero((char *)BD(t1),len*sizeof(int));
! 629: if ( lo )
! 630: bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
! 631: if ( hi )
! 632: bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
! 633:
! 634: addn(hi,lo,&mid1);
! 635: sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
! 636: kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
! 637: if ( sgn > 0 )
! 638: addn(mid1,mid2,&mid);
! 639: else
! 640: subn(mid1,mid2,&mid);
! 641: if ( mid ) {
! 642: len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
! 643: bzero((char *)BD(t2),len*sizeof(int));
! 644: bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
! 645: addn(t1,t2,nr);
! 646: } else
! 647: *nr = t1;
1.1 noro 648: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Q.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine