=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/int.c,v retrieving revision 1.1.1.1 retrieving revision 1.14 diff -u -p -r1.1.1.1 -r1.14 --- OpenXM_contrib2/asir2000/builtin/int.c 1999/12/03 07:39:07 1.1.1.1 +++ OpenXM_contrib2/asir2000/builtin/int.c 2018/03/29 01:32:50 1.14 @@ -1,4 +1,52 @@ -/* $OpenXM: OpenXM/src/asir99/builtin/int.c,v 1.1.1.1 1999/11/10 08:12:25 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/builtin/int.c,v 1.13 2015/08/14 13:51:54 fujimoto Exp $ +*/ #include "ca.h" #include "parse.h" #include "base.h" @@ -9,66 +57,61 @@ void Pup2_init_eg(), Pup2_show_eg(); void Piqr(), Pprime(), Plprime(), Pinttorat(); void Piand(), Pior(), Pixor(), Pishift(); void Pisqrt(); -void iand(), ior(), ixor(); -void isqrt(); void Plrandom(); void Pset_upkara(), Pset_uptkara(), Pset_up2kara(), Pset_upfft(); void Pmt_save(), Pmt_load(); void Psmall_jacobi(); void Pdp_set_mpi(); +void Pntoint32(),Pint32ton(); -#ifdef HMEXT void Pigcdbin(), Pigcdbmod(), PigcdEuc(), Pigcdacc(), Pigcdcntl(); void Pihex(); void Pimaxrsh(), Pilen(); void Ptype_t_NB(); -#endif /* HMEXT */ - struct ftab int_tab[] = { - {"dp_set_mpi",Pdp_set_mpi,-1}, - {"isqrt",Pisqrt,1}, - {"idiv",Pidiv,2}, - {"irem",Pirem,2}, - {"iqr",Piqr,2}, - {"igcd",Pigcd,-2}, - {"ilcm",Pilcm,2}, - {"up2_inv",Pup2_inv,2}, - {"up2_init_eg",Pup2_init_eg,0}, - {"up2_show_eg",Pup2_show_eg,0}, - {"type_t_NB",Ptype_t_NB,2}, - {"gf2nton",Pgf2nton,1}, - {"ntogf2n",Pntogf2n,1}, - {"set_upkara",Pset_upkara,-1}, - {"set_uptkara",Pset_uptkara,-1}, - {"set_up2kara",Pset_up2kara,-1}, - {"set_upfft",Pset_upfft,-1}, - {"inv",Pinv,2}, - {"inttorat",Pinttorat,3}, - {"fac",Pfac,1}, - {"prime",Pprime,1}, - {"lprime",Plprime,1}, - {"random",Prandom,-1}, - {"lrandom",Plrandom,1}, - {"iand",Piand,2}, - {"ior",Pior,2}, - {"ixor",Pixor,2}, - {"ishift",Pishift,2}, - {"small_jacobi",Psmall_jacobi,2}, -#ifdef HMEXT - {"igcdbin",Pigcdbin,2}, /* HM@CCUT extension */ - {"igcdbmod",Pigcdbmod,2}, /* HM@CCUT extension */ - {"igcdeuc",PigcdEuc,2}, /* HM@CCUT extension */ - {"igcdacc",Pigcdacc,2}, /* HM@CCUT extension */ - {"igcdcntl",Pigcdcntl,-1}, /* HM@CCUT extension */ - {"ihex",Pihex,1}, /* HM@CCUT extension */ - {"imaxrsh",Pimaxrsh,1}, /* HM@CCUT extension */ - {"ilen",Pilen,1}, /* HM@CCUT extension */ -#endif /* HMEXT */ - {"mt_save",Pmt_save,1}, - {"mt_load",Pmt_load,1}, - {0,0,0}, + {"dp_set_mpi",Pdp_set_mpi,-1}, + {"isqrt",Pisqrt,1}, + {"idiv",Pidiv,2}, + {"irem",Pirem,2}, + {"iqr",Piqr,2}, + {"igcd",Pigcd,-2}, + {"ilcm",Pilcm,2}, + {"up2_inv",Pup2_inv,2}, + {"up2_init_eg",Pup2_init_eg,0}, + {"up2_show_eg",Pup2_show_eg,0}, + {"type_t_NB",Ptype_t_NB,2}, + {"gf2nton",Pgf2nton,1}, + {"ntogf2n",Pntogf2n,1}, + {"set_upkara",Pset_upkara,-1}, + {"set_uptkara",Pset_uptkara,-1}, + {"set_up2kara",Pset_up2kara,-1}, + {"set_upfft",Pset_upfft,-1}, + {"inv",Pinv,2}, + {"inttorat",Pinttorat,3}, + {"fac",Pfac,1}, + {"prime",Pprime,1}, + {"lprime",Plprime,1}, + {"random",Prandom,-1}, + {"lrandom",Plrandom,1}, + {"iand",Piand,2}, + {"ior",Pior,2}, + {"ixor",Pixor,2}, + {"ishift",Pishift,2}, + {"small_jacobi",Psmall_jacobi,2}, + + {"igcdbin",Pigcdbin,2}, /* HM@CCUT extension */ + {"igcdbmod",Pigcdbmod,2}, /* HM@CCUT extension */ + {"igcdeuc",PigcdEuc,2}, /* HM@CCUT extension */ + {"igcdacc",Pigcdacc,2}, /* HM@CCUT extension */ + {"igcdcntl",Pigcdcntl,-1}, /* HM@CCUT extension */ + + {"mt_save",Pmt_save,1}, + {"mt_load",Pmt_load,1}, + {"ntoint32",Pntoint32,1}, + {"int32ton",Pint32ton,1}, + {0,0,0}, }; static int is_prime_small(unsigned int); @@ -76,94 +119,115 @@ static unsigned int gcd_small(unsigned int,unsigned in int TypeT_NB_check(unsigned int, unsigned int); int mpi_mag; -void Pdp_set_mpi(arg,rp) -NODE arg; -Q *rp; +void Pntoint32(NODE arg,USINT *rp) { - if ( arg ) { - asir_assert(ARG0(arg),O_N,"dp_set_mpi"); - mpi_mag = QTOS((Q)ARG0(arg)); - } - STOQ(mpi_mag,*rp); + Q q; + unsigned int t; + + asir_assert(ARG0(arg),O_N,"ntoint32"); + q = (Q)ARG0(arg); + if ( !q ) { + MKUSINT(*rp,0); + return; + } + if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 ) + error("ntoint32 : invalid argument"); + t = BD(NM(q))[0]; + if ( SGN(q) < 0 ) + t = -(int)t; + MKUSINT(*rp,t); } -void Psmall_jacobi(arg,rp) -NODE arg; -Q *rp; +void Pint32ton(NODE arg,Q *rp) { - Q a,m; - int a0,m0,s; + int t; - a = (Q)ARG0(arg); - m = (Q)ARG1(arg); - asir_assert(a,O_N,"small_jacobi"); - asir_assert(m,O_N,"small_jacobi"); - if ( !a ) - *rp = ONE; - else if ( !m || !INT(m) || !INT(a) - || PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) ) - error("small_jacobi : invalid input"); - else { - a0 = QTOS(a); m0 = QTOS(m); - s = small_jacobi(a0,m0); - STOQ(s,*rp); - } + asir_assert(ARG0(arg),O_USINT,"int32ton"); + t = (int)BDY((USINT)ARG0(arg)); + STOQ(t,*rp); } -int small_jacobi(a,m) -int a,m; +void Pdp_set_mpi(NODE arg,Q *rp) { - int m4,m8,a4,j1,i,s; + if ( arg ) { + asir_assert(ARG0(arg),O_N,"dp_set_mpi"); + mpi_mag = QTOS((Q)ARG0(arg)); + } + STOQ(mpi_mag,*rp); +} - a %= m; - if ( a == 0 || a == 1 ) - return 1; - else if ( a < 0 ) { - j1 = small_jacobi(-a,m); - m4 = m%4; - return m4==1?j1:-j1; - } else { - for ( i = 0; a && !(a&1); i++, a >>= 1 ); - if ( i&1 ) { - m8 = m%8; - s = (m8==1||m8==7) ? 1 : -1; - } else - s = 1; - /* a, m are odd */ - j1 = small_jacobi(m%a,a); - m4 = m%4; a4 = a%4; - s *= (m4==1||a4==1) ? 1 : -1; - return j1*s; - } +void Psmall_jacobi(NODE arg,Q *rp) +{ + Q a,m; + int a0,m0,s; + + a = (Q)ARG0(arg); + m = (Q)ARG1(arg); + asir_assert(a,O_N,"small_jacobi"); + asir_assert(m,O_N,"small_jacobi"); + if ( !a ) + *rp = ONE; + else if ( !m || !INT(m) || !INT(a) + || PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) ) + error("small_jacobi : invalid input"); + else { + a0 = QTOS(a); m0 = QTOS(m); + s = small_jacobi(a0,m0); + STOQ(s,*rp); + } } -void Ptype_t_NB(arg,rp) -NODE arg; -Q *rp; +int small_jacobi(int a,int m) { - if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)))) - *rp = ONE; - else - *rp = 0; + int m4,m8,a4,j1,i,s; + + a %= m; + if ( a == 0 || a == 1 ) + return 1; + else if ( a < 0 ) { + j1 = small_jacobi(-a,m); + m4 = m%4; + return m4==1?j1:-j1; + } else { + for ( i = 0; a && !(a&1); i++, a >>= 1 ); + if ( i&1 ) { + m8 = m%8; + s = (m8==1||m8==7) ? 1 : -1; + } else + s = 1; + /* a, m are odd */ + j1 = small_jacobi(m%a,a); + m4 = m%4; a4 = a%4; + s *= (m4==1||a4==1) ? 1 : -1; + return j1*s; + } } +void Ptype_t_NB(NODE arg,Q *rp) +{ + if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg)))) + *rp = ONE; + else + *rp = 0; +} + int TypeT_NB_check(unsigned int m, unsigned int t) { - unsigned int p,k,u,h,d; + unsigned int p,k,u,h,d; - if ( !(m%8) ) - return 0; - p = t*m+1; - if ( !is_prime_small(p) ) - return 0; - for ( k = 1, u = 2%p; ; k++ ) - if ( u == 1 ) - break; - else - u = (2*u)%p; - h = t*m/k; - d = gcd_small(h,m); - return d == 1 ? 1 : 0; + if ( !(m%8) ) + return 0; + p = t*m+1; + if ( !is_prime_small(p) ) + return 0; + for ( k = 1, u = 2%p; ; k++ ) + if ( u == 1 ) + break; + else + u = (2*u)%p; + h = t*m/k; + d = gcd_small(h,m); + return d == 1 ? 1 : 0; } /* @@ -174,18 +238,18 @@ int TypeT_NB_check(unsigned int m, unsigned int t) static int is_prime_small(unsigned int a) { - unsigned int b,t,i; + unsigned int b,t,i; - if ( !(a%2) ) return 0; - for ( t = a, i = 0; t; i++, t >>= 1 ); - /* b >= sqrt(a) */ - b = 1<<((i+1)/2); + if ( !(a%2) ) return 0; + for ( t = a, i = 0; t; i++, t >>= 1 ); + /* b >= sqrt(a) */ + b = 1<<((i+1)/2); - /* divisibility test by all odd numbers <= b */ - for ( i = 3; i <= b; i += 2 ) - if ( !(a%i) ) - return 0; - return 1; + /* divisibility test by all odd numbers <= b */ + for ( i = 3; i <= b; i += 2 ) + if ( !(a%i) ) + return 0; + return 1; } /* @@ -197,595 +261,574 @@ static int is_prime_small(unsigned int a) static unsigned int gcd_small(unsigned int a,unsigned int b) { - unsigned int t; + unsigned int t; - if ( b > a ) { - t = a; a = b; b = t; - } - /* Euclid's algorithm */ - while ( 1 ) - if ( !(t = a%b) ) return b; - else { - a = b; b = t; - } + if ( b > a ) { + t = a; a = b; b = t; + } + /* Euclid's algorithm */ + while ( 1 ) + if ( !(t = a%b) ) return b; + else { + a = b; b = t; + } } -void Pmt_save(arg,rp) -NODE arg; -Q *rp; +void Pmt_save(NODE arg,Q *rp) { - int ret; + int ret; - ret = mt_save(BDY((STRING)ARG0(arg))); - STOQ(ret,*rp); + ret = mt_save(BDY((STRING)ARG0(arg))); + STOQ(ret,*rp); } -void Pmt_load(arg,rp) -NODE arg; -Q *rp; +void Pmt_load(NODE arg,Q *rp) { - int ret; + int ret; - ret = mt_load(BDY((STRING)ARG0(arg))); - STOQ(ret,*rp); + ret = mt_load(BDY((STRING)ARG0(arg))); + STOQ(ret,*rp); } -void Pisqrt(arg,rp) -NODE arg; -Q *rp; +void Pisqrt(NODE arg,Q *rp) { - Q a; - N r; + Q a; + N r; - a = (Q)ARG0(arg); - asir_assert(a,O_N,"isqrt"); - if ( !a ) - *rp = 0; - else if ( SGN(a) < 0 ) - error("isqrt : negative argument"); - else { - isqrt(NM(a),&r); - NTOQ(r,1,*rp); - } + a = (Q)ARG0(arg); + asir_assert(a,O_N,"isqrt"); + if ( !a ) + *rp = 0; + else if ( SGN(a) < 0 ) + error("isqrt : negative argument"); + else { + isqrt(NM(a),&r); + NTOQ(r,1,*rp); + } } -void Pidiv(arg,rp) -NODE arg; -Obj *rp; +void Pidiv(NODE arg,Obj *rp) { - N q,r; - Q a; - Q dnd,dvr; - - dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); - asir_assert(dnd,O_N,"idiv"); - asir_assert(dvr,O_N,"idiv"); - if ( !dvr ) - error("idiv: division by 0"); - else if ( !dnd ) - *rp = 0; - else { - divn(NM(dnd),NM(dvr),&q,&r); - NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a; - } + N q,r; + Q a; + Q dnd,dvr; + + dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); + asir_assert(dnd,O_N,"idiv"); + asir_assert(dvr,O_N,"idiv"); + if ( !dvr ) + error("idiv: division by 0"); + else if ( !dnd ) + *rp = 0; + else { + divn(NM(dnd),NM(dvr),&q,&r); + NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a; + } } -void Pirem(arg,rp) -NODE arg; -Obj *rp; +void Pirem(NODE arg,Obj *rp) { - N q,r; - Q a; - Q dnd,dvr; - - dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); - asir_assert(dnd,O_N,"irem"); - asir_assert(dvr,O_N,"irem"); - if ( !dvr ) - error("irem: division by 0"); - else if ( !dnd ) - *rp = 0; - else { - divn(NM(dnd),NM(dvr),&q,&r); - NTOQ(r,SGN(dnd),a); *rp = (Obj)a; - } + N q,r; + Q a; + Q dnd,dvr; + + dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); + asir_assert(dnd,O_N,"irem"); + asir_assert(dvr,O_N,"irem"); + if ( !dvr ) + error("irem: division by 0"); + else if ( !dnd ) + *rp = 0; + else { + divn(NM(dnd),NM(dvr),&q,&r); + NTOQ(r,SGN(dnd),a); *rp = (Obj)a; + } } -void Piqr(arg,rp) -NODE arg; -LIST *rp; +void Piqr(NODE arg,LIST *rp) { - N q,r; - Q a,b; - Q dnd,dvr; - NODE node1,node2; + N q,r; + Q a,b; + Q dnd,dvr; + NODE node1,node2; - dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); - if ( !dvr ) - error("iqr: division by 0"); - else if ( !dnd ) - a = b = 0; - else if ( OID(dnd) == O_VECT ) { - iqrv((VECT)dnd,dvr,rp); return; - } else { - asir_assert(dnd,O_N,"iqr"); - asir_assert(dvr,O_N,"iqr"); - divn(NM(dnd),NM(dvr),&q,&r); - NTOQ(q,SGN(dnd)*SGN(dvr),a); - NTOQ(r,SGN(dnd),b); - } - MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1); + dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg); + if ( !dvr ) + error("iqr: division by 0"); + else if ( !dnd ) + a = b = 0; + else if ( OID(dnd) == O_VECT ) { + iqrv((VECT)dnd,dvr,rp); return; + } else { + asir_assert(dnd,O_N,"iqr"); + asir_assert(dvr,O_N,"iqr"); + divn(NM(dnd),NM(dvr),&q,&r); + NTOQ(q,SGN(dnd)*SGN(dvr),a); + NTOQ(r,SGN(dnd),b); + } + MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1); } -void Pinttorat(arg,rp) -NODE arg; -LIST *rp; +void Pinttorat(NODE arg,LIST *rp) { - Q cq,qq,t,u1,v1,r1,nm; - N m,b,q,r,c,u2,v2,r2; - NODE node1,node2; - P p; + Q cq,qq,t,u1,v1,r1,nm; + N m,b,q,r,c,u2,v2,r2; + NODE node1,node2; + P p; - asir_assert(ARG0(arg),O_N,"inttorat"); - asir_assert(ARG1(arg),O_N,"inttorat"); - asir_assert(ARG2(arg),O_N,"inttorat"); - cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg)); - if ( !cq ) { - MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1); - return; - } - divn(NM(cq),m,&q,&r); - if ( !r ) { - MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1); - return; - } else if ( SGN(cq) < 0 ) { - subn(m,r,&c); - } else - c = r; - u1 = 0; v1 = ONE; u2 = m; v2 = c; - while ( cmpn(v2,b) >= 0 ) { - divn(u2,v2,&q,&r2); u2 = v2; v2 = r2; - NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1; - } - if ( cmpn(NM(v1),b) >= 0 ) - *rp = 0; - else { - if ( SGN(v1) < 0 ) { - chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm); - } else - NTOQ(v2,1,nm); - MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1); - } + asir_assert(ARG0(arg),O_N,"inttorat"); + asir_assert(ARG1(arg),O_N,"inttorat"); + asir_assert(ARG2(arg),O_N,"inttorat"); + cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg)); + if ( !cq ) { + MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1); + return; + } + divn(NM(cq),m,&q,&r); + if ( !r ) { + MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1); + return; + } else if ( SGN(cq) < 0 ) { + subn(m,r,&c); + } else + c = r; + u1 = 0; v1 = ONE; u2 = m; v2 = c; + while ( cmpn(v2,b) >= 0 ) { + divn(u2,v2,&q,&r2); u2 = v2; v2 = r2; + NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1; + } + if ( cmpn(NM(v1),b) >= 0 ) + *rp = 0; + else { + if ( SGN(v1) < 0 ) { + chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm); + } else + NTOQ(v2,1,nm); + MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1); + } } -void Pigcd(arg,rp) -NODE arg; -Q *rp; +void Pigcd(NODE arg,Q *rp) { - N g; - Q n1,n2,a; - - if ( argc(arg) == 1 ) { - igcdv((VECT)ARG0(arg),rp); - return; - } - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"igcd"); - asir_assert(n2,O_N,"igcd"); - if ( !n1 ) - *rp = n2; - else if ( !n2 ) - *rp = n1; - else { - gcdn(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = a; - } + N g; + Q n1,n2,a; + + if ( argc(arg) == 1 ) { + igcdv((VECT)ARG0(arg),rp); + return; + } + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"igcd"); + asir_assert(n2,O_N,"igcd"); + if ( !n1 ) + *rp = n2; + else if ( !n2 ) + *rp = n1; + else { + gcdn(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = a; + } } -int comp_n(a,b) -N *a,*b; +int comp_n(N *a,N *b) { - return cmpn(*a,*b); + return cmpn(*a,*b); } -void iqrv(a,dvr,rp) -VECT a; -Q dvr; -LIST *rp; +void iqrv(VECT a,Q dvr,LIST *rp) { - int i,n; - VECT q,r; - Q dnd,qi,ri; - Q *b; - N qn,rn; - NODE n0,n1; + int i,n; + VECT q,r; + Q dnd,qi,ri; + Q *b; + N qn,rn; + NODE n0,n1; - if ( !dvr ) - error("iqrv: division by 0"); - n = a->len; b = (Q *)BDY(a); - MKVECT(q,n); MKVECT(r,n); - for ( i = 0; i < n; i++ ) { - dnd = b[i]; - if ( !dnd ) { - qi = ri = 0; - } else { - divn(NM(dnd),NM(dvr),&qn,&rn); - NTOQ(qn,SGN(dnd)*SGN(dvr),qi); - NTOQ(rn,SGN(dnd),ri); - } - BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri; - } - MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0); + if ( !dvr ) + error("iqrv: division by 0"); + n = a->len; b = (Q *)BDY(a); + MKVECT(q,n); MKVECT(r,n); + for ( i = 0; i < n; i++ ) { + dnd = b[i]; + if ( !dnd ) { + qi = ri = 0; + } else { + divn(NM(dnd),NM(dvr),&qn,&rn); + NTOQ(qn,SGN(dnd)*SGN(dvr),qi); + NTOQ(rn,SGN(dnd),ri); + } + BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri; + } + MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0); } -void igcdv(a,rp) -VECT a; -Q *rp; +/* + * gcd = GCD(a,b), ca = a/g, cb = b/g + */ + +void igcd_cofactor(Q a,Q b,Q *gcd,Q *ca,Q *cb) { - int i,j,n,nz; - N g,gt,q,r; - N *c; + N gn,tn; + + if ( !a ) { + if ( !b ) + error("igcd_cofactor : invalid input"); + else { + *ca = 0; + *cb = ONE; + *gcd = b; + } + } else if ( !b ) { + *ca = ONE; + *cb = 0; + *gcd = a; + } else { + gcdn(NM(a),NM(b),&gn); NTOQ(gn,1,*gcd); + divsn(NM(a),gn,&tn); NTOQ(tn,SGN(a),*ca); + divsn(NM(b),gn,&tn); NTOQ(tn,SGN(b),*cb); + } +} - n = a->len; - c = (N *)ALLOCA(n*sizeof(N)); - for ( i = 0; i < n; i++ ) - c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0; - qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n); - for ( ; n && ! *c; n--, c++ ); +void igcdv(VECT a,Q *rp) +{ + int i,j,n,nz; + N g,gt,q,r; + N *c; - if ( !n ) { - *rp = 0; return; - } else if ( n == 1 ) { - NTOQ(*c,1,*rp); return; - } - gcdn(c[0],c[1],&g); - for ( i = 2; i < n; i++ ) { - divn(c[i],g,&q,&r); - gcdn(g,r,>); - if ( !cmpn(g,gt) ) { - for ( j = i+1, nz = 0; j < n; j++ ) { - divn(c[j],g,&q,&r); c[j] = r; - if ( r ) - nz = 1; - } - } else - g = gt; - } - NTOQ(g,1,*rp); + n = a->len; + c = (N *)ALLOCA(n*sizeof(N)); + for ( i = 0; i < n; i++ ) + c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0; + qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n); + for ( ; n && ! *c; n--, c++ ); + + if ( !n ) { + *rp = 0; return; + } else if ( n == 1 ) { + NTOQ(*c,1,*rp); return; + } + gcdn(c[0],c[1],&g); +#if 0 + for ( i = 2; i < n; i++ ) { + divn(c[i],g,&q,&r); + gcdn(g,r,>); + if ( !cmpn(g,gt) ) { + for ( j = i+1, nz = 0; j < n; j++ ) { + divn(c[j],g,&q,&r); c[j] = r; + if ( r ) + nz = 1; + } + } else + g = gt; + } +#else + for ( i = 2; i < n; i++ ) { + gcdn(g,c[i],>); g = gt; + } +#endif + NTOQ(g,1,*rp); } -void igcdv_estimate(a,rp) -VECT a; -Q *rp; +void igcdv_estimate(VECT a,Q *rp) { - int n,i,m; - N s,t,u,g; - Q *q; + int n,i,m; + N s,t,u,g; + Q *q; - n = a->len; q = (Q *)a->body; - if ( n == 1 ) { - NTOQ(NM(q[0]),1,*rp); return; - } + n = a->len; q = (Q *)a->body; + if ( n == 1 ) { + NTOQ(NM(q[0]),1,*rp); return; + } - m = n/2; - for ( i = 0 , s = 0; i < m; i++ ) { - addn(s,NM(q[i]),&u); s = u; - } - for ( t = 0; i < n; i++ ) { - addn(t,NM(q[i]),&u); t = u; - } - gcdn(s,t,&g); NTOQ(g,1,*rp); + m = n/2; + for ( i = 0 , s = 0; i < m; i++ ) { + addn(s,NM(q[i]),&u); s = u; + } + for ( t = 0; i < n; i++ ) { + addn(t,NM(q[i]),&u); t = u; + } + gcdn(s,t,&g); NTOQ(g,1,*rp); } -void Pilcm(arg,rp) -NODE arg; -Obj *rp; +void Pilcm(NODE arg,Obj *rp) { - N g,q,r,l; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"ilcm"); - asir_assert(n2,O_N,"ilcm"); - if ( !n1 || !n2 ) - *rp = 0; - else { - gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r); - muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a; - } + N g,q,r,l; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"ilcm"); + asir_assert(n2,O_N,"ilcm"); + if ( !n1 || !n2 ) + *rp = 0; + else { + gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r); + muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a; + } } -void Piand(arg,rp) -NODE arg; -Q *rp; +void Piand(NODE arg,Q *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"iand"); - asir_assert(n2,O_N,"iand"); - if ( !n1 || !n2 ) - *rp = 0; - else { - iand(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"iand"); + asir_assert(n2,O_N,"iand"); + if ( !n1 || !n2 ) + *rp = 0; + else { + iand(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = a; + } } -void Pior(arg,rp) -NODE arg; -Q *rp; +void Pior(NODE arg,Q *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"ior"); - asir_assert(n2,O_N,"ior"); - if ( !n1 ) - *rp = n2; - else if ( !n2 ) - *rp = n1; - else { - ior(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"ior"); + asir_assert(n2,O_N,"ior"); + if ( !n1 ) + *rp = n2; + else if ( !n2 ) + *rp = n1; + else { + ior(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = a; + } } -void Pixor(arg,rp) -NODE arg; -Q *rp; +void Pixor(NODE arg,Q *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"ixor"); - asir_assert(n2,O_N,"ixor"); - if ( !n1 ) - *rp = n2; - else if ( !n2 ) - *rp = n1; - else { - ixor(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"ixor"); + asir_assert(n2,O_N,"ixor"); + if ( !n1 ) + *rp = n2; + else if ( !n2 ) + *rp = n1; + else { + ixor(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = a; + } } -void Pishift(arg,rp) -NODE arg; -Q *rp; +void Pishift(NODE arg,Q *rp) { - N g; - Q n1,s,a; - - n1 = (Q)ARG0(arg); s = (Q)ARG1(arg); - asir_assert(n1,O_N,"ixor"); - asir_assert(s,O_N,"ixor"); - if ( !n1 ) - *rp = 0; - else if ( !s ) - *rp = n1; - else { - bshiftn(NM(n1),QTOS(s),&g); - NTOQ(g,1,a); *rp = a; - } + N g; + Q n1,s,a; + + n1 = (Q)ARG0(arg); s = (Q)ARG1(arg); + asir_assert(n1,O_N,"ixor"); + asir_assert(s,O_N,"ixor"); + if ( !n1 ) + *rp = 0; + else if ( !s ) + *rp = n1; + else { + bshiftn(NM(n1),QTOS(s),&g); + NTOQ(g,1,a); *rp = a; + } } -void isqrt(a,r) -N a,*r; +void isqrt(N a,N *r) { - int k; - N x,t,x2,xh,quo,rem; + int k; + N x,t,x2,xh,quo,rem; - if ( !a ) - *r = 0; - else if ( UNIN(a) ) - *r = ONEN; - else { - k = n_bits(a); /* a <= 2^k-1 */ - bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */ - while ( 1 ) { - pwrn(x,2,&t); - if ( cmpn(t,a) <= 0 ) { - *r = x; return; - } else { - if ( BD(x)[0] & 1 ) - addn(x,a,&t); - else - t = a; - bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem); - bshiftn(x,1,&xh); addn(quo,xh,&x); - } - } - } + if ( !a ) + *r = 0; + else if ( UNIN(a) ) + *r = ONEN; + else { + k = n_bits(a); /* a <= 2^k-1 */ + bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */ + while ( 1 ) { + pwrn(x,2,&t); + if ( cmpn(t,a) <= 0 ) { + *r = x; return; + } else { + if ( BD(x)[0] & 1 ) + addn(x,a,&t); + else + t = a; + bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem); + bshiftn(x,1,&xh); addn(quo,xh,&x); + } + } + } } -void iand(n1,n2,r) -N n1,n2,*r; +void iand(N n1,N n2,N *r) { - int d1,d2,d,i; - N nr; - int *p1,*p2,*pr; + int d1,d2,d,i; + N nr; + int *p1,*p2,*pr; - d1 = PL(n1); d2 = PL(n2); - d = MIN(d1,d2); - nr = NALLOC(d); - for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ ) - pr[i] = p1[i] & p2[i]; - for ( i = d-1; i >= 0 && !pr[i]; i-- ); - if ( i < 0 ) - *r = 0; - else { - PL(nr) = i+1; *r = nr; - } + d1 = PL(n1); d2 = PL(n2); + d = MIN(d1,d2); + nr = NALLOC(d); + for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ ) + pr[i] = p1[i] & p2[i]; + for ( i = d-1; i >= 0 && !pr[i]; i-- ); + if ( i < 0 ) + *r = 0; + else { + PL(nr) = i+1; *r = nr; + } } -void ior(n1,n2,r) -N n1,n2,*r; +void ior(N n1,N n2,N *r) { - int d1,d2,i; - N nr; - int *p1,*p2,*pr; + int d1,d2,i; + N nr; + int *p1,*p2,*pr; - if ( PL(n1) < PL(n2) ) { - nr = n1; n1 = n2; n2 = nr; - } - d1 = PL(n1); d2 = PL(n2); - *r = nr = NALLOC(d1); - for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ ) - pr[i] = p1[i] | p2[i]; - for ( ; i < d1; i++ ) - pr[i] = p1[i]; - for ( i = d1-1; i >= 0 && !pr[i]; i-- ); - if ( i < 0 ) - *r = 0; - else { - PL(nr) = i+1; *r = nr; - } + if ( PL(n1) < PL(n2) ) { + nr = n1; n1 = n2; n2 = nr; + } + d1 = PL(n1); d2 = PL(n2); + *r = nr = NALLOC(d1); + for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ ) + pr[i] = p1[i] | p2[i]; + for ( ; i < d1; i++ ) + pr[i] = p1[i]; + for ( i = d1-1; i >= 0 && !pr[i]; i-- ); + if ( i < 0 ) + *r = 0; + else { + PL(nr) = i+1; *r = nr; + } } -void ixor(n1,n2,r) -N n1,n2,*r; +void ixor(N n1,N n2,N *r) { - int d1,d2,i; - N nr; - int *p1,*p2,*pr; + int d1,d2,i; + N nr; + int *p1,*p2,*pr; - if ( PL(n1) < PL(n2) ) { - nr = n1; n1 = n2; n2 = nr; - } - d1 = PL(n1); d2 = PL(n2); - *r = nr = NALLOC(d1); - for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ ) - pr[i] = p1[i] ^ p2[i]; - for ( ; i < d1; i++ ) - pr[i] = p1[i]; - for ( i = d1-1; i >= 0 && !pr[i]; i-- ); - if ( i < 0 ) - *r = 0; - else { - PL(nr) = i+1; *r = nr; - } + if ( PL(n1) < PL(n2) ) { + nr = n1; n1 = n2; n2 = nr; + } + d1 = PL(n1); d2 = PL(n2); + *r = nr = NALLOC(d1); + for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ ) + pr[i] = p1[i] ^ p2[i]; + for ( ; i < d1; i++ ) + pr[i] = p1[i]; + for ( i = d1-1; i >= 0 && !pr[i]; i-- ); + if ( i < 0 ) + *r = 0; + else { + PL(nr) = i+1; *r = nr; + } } -void Pup2_init_eg(rp) -Obj *rp; +void Pup2_init_eg(Obj *rp) { - up2_init_eg(); - *rp = 0; + up2_init_eg(); + *rp = 0; } -void Pup2_show_eg(rp) -Obj *rp; +void Pup2_show_eg(Obj *rp) { - up2_show_eg(); - *rp = 0; + up2_show_eg(); + *rp = 0; } -void Pgf2nton(arg,rp) -NODE arg; -Q *rp; +void Pgf2nton(NODE arg,Q *rp) { - if ( !ARG0(arg) ) - *rp = 0; - else - up2ton(((GF2N)ARG0(arg))->body,rp); + if ( !ARG0(arg) ) + *rp = 0; + else + up2ton(((GF2N)ARG0(arg))->body,rp); } -void Pntogf2n(arg,rp) -NODE arg; -GF2N *rp; +void Pntogf2n(NODE arg,GF2N *rp) { - UP2 t; + UP2 t; - ntoup2(ARG0(arg),&t); - MKGF2N(t,*rp); + ntoup2(ARG0(arg),&t); + MKGF2N(t,*rp); } -void Pup2_inv(arg,rp) -NODE arg; -P *rp; +void Pup2_inv(NODE arg,P *rp) { - UP2 a,b,t; + UP2 a,b,t; - ptoup2(ARG0(arg),&a); - ptoup2(ARG1(arg),&b); - invup2(a,b,&t); - up2top(t,rp); + ptoup2(ARG0(arg),&a); + ptoup2(ARG1(arg),&b); + invup2(a,b,&t); + up2top(t,rp); } -void Pinv(arg,rp) -NODE arg; -Num *rp; +void Pinv(NODE arg,Num *rp) { - Num n; - Q mod; - MQ r; - int inv; + Num n; + Q mod; + MQ r; + int inv; - n = (Num)ARG0(arg); mod = (Q)ARG1(arg); - asir_assert(n,O_N,"inv"); - asir_assert(mod,O_N,"inv"); - if ( !n || !mod ) - error("inv: invalid input"); - else - switch ( NID(n) ) { - case N_Q: - invl((Q)n,mod,(Q *)rp); - break; - case N_M: - inv = invm(CONT((MQ)n),QTOS(mod)); - STOMQ(inv,r); - *rp = (Num)r; - break; - default: - error("inv: invalid input"); - } + n = (Num)ARG0(arg); mod = (Q)ARG1(arg); + asir_assert(n,O_N,"inv"); + asir_assert(mod,O_N,"inv"); + if ( !n || !mod ) + error("inv: invalid input"); + else + switch ( NID(n) ) { + case N_Q: + invl((Q)n,mod,(Q *)rp); + break; + case N_M: + inv = invm(CONT((MQ)n),QTOS(mod)); + STOMQ(inv,r); + *rp = (Num)r; + break; + default: + error("inv: invalid input"); + } } -void Pfac(arg,rp) -NODE arg; -Q *rp; +void Pfac(NODE arg,Q *rp) { - asir_assert(ARG0(arg),O_N,"fac"); - factorial(QTOS((Q)ARG0(arg)),rp); + asir_assert(ARG0(arg),O_N,"fac"); + factorial(QTOS((Q)ARG0(arg)),rp); } -void Plrandom(arg,rp) -NODE arg; -Q *rp; +void Plrandom(NODE arg,Q *rp) { - N r; + N r; - asir_assert(ARG0(arg),O_N,"lrandom"); - randomn(QTOS((Q)ARG0(arg)),&r); - NTOQ(r,1,*rp); + asir_assert(ARG0(arg),O_N,"lrandom"); + randomn(QTOS((Q)ARG0(arg)),&r); + NTOQ(r,1,*rp); } -void Prandom(arg,rp) -NODE arg; -Q *rp; +void Prandom(NODE arg,Q *rp) { - unsigned int r; + unsigned int r; #if 0 #if defined(_PA_RISC1_1) - r = mrand48()&BMASK; + r = mrand48()&BMASK; #else - if ( arg ) - srandom(QTOS((Q)ARG0(arg))); - r = random()&BMASK; + if ( arg ) + srandom(QTOS((Q)ARG0(arg))); + r = random()&BMASK; #endif #endif - if ( arg ) - mt_sgenrand(QTOS((Q)ARG0(arg))); - r = mt_genrand(); - UTOQ(r,*rp); + if ( arg ) + mt_sgenrand(QTOS((Q)ARG0(arg))); + r = mt_genrand(); + UTOQ(r,*rp); } -#if defined(VISUAL) || defined(THINK_C) +#if defined(VISUAL) || defined(__MINGW32__) void srandom(unsigned int); static unsigned int R_Next; @@ -796,178 +839,157 @@ unsigned int random() { return R_Next = (R_Next * 1103515245 + 12345); } -void srandom(s) -unsigned int s; +void srandom(unsigned int s) { - if ( s ) - R_Next = s; + if ( s ) + R_Next = s; else if ( !R_Next ) R_Next = 1; } #endif -void Pprime(arg,rp) -NODE arg; -Q *rp; +void Pprime(NODE arg,Q *rp) { - int i; + int i; - asir_assert(ARG0(arg),O_N,"prime"); - i = QTOS((Q)ARG0(arg)); - if ( i < 0 || i >= 1900 ) - *rp = 0; - else - STOQ(sprime[i],*rp); + asir_assert(ARG0(arg),O_N,"prime"); + i = QTOS((Q)ARG0(arg)); + if ( i < 0 || i >= 1900 ) + *rp = 0; + else + STOQ(sprime[i],*rp); } -void Plprime(arg,rp) -NODE arg; -Q *rp; +void Plprime(NODE arg,Q *rp) { - int i; + int i,p; - asir_assert(ARG0(arg),O_N,"lprime"); - i = QTOS((Q)ARG0(arg)); - if ( i < 0 || i >= 1000 ) - *rp = 0; - else - STOQ(lprime[i],*rp); + asir_assert(ARG0(arg),O_N,"lprime"); + i = QTOS((Q)ARG0(arg)); + if ( i < 0 ) + *rp = 0; + else + p = get_lprime(i); + STOQ(p,*rp); } extern int up_kara_mag, up_tkara_mag, up_fft_mag; -void Pset_upfft(arg,rp) -NODE arg; -Q *rp; +void Pset_upfft(NODE arg,Q *rp) { - if ( arg ) { - asir_assert(ARG0(arg),O_N,"set_upfft"); - up_fft_mag = QTOS((Q)ARG0(arg)); - } - STOQ(up_fft_mag,*rp); + if ( arg ) { + asir_assert(ARG0(arg),O_N,"set_upfft"); + up_fft_mag = QTOS((Q)ARG0(arg)); + } + STOQ(up_fft_mag,*rp); } -void Pset_upkara(arg,rp) -NODE arg; -Q *rp; +void Pset_upkara(NODE arg,Q *rp) { - if ( arg ) { - asir_assert(ARG0(arg),O_N,"set_upkara"); - up_kara_mag = QTOS((Q)ARG0(arg)); - } - STOQ(up_kara_mag,*rp); + if ( arg ) { + asir_assert(ARG0(arg),O_N,"set_upkara"); + up_kara_mag = QTOS((Q)ARG0(arg)); + } + STOQ(up_kara_mag,*rp); } -void Pset_uptkara(arg,rp) -NODE arg; -Q *rp; +void Pset_uptkara(NODE arg,Q *rp) { - if ( arg ) { - asir_assert(ARG0(arg),O_N,"set_uptkara"); - up_tkara_mag = QTOS((Q)ARG0(arg)); - } - STOQ(up_tkara_mag,*rp); + if ( arg ) { + asir_assert(ARG0(arg),O_N,"set_uptkara"); + up_tkara_mag = QTOS((Q)ARG0(arg)); + } + STOQ(up_tkara_mag,*rp); } extern int up2_kara_mag; -void Pset_up2kara(arg,rp) -NODE arg; -Q *rp; +void Pset_up2kara(NODE arg,Q *rp) { - if ( arg ) { - asir_assert(ARG0(arg),O_N,"set_up2kara"); - up2_kara_mag = QTOS((Q)ARG0(arg)); - } - STOQ(up2_kara_mag,*rp); + if ( arg ) { + asir_assert(ARG0(arg),O_N,"set_up2kara"); + up2_kara_mag = QTOS((Q)ARG0(arg)); + } + STOQ(up2_kara_mag,*rp); } -#ifdef HMEXT -void Pigcdbin(arg,rp) -NODE arg; -Obj *rp; +void Pigcdbin(NODE arg,Obj *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"igcd"); - asir_assert(n2,O_N,"igcd"); - if ( !n1 ) - *rp = (Obj)n2; - else if ( !n2 ) - *rp = (Obj)n1; - else { - gcdbinn(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = (Obj)a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"igcd"); + asir_assert(n2,O_N,"igcd"); + if ( !n1 ) + *rp = (Obj)n2; + else if ( !n2 ) + *rp = (Obj)n1; + else { + gcdbinn(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = (Obj)a; + } } -void Pigcdbmod(arg,rp) -NODE arg; -Obj *rp; +void Pigcdbmod(NODE arg,Obj *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"igcdbmod"); - asir_assert(n2,O_N,"igcdbmod"); - if ( !n1 ) - *rp = (Obj)n2; - else if ( !n2 ) - *rp = (Obj)n1; - else { - gcdbmodn(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = (Obj)a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"igcdbmod"); + asir_assert(n2,O_N,"igcdbmod"); + if ( !n1 ) + *rp = (Obj)n2; + else if ( !n2 ) + *rp = (Obj)n1; + else { + gcdbmodn(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = (Obj)a; + } } -void Pigcdacc(arg,rp) -NODE arg; -Obj *rp; +void Pigcdacc(NODE arg,Obj *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"igcdacc"); - asir_assert(n2,O_N,"igcdacc"); - if ( !n1 ) - *rp = (Obj)n2; - else if ( !n2 ) - *rp = (Obj)n1; - else { - gcdaccn(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = (Obj)a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"igcdacc"); + asir_assert(n2,O_N,"igcdacc"); + if ( !n1 ) + *rp = (Obj)n2; + else if ( !n2 ) + *rp = (Obj)n1; + else { + gcdaccn(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = (Obj)a; + } } -void PigcdEuc(arg,rp) -NODE arg; -Obj *rp; +void PigcdEuc(NODE arg,Obj *rp) { - N g; - Q n1,n2,a; - - n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); - asir_assert(n1,O_N,"igcdbmod"); - asir_assert(n2,O_N,"igcdbmod"); - if ( !n1 ) - *rp = (Obj)n2; - else if ( !n2 ) - *rp = (Obj)n1; - else { - gcdEuclidn(NM(n1),NM(n2),&g); - NTOQ(g,1,a); *rp = (Obj)a; - } + N g; + Q n1,n2,a; + + n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg); + asir_assert(n1,O_N,"igcdbmod"); + asir_assert(n2,O_N,"igcdbmod"); + if ( !n1 ) + *rp = (Obj)n2; + else if ( !n2 ) + *rp = (Obj)n1; + else { + gcdEuclidn(NM(n1),NM(n2),&g); + NTOQ(g,1,a); *rp = (Obj)a; + } } extern int igcd_algorithm; - /* == 0 : Euclid, - * == 1 : binary, - * == 2 : bmod, - * >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm, + /* == 0 : Euclid, + * == 1 : binary, + * == 2 : bmod, + * >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm, */ extern int igcd_thre_inidiv; /* @@ -982,111 +1004,52 @@ extern int igcdacc_thre; * > igcdacc_thre, "bmod" reduction is done. */ -void Pigcdcntl(arg,rp) -NODE arg; -Obj *rp; +void Pigcdcntl(NODE arg,Obj *rp) { - Obj p; - Q a; - int k, m; + Obj p; + Q a; + int k, m; - if ( arg ) { - p = (Obj)ARG0(arg); - if ( !p ) { - igcd_algorithm = 0; - *rp = p; - return; - } else if ( OID(p) == O_N ) { - k = QTOS((Q)p); - a = (Q)p; - if ( k >= 0 ) igcd_algorithm = k; - else if ( k == -1 ) { - ret_thre: - k = - igcd_thre_inidiv - igcdacc_thre*10000; - STOQ(k,a); - *rp = (Obj)a; - return; - } else { - if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m; - if ( (m = -k/10000) != 0 ) igcdacc_thre = m; - goto ret_thre; - } - } else if ( OID(p) == O_STR ) { - char *n = BDY((STRING) p); + if ( arg ) { + p = (Obj)ARG0(arg); + if ( !p ) { + igcd_algorithm = 0; + *rp = p; + return; + } else if ( OID(p) == O_N ) { + k = QTOS((Q)p); + a = (Q)p; + if ( k >= 0 ) igcd_algorithm = k; + else if ( k == -1 ) { + ret_thre: + k = - igcd_thre_inidiv - igcdacc_thre*10000; + STOQ(k,a); + *rp = (Obj)a; + return; + } else { + if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m; + if ( (m = -k/10000) != 0 ) igcdacc_thre = m; + goto ret_thre; + } + } else if ( OID(p) == O_STR ) { + char *n = BDY((STRING) p); - if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" ) - || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) ) - k = igcd_algorithm = 1; - else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) ) - igcd_algorithm = 2; - else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" ) - || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) ) - igcd_algorithm = 0; - else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" ) - || !strcmp( n, "gen" ) || !strcmp( n, "Gen" ) - || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) ) - igcd_algorithm = 3; - else error( "igcdhow : invalid algorithm specification" ); - } - } - STOQ(igcd_algorithm,a); - *rp = (Obj)a; - return; + if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" ) + || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) ) + k = igcd_algorithm = 1; + else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) ) + igcd_algorithm = 2; + else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" ) + || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) ) + igcd_algorithm = 0; + else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" ) + || !strcmp( n, "gen" ) || !strcmp( n, "Gen" ) + || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) ) + igcd_algorithm = 3; + else error( "igcdhow : invalid algorithm specification" ); + } + } + STOQ(igcd_algorithm,a); + *rp = (Obj)a; + return; } - - -void Pimaxrsh(arg,rp) -NODE arg; -LIST *rp; -{ - N n1, n2; - Q q; - NODE node1, node2; - int bits; - N maxrshn(); - - q = (Q)ARG0(arg); - asir_assert(q,O_N,"imaxrsh"); - if ( !q ) n1 = n2 = 0; - else { - n1 = maxrshn( NM(q), &bits ); - STON( bits, n2 ); - } - NTOQ( n2, 1, q ); - MKNODE( node2, q, 0 ); - NTOQ( n1, 1, q ); - MKNODE( node1, q, node2 ); - MKLIST( *rp, node1 ); -} -void Pilen(arg,rp) -NODE arg; -Obj *rp; -{ - Q q; - int l; - - q = (Q)ARG0(arg); - asir_assert(q,O_N,"ilenh"); - if ( !q ) l = 0; - else l = PL(NM(q)); - STOQ(l,q); - *rp = (Obj)q; -} - -void Pihex(arg,rp) -NODE arg; -Obj *rp; -{ - Q n1; - - n1 = (Q)ARG0(arg); - asir_assert(n1,O_N,"ihex"); - if ( n1 ) { - int i, l = PL(NM(n1)), *p = BD(NM(n1)); - - for ( i = 0; i < l; i++ ) printf( " 0x%08x", p[i] ); - printf( "\n" ); - } else printf( " 0x%08x\n", 0 ); - *rp = (Obj)n1; -} -#endif /* HMEXT */