=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Q.c,v retrieving revision 1.1.1.1 retrieving revision 1.10 diff -u -p -r1.1.1.1 -r1.10 --- OpenXM_contrib2/asir2000/engine/Q.c 1999/12/03 07:39:08 1.1.1.1 +++ OpenXM_contrib2/asir2000/engine/Q.c 2018/03/29 01:32:51 1.10 @@ -1,504 +1,648 @@ -/* $OpenXM: OpenXM/src/asir99/engine/Q.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, + * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. + * + * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.9 2002/08/14 04:49:52 noro Exp $ +*/ #include "ca.h" #include "base.h" #include "inline.h" -void addq(n1,n2,nr) -Q n1,n2,*nr; +void addq(Q n1,Q n2,Q *nr) { - N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m; - int sgn; - unsigned t,t1; + N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m; + int sgn; + unsigned t,t1; - if ( !n1 ) - *nr = n2; - else if ( !n2 ) - *nr = n1; - else if ( INT(n1) ) - if ( INT(n2) ) { - nm1 = NM(n1); nm2 = NM(n2); - if ( SGN(n1) == SGN(n2) ) { - if ( PL(nm1) == 1 && PL(nm2) == 1 ) { - t1 = BD(nm1)[0]; t = t1+BD(nm2)[0]; - if ( t < t1 ) { - m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1; - } else - UTON(t,m); - } else - addn(NM(n1),NM(n2),&m); - NTOQ(m,SGN(n1),*nr); - } else { - if ( PL(nm1) == 1 && PL(nm2) == 1 ) { - t1 = BD(nm1)[0]; t = t1-BD(nm2)[0]; - if ( !t ) - sgn = 0; - else if ( t > t1 ) { - sgn = -1; t = -((int)t); UTON(t,m); - } else { - sgn = 1; UTON(t,m); - } - } else - sgn = subn(NM(n1),NM(n2),&m); - if ( sgn ) { - if ( SGN(n1) == sgn ) - NTOQ(m,1,*nr); - else - NTOQ(m,-1,*nr); - } else - *nr = 0; - } - } else { - kmuln(NM(n1),DN(n2),&m); - if ( SGN(n1) == SGN(n2) ) { - sgn = SGN(n1); addn(m,NM(n2),&nm); - } else - sgn = SGN(n1)*subn(m,NM(n2),&nm); - if ( sgn ) { - dn = DN(n2); NDTOQ(nm,dn,sgn,*nr); - } else - *nr = 0; - } - else if ( INT(n2) ) { - kmuln(NM(n2),DN(n1),&m); - if ( SGN(n1) == SGN(n2) ) { - sgn = SGN(n1); addn(m,NM(n1),&nm); - } else - sgn = SGN(n1)*subn(NM(n1),m,&nm); - if ( sgn ) { - dn = DN(n1); NDTOQ(nm,dn,sgn,*nr); - } else - *nr = 0; - } else { - gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2); - kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2); - if ( SGN(n1) == SGN(n2) ) { - sgn = SGN(n1); addn(nm1,nm2,&nm3); - } else - sgn = SGN(n1)*subn(nm1,nm2,&nm3); - if ( sgn ) { - gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0); - kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn); - if ( UNIN(dn) ) - NTOQ(nm,sgn,*nr); - else - NDTOQ(nm,dn,sgn,*nr); - } else - *nr = 0; - } + if ( !n1 ) + *nr = n2; + else if ( !n2 ) + *nr = n1; + else if ( INT(n1) ) + if ( INT(n2) ) { + nm1 = NM(n1); nm2 = NM(n2); + if ( SGN(n1) == SGN(n2) ) { + if ( PL(nm1) == 1 && PL(nm2) == 1 ) { + t1 = BD(nm1)[0]; t = t1+BD(nm2)[0]; + if ( t < t1 ) { + m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1; + } else + UTON(t,m); + } else + addn(NM(n1),NM(n2),&m); + NTOQ(m,SGN(n1),*nr); + } else { + if ( PL(nm1) == 1 && PL(nm2) == 1 ) { + t1 = BD(nm1)[0]; t = t1-BD(nm2)[0]; + if ( !t ) + sgn = 0; + else if ( t > t1 ) { + sgn = -1; t = -((int)t); UTON(t,m); + } else { + sgn = 1; UTON(t,m); + } + } else + sgn = subn(NM(n1),NM(n2),&m); + if ( sgn ) { + if ( SGN(n1) == sgn ) + NTOQ(m,1,*nr); + else + NTOQ(m,-1,*nr); + } else + *nr = 0; + } + } else { + kmuln(NM(n1),DN(n2),&m); + if ( SGN(n1) == SGN(n2) ) { + sgn = SGN(n1); addn(m,NM(n2),&nm); + } else + sgn = SGN(n1)*subn(m,NM(n2),&nm); + if ( sgn ) { + dn = DN(n2); NDTOQ(nm,dn,sgn,*nr); + } else + *nr = 0; + } + else if ( INT(n2) ) { + kmuln(NM(n2),DN(n1),&m); + if ( SGN(n1) == SGN(n2) ) { + sgn = SGN(n1); addn(m,NM(n1),&nm); + } else + sgn = SGN(n1)*subn(NM(n1),m,&nm); + if ( sgn ) { + dn = DN(n1); NDTOQ(nm,dn,sgn,*nr); + } else + *nr = 0; + } else { + gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2); + kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2); + if ( SGN(n1) == SGN(n2) ) { + sgn = SGN(n1); addn(nm1,nm2,&nm3); + } else + sgn = SGN(n1)*subn(nm1,nm2,&nm3); + if ( sgn ) { + gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0); + kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn); + if ( UNIN(dn) ) + NTOQ(nm,sgn,*nr); + else + NDTOQ(nm,dn,sgn,*nr); + } else + *nr = 0; + } } -void subq(n1,n2,nr) -Q n1,n2,*nr; +void subq(Q n1,Q n2,Q *nr) { - Q m; + Q m; - if ( !n1 ) - if ( !n2 ) - *nr = 0; - else { - DUPQ(n2,*nr); SGN(*nr) = -SGN(n2); - } - else if ( !n2 ) - *nr = n1; - else if ( n1 == n2 ) - *nr = 0; - else { - DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr); - } + if ( !n1 ) + if ( !n2 ) + *nr = 0; + else { + DUPQ(n2,*nr); SGN(*nr) = -SGN(n2); + } + else if ( !n2 ) + *nr = n1; + else if ( n1 == n2 ) + *nr = 0; + else { + DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr); + } } -void mulq(n1,n2,nr) -Q n1,n2,*nr; +void mulq(Q n1,Q n2,Q *nr) { - N nm,nm1,nm2,dn,dn1,dn2,g; - int sgn; - unsigned u,l; + N nm,nm1,nm2,dn,dn1,dn2,g; + int sgn; + unsigned u,l; - if ( !n1 || !n2 ) *nr = 0; - else if ( INT(n1) ) - if ( INT(n2) ) { - nm1 = NM(n1); nm2 = NM(n2); - if ( PL(nm1) == 1 && PL(nm2) == 1 ) { - DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l) - if ( u ) { - nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u; - } else { - nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l; - } - } else - kmuln(nm1,nm2,&nm); - sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr); - } else { - gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn); - kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2); - if ( UNIN(dn) ) - NTOQ(nm,sgn,*nr); - else - NDTOQ(nm,dn,sgn,*nr); - } - else if ( INT(n2) ) { - gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn); - kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2); - if ( UNIN(dn) ) - NTOQ(nm,sgn,*nr); - else - NDTOQ(nm,dn,sgn,*nr); - } else { - gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2); - gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2); - kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2); - if ( UNIN(dn) ) - NTOQ(nm,sgn,*nr); - else - NDTOQ(nm,dn,sgn,*nr); - } + if ( !n1 || !n2 ) *nr = 0; + else if ( INT(n1) ) + if ( INT(n2) ) { + nm1 = NM(n1); nm2 = NM(n2); + if ( PL(nm1) == 1 && PL(nm2) == 1 ) { + DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l) + if ( u ) { + nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u; + } else { + nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l; + } + } else + kmuln(nm1,nm2,&nm); + sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr); + } else { + gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn); + kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2); + if ( UNIN(dn) ) + NTOQ(nm,sgn,*nr); + else + NDTOQ(nm,dn,sgn,*nr); + } + else if ( INT(n2) ) { + gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn); + kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2); + if ( UNIN(dn) ) + NTOQ(nm,sgn,*nr); + else + NDTOQ(nm,dn,sgn,*nr); + } else { + gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2); + gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2); + kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2); + if ( UNIN(dn) ) + NTOQ(nm,sgn,*nr); + else + NDTOQ(nm,dn,sgn,*nr); + } } -void divq(n1,n2,nq) -Q n1,n2,*nq; +void divq(Q n1,Q n2,Q *nq) { - Q m; + Q m; - if ( !n2 ) { - error("division by 0"); - *nq = 0; - return; - } else if ( !n1 ) - *nq = 0; - else if ( n1 == n2 ) - *nq = ONE; - else { - invq(n2,&m); mulq(n1,m,nq); - } + if ( !n2 ) { + error("division by 0"); + *nq = 0; + return; + } else if ( !n1 ) + *nq = 0; + else if ( n1 == n2 ) + *nq = ONE; + else { + invq(n2,&m); mulq(n1,m,nq); + } } -void invq(n,nr) -Q n,*nr; +void divsq(Q n1,Q n2,Q *nq) { - N nm,dn; + int sgn; + N tn; - if ( INT(n) ) - if ( UNIN(NM(n)) ) - if ( SGN(n) > 0 ) - *nr = ONE; - else { - nm = ONEN; NTOQ(nm,SGN(n),*nr); - } - else { - nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr); - } - else if ( UNIN(NM(n)) ) { - nm = DN(n); NTOQ(nm,SGN(n),*nr); - } else { - nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr); - } + if ( !n2 ) { + error("division by 0"); + *nq = 0; + return; + } else if ( !n1 ) + *nq = 0; + else if ( n1 == n2 ) + *nq = ONE; + else { + divsn(NM(n1),NM(n2),&tn); + sgn = SGN(n1)*SGN(n2); + NTOQ(tn,sgn,*nq); + } } -void chsgnq(n,nr) -Q n,*nr; +void invq(Q n,Q *nr) { - Q t; + N nm,dn; - if ( !n ) - *nr = 0; - else { - DUPQ(n,t); SGN(t) = -SGN(t); *nr = t; - } + if ( INT(n) ) + if ( UNIN(NM(n)) ) + if ( SGN(n) > 0 ) + *nr = ONE; + else { + nm = ONEN; NTOQ(nm,SGN(n),*nr); + } + else { + nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr); + } + else if ( UNIN(NM(n)) ) { + nm = DN(n); NTOQ(nm,SGN(n),*nr); + } else { + nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr); + } } -void pwrq(n1,n,nr) -Q n1,n,*nr; +void chsgnq(Q n,Q *nr) { - N nm,dn; - int sgn; + Q t; - if ( !n || UNIQ(n1) ) - *nr = ONE; - else if ( !n1 ) - *nr = 0; - else if ( !INT(n) ) { - error("can't calculate fractional power."); *nr = 0; - } else if ( MUNIQ(n1) ) - *nr = BD(NM(n))[0]%2 ? n1 : ONE; - else if ( PL(NM(n)) > 1 ) { - error("exponent too big."); *nr = 0; - } else { - sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1); - pwrn(NM(n1),BD(NM(n))[0],&nm); - if ( INT(n1) ) { - if ( UNIN(nm) ) - if ( sgn == 1 ) - *nr = ONE; - else - NTOQ(ONEN,sgn,*nr); - else if ( SGN(n) > 0 ) - NTOQ(nm,sgn,*nr); - else - NDTOQ(ONEN,nm,sgn,*nr); - } else { - pwrn(DN(n1),BD(NM(n))[0],&dn); - if ( SGN(n) > 0 ) - NDTOQ(nm,dn,sgn,*nr); - else if ( UNIN(nm) ) - NTOQ(dn,sgn,*nr); - else - NDTOQ(dn,nm,sgn,*nr); - } - } + if ( !n ) + *nr = 0; + else { + DUPQ(n,t); SGN(t) = -SGN(t); *nr = t; + } } -int cmpq(q1,q2) -Q q1,q2; +void pwrq(Q n1,Q n,Q *nr) { - int sgn; - N t,s; + N nm,dn; + int sgn; - if ( !q1 ) - if ( !q2 ) - return ( 0 ); - else - return ( -SGN(q2) ); - else if ( !q2 ) - return ( SGN(q1) ); - else if ( SGN(q1) != SGN(q2) ) - return ( SGN(q1) ); - else if ( INT(q1) ) - if ( INT(q2) ) { - sgn = cmpn(NM(q1),NM(q2)); - if ( !sgn ) - return ( 0 ); - else - return ( SGN(q1)==sgn?1:-1 ); - } else { - kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2)); - return ( SGN(q1) * sgn ); - } - else if ( INT(q2) ) { - kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t); - return ( SGN(q1) * sgn ); - } else { - kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s); - return ( SGN(q1) * sgn ); - } + if ( !n || UNIQ(n1) ) + *nr = ONE; + else if ( !n1 ) + *nr = 0; + else if ( !INT(n) ) { + error("can't calculate fractional power."); *nr = 0; + } else if ( MUNIQ(n1) ) + *nr = BD(NM(n))[0]%2 ? n1 : ONE; + else if ( PL(NM(n)) > 1 ) { + error("exponent too big."); *nr = 0; + } else { + sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1); + pwrn(NM(n1),BD(NM(n))[0],&nm); + if ( INT(n1) ) { + if ( UNIN(nm) ) + if ( sgn == 1 ) + *nr = ONE; + else + NTOQ(ONEN,sgn,*nr); + else if ( SGN(n) > 0 ) + NTOQ(nm,sgn,*nr); + else + NDTOQ(ONEN,nm,sgn,*nr); + } else { + pwrn(DN(n1),BD(NM(n))[0],&dn); + if ( SGN(n) > 0 ) + NDTOQ(nm,dn,sgn,*nr); + else if ( UNIN(nm) ) + NTOQ(dn,sgn,*nr); + else + NDTOQ(dn,nm,sgn,*nr); + } + } } -void remq(n,m,nr) -Q n,m; -Q *nr; +int cmpq(Q q1,Q q2) { - N q,r,s; + int sgn; + N t,s; - if ( !n ) { - *nr = 0; - return; - } - divn(NM(n),NM(m),&q,&r); - if ( !r ) - *nr = 0; - else if ( SGN(n) > 0 ) - NTOQ(r,1,*nr); - else { - subn(NM(m),r,&s); NTOQ(s,1,*nr); - } + if ( !q1 ) + if ( !q2 ) + return ( 0 ); + else + return ( -SGN(q2) ); + else if ( !q2 ) + return ( SGN(q1) ); + else if ( SGN(q1) != SGN(q2) ) + return ( SGN(q1) ); + else if ( INT(q1) ) + if ( INT(q2) ) { + sgn = cmpn(NM(q1),NM(q2)); + if ( !sgn ) + return ( 0 ); + else + return ( SGN(q1)==sgn?1:-1 ); + } else { + kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2)); + return ( SGN(q1) * sgn ); + } + else if ( INT(q2) ) { + kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t); + return ( SGN(q1) * sgn ); + } else { + kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s); + return ( SGN(q1) * sgn ); + } } -void mkbc(n,t) -int n; -Q *t; +void remq(Q n,Q m,Q *nr) { - int i; - N c,d; + N q,r,s; - for ( t[0] = ONE, i = 1; i <= n/2; i++ ) { - c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1; - kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]); - } - for ( ; i <= n; i++ ) - t[i] = t[n-i]; + if ( !n ) { + *nr = 0; + return; + } + divn(NM(n),NM(m),&q,&r); + if ( !r ) + *nr = 0; + else if ( SGN(n) > 0 ) + NTOQ(r,1,*nr); + else { + subn(NM(m),r,&s); NTOQ(s,1,*nr); + } } -void factorial(n,r) -int n; -Q *r; +/* t = [nC0 nC1 ... nCn] */ + +void mkbc(int n,Q *t) { - Q t,iq,s; - unsigned int i,m,m0; - unsigned int c; + int i; + N c,d; - if ( !n ) - *r = ONE; - else if ( n < 0 ) - *r = 0; - else { - for ( i = 1, t = ONE; (int)i <= n; ) { - for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) { - m0 = m; DM(m0,i,c,m) - } - if ( c ) { - m = m0; i--; - } - UTOQ(m,iq); mulq(t,iq,&s); t = s; - } - *r = t; - } + for ( t[0] = ONE, i = 1; i <= n/2; i++ ) { + c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1; + kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]); + } + for ( ; i <= n; i++ ) + t[i] = t[n-i]; } -void invl(a,mod,ar) -Q a,mod,*ar; +/* + * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+... + * + * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)] + * where W(k,l,i) = i! * kCi * lCi + */ + +void mkwc(int k,int l,Q *t) { - Q f1,f2,a1,a2,q,m,s; - N qn,rn; + int i,n,up,low; + N nm,d,c; - a1 = ONE; a2 = 0; - DUPQ(a,f1); SGN(f1) = 1; f2 = mod; + n = MIN(k,l); + for ( t[0] = ONE, i = 1; i <= n; i++ ) { + DM(k-i+1,l-i+1,up,low); + if ( up ) { + nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up; + } else { + nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low; + } + kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]); + } +} - while ( 1 ) { - divn(NM(f1),NM(f2),&qn,&rn); - if ( !qn ) - q = 0; - else - NTOQ(qn,1,q); +/* mod m table */ +/* XXX : should be optimized */ - f1 = f2; - if ( !rn ) - break; - else - NTOQ(rn,1,f2); +void mkwcm(int k,int l,int m,int *t) +{ + int i,n; + Q *s; - mulq(a2,q,&m); subq(a1,m,&s); a1 = a2; - if ( !s ) - a2 = 0; - else - remq(s,mod,&a2); - } - if ( SGN(a) < 0 ) - chsgnq(a2,&m); - else - m = a2; + n = MIN(k,l); + s = (Q *)ALLOCA((n+1)*sizeof(Q)); + mkwc(k,l,s); + for ( i = 0; i <= n; i++ ) { + t[i] = rem(NM(s[i]),m); + } +} - if ( SGN(m) >= 0 ) - *ar = m; - else - addq(m,mod,ar); +#if 0 +void factorial(int n,Q *r) +{ + Q t,iq,s; + unsigned int i,m,m0; + unsigned int c; + + if ( !n ) + *r = ONE; + else if ( n < 0 ) + *r = 0; + else { + for ( i = 1, t = ONE; (int)i <= n; ) { + for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) { + m0 = m; DM(m0,i,c,m) + } + if ( c ) { + m = m0; i--; + } + UTOQ(m,iq); mulq(t,iq,&s); t = s; + } + *r = t; + } } +#endif +void partial_factorial(int s,int e,N *r); + +void factorial(int n,Q *r) +{ + N nm; + + if ( !n ) + *r = ONE; + else if ( n < 0 ) + *r = 0; + else { + partial_factorial(1,n,&nm); + NTOQ(nm,1,*r); + } +} + +/* s*(s+1)*...*e */ +void partial_factorial(int s,int e,N *r) +{ + unsigned int len,i,m,m0,c; + unsigned int *p,*p1,*p2; + N nm,r1,r2; + + /* XXX */ + if ( e-s > 2000 ) { + m = (e+s)/2; + partial_factorial(s,m,&r1); + partial_factorial(m+1,e,&r2); + kmuln(r1,r2,r); + return; + } + /* find i s.t. 2^(i-1) < m <= 2^i */ + for ( i = 0, m = 1; m < e; m <<=1, i++ ); + /* XXX estimate of word length of the result */ + len = (i*(e-s+1)+1+31)/32; + p = ALLOCA((len+1)*sizeof(int)); + p1 = ALLOCA((len+1)*sizeof(int)); + p[0] = s; + len = 1; + for ( i = s+1; (int)i <= e; ) { + for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) { + m0 = m; DM(m0,i,c,m) + } + if ( c ) { + m = m0; i--; + } + bzero(p1,(len+1)*sizeof(int)); + muln_1(p,len,m,p1); + if ( p1[len] ) + len++; + p2 = p; p = p1; p1 = p2; + } + *r = nm = NALLOC(len); + bcopy(p,BD(nm),sizeof(int)*len); + PL(nm) = len; +} + +void invl(Q a,Q mod,Q *ar) +{ + Q f1,f2,a1,a2,q,m,s; + N qn,rn; + + a1 = ONE; a2 = 0; + DUPQ(a,f1); SGN(f1) = 1; f2 = mod; + + while ( 1 ) { + divn(NM(f1),NM(f2),&qn,&rn); + if ( !qn ) + q = 0; + else + NTOQ(qn,1,q); + + f1 = f2; + if ( !rn ) + break; + else + NTOQ(rn,1,f2); + + mulq(a2,q,&m); subq(a1,m,&s); a1 = a2; + if ( !s ) + a2 = 0; + else + remq(s,mod,&a2); + } + if ( SGN(a) < 0 ) + chsgnq(a2,&m); + else + m = a2; + + if ( SGN(m) >= 0 ) + *ar = m; + else + addq(m,mod,ar); +} + int kara_mag=100; -void kmuln(n1,n2,nr) -N n1,n2,*nr; +void kmuln(N n1,N n2,N *nr) { - N n,t,s,m,carry; - int d,d1,d2,len,i,l; - int *r,*r0; + N n,t,s,m,carry; + int d,d1,d2,len,i,l; + int *r,*r0; - if ( !n1 || !n2 ) { - *nr = 0; return; - } - d1 = PL(n1); d2 = PL(n2); - if ( (d1 < kara_mag) || (d2 < kara_mag) ) { - muln(n1,n2,nr); return; - } - if ( d1 < d2 ) { - n = n1; n1 = n2; n2 = n; - d = d1; d1 = d2; d2 = d; - } - if ( d2 > (d1+1)/2 ) { - kmulnmain(n1,n2,nr); return; - } - d = (d1/d2)+((d1%d2)!=0); - len = (d+1)*d2; - r0 = (int *)ALLOCA(len*sizeof(int)); - bzero((char *)r0,len*sizeof(int)); - for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) { - extractn(n1,i*d2,d2,&m); - if ( m ) { - kmuln(m,n2,&t); - addn(t,carry,&s); - copyn(s,d2,r); - extractn(s,d2,d2,&carry); - } else { - copyn(carry,d2,r); - carry = 0; - } - } - copyn(carry,d2,r); - for ( l = len - 1; !r0[l]; l-- ); - l++; - *nr = t = NALLOC(l); PL(t) = l; - bcopy((char *)r0,(char *)BD(t),l*sizeof(int)); + if ( !n1 || !n2 ) { + *nr = 0; return; + } + d1 = PL(n1); d2 = PL(n2); + if ( (d1 < kara_mag) || (d2 < kara_mag) ) { + muln(n1,n2,nr); return; + } + if ( d1 < d2 ) { + n = n1; n1 = n2; n2 = n; + d = d1; d1 = d2; d2 = d; + } + if ( d2 > (d1+1)/2 ) { + kmulnmain(n1,n2,nr); return; + } + d = (d1/d2)+((d1%d2)!=0); + len = (d+1)*d2; + r0 = (int *)ALLOCA(len*sizeof(int)); + bzero((char *)r0,len*sizeof(int)); + for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) { + extractn(n1,i*d2,d2,&m); + if ( m ) { + kmuln(m,n2,&t); + addn(t,carry,&s); + copyn(s,d2,r); + extractn(s,d2,d2,&carry); + } else { + copyn(carry,d2,r); + carry = 0; + } + } + copyn(carry,d2,r); + for ( l = len - 1; !r0[l]; l-- ); + l++; + *nr = t = NALLOC(l); PL(t) = l; + bcopy((char *)r0,(char *)BD(t),l*sizeof(int)); } -void extractn(n,index,len,nr) -N n; -int index,len; -N *nr; +void extractn(N n,int index,int len,N *nr) { - unsigned int *m; - int l; - N t; + unsigned int *m; + int l; + N t; - if ( !n ) { - *nr = 0; return; - } - m = BD(n)+index; - if ( (l = PL(n)-index) >= len ) { - for ( l = len - 1; (l >= 0) && !m[l]; l-- ); - l++; - } - if ( l <= 0 ) { - *nr = 0; return; - } else { - *nr = t = NALLOC(l); PL(t) = l; - bcopy((char *)m,(char *)BD(t),l*sizeof(int)); - } + if ( !n ) { + *nr = 0; return; + } + m = BD(n)+index; + if ( (l = PL(n)-index) >= len ) { + for ( l = len - 1; (l >= 0) && !m[l]; l-- ); + l++; + } + if ( l <= 0 ) { + *nr = 0; return; + } else { + *nr = t = NALLOC(l); PL(t) = l; + bcopy((char *)m,(char *)BD(t),l*sizeof(int)); + } } -void copyn(n,len,p) -N n; -int len; -int *p; +void copyn(N n,int len,int *p) { - if ( n ) - bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int)); + if ( n ) + bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int)); } -void dupn(n,p) -N n; -N p; +void dupn(N n,N p) { - if ( n ) - bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int)); + if ( n ) + bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int)); } -void kmulnmain(n1,n2,nr) -N n1,n2,*nr; +void kmulnmain(N n1,N n2,N *nr) { - int d1,d2,h,sgn,sgn1,sgn2,len; - N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; + int d1,d2,h,sgn,sgn1,sgn2,len; + N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; - d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2; - extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi); - extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi); - kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo); - len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len; - bzero((char *)BD(t1),len*sizeof(int)); - if ( lo ) - bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int)); - if ( hi ) - bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int)); + d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2; + extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi); + extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi); + kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo); + len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len; + bzero((char *)BD(t1),len*sizeof(int)); + if ( lo ) + bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int)); + if ( hi ) + bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int)); - addn(hi,lo,&mid1); - sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2); - kmuln(s1,s2,&mid2); sgn = sgn1*sgn2; - if ( sgn > 0 ) - addn(mid1,mid2,&mid); - else - subn(mid1,mid2,&mid); - if ( mid ) { - len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len; - bzero((char *)BD(t2),len*sizeof(int)); - bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int)); - addn(t1,t2,nr); - } else - *nr = t1; + addn(hi,lo,&mid1); + sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2); + kmuln(s1,s2,&mid2); sgn = sgn1*sgn2; + if ( sgn > 0 ) + addn(mid1,mid2,&mid); + else + subn(mid1,mid2,&mid); + if ( mid ) { + len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len; + bzero((char *)BD(t2),len*sizeof(int)); + bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int)); + addn(t1,t2,nr); + } else + *nr = t1; }