/*
* 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/asir2018/builtin/poly.c,v 1.4 2020/10/06 06:31:19 noro Exp $
*/
#include "ca.h"
#include "parse.h"
#include "base.h"
void Pranp();
void Pheadsgn();
void Pmul_trunc(),Pquo_trunc();
void Pumul(),Pumul_ff(),Pusquare(),Pusquare_ff(),Putmul(),Putmul_ff();
void Pkmul(),Pksquare(),Pktmul();
void Pord(), Premove_vars(), Pcoef0(), Pcoef(), Pdeg(), Pmindeg(), Psetmod();
void Pcoef_gf2n();
void getcoef(), getdeglist(), mergedeglist(), change_mvar(), restore_mvar();
void Pp_mag(),Pmaxblen();
void Pmergelist(), Pch_mv(), Pre_mv(), Pdeglist();
void Pptomp(),Pmptop();
void Pptolmp(),Plmptop();
void Pptosfp(),Psfptop(),Psf_galois_action(),Psf_embed(),Psf_find_root();
void Psf_minipoly(),Psf_log(),Psfptopsfp();
void Pptogf2n(),Pgf2ntop(),Pgf2ntovect();
void Pptogfpn(),Pgfpntop();
void Pfind_root_gf2n();
void Pumul_specialmod(),Pusquare_specialmod(),Putmul_specialmod();
void Pureverse(),Putrunc(),Pudecomp(),Purembymul(),Purembymul_precomp();
void Puinvmod(),Purevinvmod();
void Ppwrmod_ff(),Ppwrtab_ff(),Pgeneric_pwrmod_ff();
void Pkpwrmod_lm(),Pkpwrtab_lm(),Pkgeneric_pwrmod_lm();
void Pkmulum();
void Pksquareum();
void Pfmultest();
void Phfmul_lm();
void Plazy_lm();
void Psetmod_ff();
void Psimp_ff();
void Pextdeg_ff();
void Pcharacteristic_ff();
void Pfield_type_ff();
void Pfield_order_ff();
void Prandom_ff();
void Ptracemod_gf2n();
void Psparsemod_gf2n();
void Pmultest_gf2n();
void Psquaretest_gf2n();
void Pinvtest_gf2n();
void Pbininv_gf2n();
void Prinvtest_gf2n();
void Pis_irred_gf2();
void Pis_irred_ddd_gf2();
void Pget_next_fft_prime();
void Puadj_coef();
void Preorder();
void Phomogeneous_part();
void Phomogeneous_deg();
void simp_ff(Obj,Obj *);
void ranp(int,UP *);
void field_order_ff(Z *);
int current_ff;
struct ftab poly_tab[] = {
{"headsgn",Pheadsgn,1},
{"quo_trunc",Pquo_trunc,2},
{"mul_trunc",Pmul_trunc,3},
{"homogeneous_deg",Phomogeneous_deg,-2},
{"homogeneous_part",Phomogeneous_part,-3},
{"reorder",Preorder,3},
{"uadj_coef",Puadj_coef,3},
{"ranp",Pranp,2},
{"p_mag",Pp_mag,1},
{"maxblen",Pmaxblen,1},
{"ord",Pord,-1},
{"remove_vars",Premove_vars,1},
{"delete_vars",Premove_vars,1},
{"coef0",Pcoef0,-3},
{"coef",Pcoef,-3},
{"coef_gf2n",Pcoef_gf2n,2},
{"deg",Pdeg,2},
{"mindeg",Pmindeg,2},
{"setmod",Psetmod,-1},
{"sparsemod_gf2n",Psparsemod_gf2n,-1},
{"setmod_ff",Psetmod_ff,-3},
{"simp_ff",Psimp_ff,1},
{"extdeg_ff",Pextdeg_ff,0},
{"characteristic_ff",Pcharacteristic_ff,0},
{"field_type_ff",Pfield_type_ff,0},
{"field_order_ff",Pfield_order_ff,0},
{"random_ff",Prandom_ff,0},
{"deglist",Pdeglist,2},
{"mergelist",Pmergelist,2},
{"ch_mv",Pch_mv,2},
{"re_mv",Pre_mv,2},
{"ptomp",Pptomp,-2},
{"mptop",Pmptop,1},
{"ptolmp",Pptolmp,1},
{"lmptop",Plmptop,1},
{"sf_galois_action",Psf_galois_action,2},
{"sf_find_root",Psf_find_root,1},
{"sf_minipoly",Psf_minipoly,2},
{"sf_embed",Psf_embed,3},
{"sf_log",Psf_log,1},
{"ptosfp",Pptosfp,1},
{"sfptop",Psfptop,1},
{"sfptopsfp",Psfptopsfp,2},
{"ptogf2n",Pptogf2n,1},
{"gf2ntop",Pgf2ntop,-2},
{"gf2ntovect",Pgf2ntovect,1},
{"ptogfpn",Pptogfpn,1},
{"gfpntop",Pgfpntop,-2},
{"kmul",Pkmul,2},
{"ksquare",Pksquare,1},
{"ktmul",Pktmul,3},
{"umul",Pumul,2},
{"usquare",Pusquare,1},
{"ureverse_inv_as_power_series",Purevinvmod,2},
{"uinv_as_power_series",Puinvmod,2},
{"umul_specialmod",Pumul_specialmod,3},
{"usquare_specialmod",Pusquare_specialmod,2},
{"utmul_specialmod",Putmul_specialmod,4},
{"utmul",Putmul,3},
{"umul_ff",Pumul_ff,2},
{"usquare_ff",Pusquare_ff,1},
{"utmul_ff",Putmul_ff,3},
/* for historical reason */
{"trunc",Putrunc,2},
{"decomp",Pudecomp,2},
{"utrunc",Putrunc,2},
{"udecomp",Pudecomp,2},
{"ureverse",Pureverse,-2},
{"urembymul",Purembymul,2},
{"urembymul_precomp",Purembymul_precomp,3},
{"lazy_lm",Plazy_lm,1},
{"lazy_ff",Plazy_lm,1},
{"pwrmod_ff",Ppwrmod_ff,1},
{"generic_pwrmod_ff",Pgeneric_pwrmod_ff,3},
{"pwrtab_ff",Ppwrtab_ff,2},
{"tracemod_gf2n",Ptracemod_gf2n,3},
{"b_find_root_gf2n",Pfind_root_gf2n,1},
{"is_irred_gf2",Pis_irred_gf2,1},
{"is_irred_ddd_gf2",Pis_irred_ddd_gf2,1},
{"kpwrmod_lm",Pkpwrmod_lm,1},
{"kgeneric_pwrmod_lm",Pkgeneric_pwrmod_lm,3},
{"kpwrtab_lm",Pkpwrtab_lm,2},
{"kmulum",Pkmulum,3},
{"ksquareum",Pksquareum,2},
{"fmultest",Pfmultest,3},
{"hfmul_lm",Phfmul_lm,2},
{"multest_gf2n",Pmultest_gf2n,2},
{"squaretest_gf2n",Psquaretest_gf2n,1},
{"bininv_gf2n",Pbininv_gf2n,2},
{"invtest_gf2n",Pinvtest_gf2n,1},
{"rinvtest_gf2n",Prinvtest_gf2n,0},
{"get_next_fft_prime",Pget_next_fft_prime,2},
{0,0,0},
};
void Pheadsgn(NODE arg,Z *rp)
{
int s;
s = headsgn((P)ARG0(arg));
STOZ(s,*rp);
}
void Pmul_trunc(NODE arg,P *rp)
{
P p1,p2,p,h;
VL vl0,vl1,vl2,tvl,vl;
VN vn;
int i,n;
p1 = (P)ARG0(arg);
p2 = (P)ARG1(arg);
get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&tvl);
p = (P)ARG2(arg);
get_vars((Obj)p,&vl0); mergev(CO,tvl,vl0,&vl);
for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
vn[i].v = tvl->v;
vn[i].n = 0;
}
vn[i].v = 0;
vn[i].n = 0;
for ( h = p, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
for ( ; vn[i].v != VR(h); i++ );
vn[i].n = ZTOS(DEG(DC(h)));
}
mulp_trunc(vl,p1,p2,vn,rp);
}
void Pquo_trunc(NODE arg,P *rp)
{
P p1,p2,p,h;
VL vl0,vl1,vl2,tvl,vl;
VN vn;
int i,n;
p1 = (P)ARG0(arg);
p2 = (P)ARG1(arg);
if ( !p1 )
*rp = 0;
else if ( NUM(p2) )
divsp(CO,p1,p2,rp);
else {
get_vars((Obj)p1,&vl1); get_vars((Obj)p2,&vl2); mergev(CO,vl1,vl2,&vl);
for ( tvl = vl, n = 0; tvl; tvl = NEXT(tvl), n++ );
vn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
for ( i = 0, tvl = vl; i < n; tvl = NEXT(tvl), i++ ) {
vn[i].v = tvl->v;
vn[i].n = 0;
}
vn[i].v = 0;
vn[i].n = 0;
for ( h = p2, i = 0; OID(h) == O_P; h = COEF(DC(h)) ) {
for ( ; vn[i].v != VR(h); i++ );
vn[i].n = ZTOS(DEG(DC(h)));
}
quop_trunc(vl,p1,p2,vn,rp);
}
}
void Phomogeneous_part(NODE arg,P *rp)
{
if ( argc(arg) == 2 )
exthp(CO,(P)ARG0(arg),ZTOS((Q)ARG1(arg)),rp);
else
exthpc_generic(CO,(P)ARG0(arg),ZTOS((Q)ARG2(arg)),
VR((P)ARG1(arg)),rp);
}
void Phomogeneous_deg(NODE arg,Z *rp)
{
int d;
if ( argc(arg) == 1 )
d = homdeg((P)ARG0(arg));
else
d = getchomdeg(VR((P)ARG1(arg)),(P)ARG0(arg));
STOZ(d,*rp);
}
/*
p1 = reorder(p,ovl,nvl) => p1 is 'sorted accoding to nvl.
*/
void Preorder(NODE arg,P *rp)
{
VL ovl,nvl,tvl;
NODE n;
for ( ovl = 0, n = BDY((LIST)ARG1(arg)); n; n = NEXT(n) ) {
if ( !ovl ) {
NEWVL(ovl); tvl = ovl;
} else {
NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
}
VR(tvl) = VR((P)BDY(n));
}
for ( nvl = 0, n = BDY((LIST)ARG2(arg)); n; n = NEXT(n) ) {
if ( !nvl ) {
NEWVL(nvl); tvl = nvl;
} else {
NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
}
VR(tvl) = VR((P)BDY(n));
}
reorderp(nvl,ovl,(P)ARG0(arg),rp);
}
/*
uadj_coef(F,M,M2)
if ( F is a non-negative integer )
return F > M2 ? F-M : M2;
else
F = CN*V^N+...+C0
return uadj_coef(CN,M,M2)*V^N+...+uadj_coef(C0,M,M2);
*/
void Puadj_coef(NODE arg,P *rp)
{
UP f,r;
Z m,m2;
ptoup((P)ARG0(arg),&f);
m = (Z)ARG1(arg);
m2 = (Z)ARG2(arg);
adj_coefup(f,m,m2,&r);
uptop(r,rp);
}
/*
get_next_fft_prime(StartIndex,Bits)
tries to find smallest Index >= StartIndex s.t.
2^(Bits-1)|FFTprime[Index]-1
return [Index,Mod] or 0 (not exist)
*/
void Pget_next_fft_prime(NODE arg,LIST *rp)
{
unsigned int mod,d;
int start,bits,i;
NODE n;
Z q,ind;
start = ZTOS((Q)ARG0(arg));
bits = ZTOS((Q)ARG1(arg));
for ( i = start; ; i++ ) {
get_fft_prime(i,(int *)&mod,(int *)&d);
if ( !mod ) {
*rp = 0; return;
}
if ( bits <= (int)d ) {
UTOZ(mod,q);
UTOZ(i,ind);
n = mknode(2,ind,q);
MKLIST(*rp,n);
return;
}
}
}
void Pranp(NODE arg,P *rp)
{
int n;
UP c;
n = ZTOS((Q)ARG0(arg));
ranp(n,&c);
if ( c ) {
up_var = VR((P)ARG1(arg));
uptop(c,rp);
} else
*rp = 0;
}
void ranp(int n,UP *nr)
{
int i;
unsigned int r;
Z q;
UP c;
*nr = c = UPALLOC(n);
for ( i = 0; i <= n; i++ ) {
r = random();
UTOZ(r,q);
c->c[i] = (Num)q;
}
for ( i = n; i >= 0 && !c->c[i]; i-- );
if ( i >= 0 )
c->d = i;
else
*nr = 0;
}
void Pmaxblen(NODE arg,Z *rp)
{
int l;
l = maxblenp(ARG0(arg));
STOZ(l,*rp);
}
void Pp_mag(NODE arg,Z *rp)
{
int l;
l = p_mag(ARG0(arg));
STOZ(l,*rp);
}
void Pord(NODE arg,LIST *listp)
{
NODE n,tn,p,opt;
char *key;
Obj value;
int overwrite=0;
LIST l;
VL vl,tvl,svl;
P t;
int i,j;
V *va;
V v;
#if 0
printf("LASTCO="); printv(CO,LASTCO->v); printf("\n");
#endif
if ( current_option ) {
for ( opt = current_option; opt; opt = NEXT(opt) ) {
p = BDY((LIST)BDY(opt));
key = BDY((STRING)BDY(p));
value = (Obj)BDY(NEXT(p));
if ( !strcmp(key,"overwrite") && value ) {
overwrite = value ? 1 : 0;
break;
}
}
}
if ( argc(arg) ) {
asir_assert(ARG0(arg),O_LIST,"ord");
for ( vl = 0, i = 0, n = BDY((LIST)ARG0(arg));
n; n = NEXT(n), i++ ) {
if ( !vl ) {
NEWVL(vl); tvl = vl;
} else {
NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
}
if ( !(t = (P)BDY(n)) || (OID(t) != O_P) )
error("ord : invalid argument");
VR(tvl) = VR(t);
}
if ( !overwrite ) {
va = (V *)ALLOCA(i*sizeof(V));
for ( j = 0, svl = vl; j < i; j++, svl = NEXT(svl) )
va[j] = VR(svl);
for ( svl = CO; svl; svl = NEXT(svl) ) {
v = VR(svl);
for ( j = 0; j < i; j++ )
if ( v == va[j] )
break;
if ( j == i ) {
if ( !vl ) {
NEWVL(vl); tvl = vl;
} else {
NEWVL(NEXT(tvl)); tvl = NEXT(tvl);
}
VR(tvl) = v;
}
}
} else {
for ( svl = vl; svl; svl = NEXT(svl) ) {
if ( svl->v->attr == (pointer)V_PF )
((PFINS)svl->v->priv)->pf->ins = 0;
}
}
if ( vl )
NEXT(tvl) = 0;
CO = vl;
update_LASTCO();
}
for ( n = 0, vl = CO; vl; vl = NEXT(vl) ) {
NEXTNODE(n,tn); MKV(VR(vl),t); BDY(tn) = (pointer)t;
}
NEXT(tn) = 0; MKLIST(l,n); *listp = l;
}
void Premove_vars(NODE arg,LIST *listp)
{
NODE l,nd,tnd;
V *v,*va;
int n,na,i,j;
VL vl,vl1;
P t;
LIST list;
asir_assert(ARG0(arg),O_LIST,"remove_vars");
l = BDY((LIST)ARG0(arg)); n = length(l);
v = (V *)ALLOCA(n*sizeof(V));
for ( i = 0; i < n; i++, l = NEXT(l) )
if ( !(t = (P)BDY(l)) || (OID(t) != O_P) )
error("ord : invalid argument");
else v[i] = VR(t);
for ( na = 0, vl = CO; vl; vl = NEXT(vl), na++ );
va = (V *)ALLOCA(na*sizeof(V));
for ( i = 0, vl = CO; i < na; i++, vl = NEXT(vl) ) va[i] = VR(vl);
for ( i = 0; i < na; i++ )
for ( j = 0; j < n; j++ ) if ( va[i] == v[j] ) va[i] = 0;
for ( vl = 0, i = na-1; i >= 0; i-- )
if ( va[i] ) {
NEWVL(vl1); VR(vl1) = va[i]; NEXT(vl1) = vl; vl = vl1;
}
CO = vl;
for ( nd = 0, vl = CO; vl; vl = NEXT(vl) ) {
NEXTNODE(nd,tnd); MKV(VR(vl),t); BDY(tnd) = (pointer)t;
}
if ( nd ) NEXT(tnd) = 0;
MKLIST(list,nd); *listp = list;
}
void Pcoef0(NODE arg,Obj *rp)
{
Obj t,n;
P s;
DCP dc;
int id;
V v;
VL vl;
if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
*rp = 0;
else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
*rp = 0;
else if ( id == O_N )
if ( !n )
*rp = t;
else
*rp = 0;
else {
if ( argc(arg) == 3 ) {
if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
reordvar(CO,v,&vl); reorderp(vl,CO,(P)t,&s);
} else
s = (P)t;
if ( VR(s) != v ) {
if ( n )
*rp = 0;
else
*rp = t;
return;
}
} else
s = (P)t;
for ( dc = DC(s); dc && cmpz(DEG(dc),(Z)n); dc = NEXT(dc) );
if ( dc )
*rp = (Obj)COEF(dc);
else
*rp = 0;
}
}
void Pcoef(NODE arg,Obj *rp)
{
Obj t,n;
P s;
DCP dc;
int id;
V v;
if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
*rp = 0;
else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
*rp = 0;
else if ( id == O_N ) {
if ( !n )
*rp = t;
else
*rp = 0;
} else {
if ( argc(arg) == 3 ) {
if ( (v = VR((P)ARG2(arg))) != VR((P)t) ) {
getcoef(CO,(P)t,v,(Z)n,(P *)rp); return;
} else
s = (P)t;
if ( VR(s) != v ) {
if ( n )
*rp = 0;
else
*rp = t;
return;
}
} else
s = (P)t;
for ( dc = DC(s); dc && cmpz(DEG(dc),(Z)n); dc = NEXT(dc) );
if ( dc )
*rp = (Obj)COEF(dc);
else
*rp = 0;
}
}
void Pcoef_gf2n(NODE arg,Obj *rp)
{
Obj t,n;
int id,d;
UP2 up2;
if ( !(t = (Obj)ARG0(arg)) || ((id = OID(ARG0(arg))) > O_P) )
*rp = 0;
else if ( (n = (Obj)ARG1(arg)) && (OID(n) > O_N) )
*rp = 0;
else if ( id == O_N && NID((Num)t) == N_GF2N ) {
d = ZTOS((Q)n);
up2 = ((GF2N)t)->body;
if ( d > degup2(up2) )
*rp = 0;
else
*rp = (Obj)(up2->b[d/BSH]&(((unsigned long)1)<<(d%BSH))?ONE:0);
} else
*rp = 0;
}
void Pdeg(NODE arg,Z *rp)
{
Obj t,v;
int d;
#if 0
if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
*rp = 0;
else if ( !(v = (Obj)ARG1(arg)) || (VR((P)v) != VR((P)t)) )
*rp = 0;
else
*rp = (Obj)DEG(DC((P)t));
#endif
if ( !(t = (Obj)ARG0(arg)) )
STOZ(-1,*rp);
else if ( OID(t) != O_P ) {
if ( OID(t) == O_N && NID(t) == N_GF2N
&& (v=(Obj)ARG1(arg)) && OID(v)== O_N && NID(v) == N_GF2N ) {
d = degup2(((GF2N)t)->body);
STOZ(d,*rp);
} else
*rp = 0;
} else
degp(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
}
void Pmindeg(NODE arg,Z *rp)
{
Obj t;
if ( !(t = (Obj)ARG0(arg)) || (OID(t) != O_P) )
*rp = 0;
else
getmindeg(VR((P)ARG1(arg)),(P)ARG0(arg),rp);
}
void Psetmod(NODE arg,Z *rp)
{
if ( arg ) {
asir_assert(ARG0(arg),O_N,"setmod");
current_mod = ZTOS((Q)ARG0(arg));
}
STOZ(current_mod,*rp);
}
void Psparsemod_gf2n(NODE arg,Z *rp)
{
int id;
if ( arg && current_mod_gf2n )
current_mod_gf2n->id = ARG0(arg)?1:0;
if ( !current_mod_gf2n )
id = -1;
else
id = current_mod_gf2n->id;
STOZ(id,*rp);
}
void Pmultest_gf2n(NODE arg,GF2N *rp)
{
GF2N a,b,c;
int i;
a = (GF2N)ARG0(arg); b = (GF2N)ARG0(arg);
for ( i = 0; i < 10000; i++ )
mulgf2n(a,b,&c);
*rp = c;
}
void Psquaretest_gf2n(NODE arg,GF2N *rp)
{
GF2N a,c;
int i;
a = (GF2N)ARG0(arg);
for ( i = 0; i < 10000; i++ )
squaregf2n(a,&c);
*rp = c;
}
void Pinvtest_gf2n(NODE arg,GF2N *rp)
{
GF2N a,c;
int i;
a = (GF2N)ARG0(arg);
for ( i = 0; i < 10000; i++ )
invgf2n(a,&c);
*rp = c;
}
void Pbininv_gf2n(NODE arg,GF2N *rp)
{
UP2 a,inv;
int n;
a = ((GF2N)ARG0(arg))->body;
n = ZTOS((Q)ARG1(arg));
type1_bin_invup2(a,n,&inv);
MKGF2N(inv,*rp);
}
void Prinvtest_gf2n(Real *rp)
{
GF2N *a;
GF2N c;
double t0,t1,r;
int i;
double get_clock();
a = (GF2N *)ALLOCA(1000*sizeof(GF2N));
for ( i = 0; i < 1000; i++ ) {
randomgf2n(&a[i]);
}
t0 = get_clock();
for ( i = 0; i < 1000; i++ )
invgf2n(a[i],&c);
t1 = get_clock();
r = (t1-t0)/1000;
MKReal(r,*rp);
}
void Pfind_root_gf2n(NODE arg,GF2N *rp)
{
#if 0
UP p;
ptoup((P)ARG0(arg),&p);
find_root_gf2n(p,rp);
#else
UP2 p;
ptoup2((P)ARG0(arg),&p);
find_root_up2(p,rp);
#endif
}
void Pis_irred_gf2(NODE arg,Z *rp)
{
UP2 t;
ptoup2(ARG0(arg),&t);
*rp = irredcheckup2(t) ? ONE : 0;
}
void Pis_irred_ddd_gf2(NODE arg,Z *rp)
{
UP2 t;
int ret;
ptoup2(ARG0(arg),&t);
ret = irredcheck_dddup2(t);
STOZ(ret,*rp);
}
void Psetmod_ff(NODE arg,Obj *rp)
{
int ac;
int d;
Obj mod,defpoly;
Z n;
UP up;
UP2 up2;
UM dp;
Z q,r;
P p,p1,y;
NODE n0,n1;
LIST list;
ac = argc(arg);
if ( ac == 1 ) {
mod = (Obj)ARG0(arg);
if ( !mod )
current_ff = FF_NOT_SET;
else {
switch ( OID(mod) ) {
case O_N:
current_ff = FF_GFP;
setmod_lm((Z)mod);
break;
case O_P:
current_ff = FF_GF2N;
setmod_gf2n((P)mod); break;
default:
error("setmod_ff : invalid argument");
}
}
} else if ( ac == 2 ) {
if ( OID(ARG0(arg)) == O_N ) {
/* small finite field; primitive root representation */
current_ff = FF_GFS;
setmod_sf(ZTOS((Q)ARG0(arg)),ZTOS((Q)ARG1(arg)));
} else {
mod = (Obj)ARG1(arg);
current_ff = FF_GFPN;
defpoly = (Obj)ARG0(arg);
if ( !mod || !defpoly )
error("setmod_ff : invalid argument");
setmod_lm((Z)mod);
setmod_gfpn((P)defpoly);
}
} else if ( ac == 3 ) {
/* finite extension of a small finite field */
current_ff = FF_GFS;
setmod_sf(ZTOS((Q)ARG0(arg)),ZTOS((Q)ARG1(arg)));
d = ZTOS((Q)ARG2(arg));
generate_defpoly_sfum(d,&dp);
setmod_gfsn(dp);
current_ff = FF_GFSN;
}
switch ( current_ff ) {
case FF_GFP:
getmod_lm(&n); *rp = (Obj)n; break;
case FF_GF2N:
getmod_gf2n(&up2); up2top(up2,&p); *rp = (Obj)p; break;
case FF_GFPN:
getmod_lm(&n);
getmod_gfpn(&up); uptop(up,&p);
MKNODE(n1,n,0); MKNODE(n0,p,n1);
MKLIST(list,n0);
*rp = (Obj)list; break;
case FF_GFS:
case FF_GFSN:
STOZ(current_gfs_p,q);
if ( current_gfs_ext )
enc_to_p(current_gfs_p,current_gfs_iton[1],
VR(current_gfs_ext),&p);
else {
if ( current_gfs_p == 2 )
r = ONE;
else if ( !current_gfs_ntoi )
r = 0;
else
STOZ(current_gfs_iton[1],r);
p = (P)r;
}
switch ( current_ff ) {
case FF_GFS:
n0 = mknode(3,q,current_gfs_ext,p);
break;
case FF_GFSN:
getmod_gfsn(&dp);
makevar("y",&y);
sfumtop(VR(y),dp,&p1);
n0 = mknode(4,q,current_gfs_ext,p,p1);
break;
}
MKLIST(list,n0);
*rp = (Obj)list; break;
default:
*rp = 0; break;
}
}
void Pextdeg_ff(Z *rp)
{
int d;
UP2 up2;
UP up;
UM dp;
switch ( current_ff ) {
case FF_GFP:
*rp = ONE; break;
case FF_GF2N:
getmod_gf2n(&up2); d = degup2(up2); STOZ(d,*rp); break;
case FF_GFPN:
getmod_gfpn(&up); STOZ(up->d,*rp); break;
case FF_GFS:
if ( !current_gfs_ext )
*rp = ONE;
else
*rp = DEG(DC(current_gfs_ext));
break;
case FF_GFSN:
getmod_gfsn(&dp);
STOZ(DEG(dp),*rp);
break;
default:
error("extdeg_ff : current_ff is not set");
}
}
void Pcharacteristic_ff(Z *rp)
{
switch ( current_ff ) {
case FF_GFP:
case FF_GFPN:
getmod_lm(rp); break;
case FF_GF2N:
STOZ(2,*rp); break;
case FF_GFS:
case FF_GFSN:
STOZ(current_gfs_p,*rp); break;
default:
error("characteristic_ff : current_ff is not set");
}
}
void Pfield_type_ff(Z *rp)
{
STOZ(current_ff,*rp);
}
void Pfield_order_ff(Z *rp)
{
field_order_ff(rp);
}
void Prandom_ff(Obj *rp)
{
LM l;
GF2N g;
GFPN p;
GFS s;
GFSN spn;
switch ( current_ff ) {
case FF_GFP:
random_lm(&l); *rp = (Obj)l; break;
case FF_GF2N:
randomgf2n(&g); *rp = (Obj)g; break;
case FF_GFPN:
randomgfpn(&p); *rp = (Obj)p; break;
case FF_GFS:
randomgfs(&s); *rp = (Obj)s; break;
case FF_GFSN:
randomgfsn(&spn); *rp = (Obj)spn; break;
default:
error("random_ff : current_ff is not set");
}
}
void Psimp_ff(NODE arg,Obj *rp)
{
simp_ff((Obj)ARG0(arg),rp);
}
void getcoef(VL vl,P p,V v,Z d,P *r)
{
P s,t,u,a,b,x;
DCP dc;
V w;
if ( !p )
*r = 0;
else if ( NUM(p) )
*r = d ? 0 : p;
else if ( (w=VR(p)) == v ) {
for ( dc = DC(p); dc && cmpz(DEG(dc),d); dc = NEXT(dc) );
*r = dc ? COEF(dc) : 0;
} else {
MKV(w,x);
for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
getcoef(vl,COEF(dc),v,d,&t);
if ( t ) {
pwrp(vl,x,DEG(dc),&u); mulp(vl,t,u,&a);
addp(vl,s,a,&b); s = b;
}
}
*r = s;
}
}
void Pdeglist(NODE arg,LIST *rp)
{
NODE d;
getdeglist((P)ARG0(arg),VR((P)ARG1(arg)),&d);
MKLIST(*rp,d);
}
void Pch_mv(NODE arg,P *rp)
{
change_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
}
void Pre_mv(NODE arg,P *rp)
{
restore_mvar(CO,(P)ARG0(arg),VR((P)ARG1(arg)),rp);
}
void change_mvar(VL vl,P p,V v,P *r)
{
Z d;
DCP dc,dc0;
NODE dl;
if ( !p || NUM(p) || (VR(p) == v) )
*r = p;
else {
getdeglist(p,v,&dl);
for ( dc0 = 0; dl; dl = NEXT(dl) ) {
NEXTDC(dc0,dc); DEG(dc) = d = (Z)BDY(dl);
getcoef(vl,p,v,d,&COEF(dc));
}
NEXT(dc) = 0; MKP(v,dc0,*r);
}
}
void restore_mvar(VL vl,P p,V v,P *r)
{
P s,u,a,b,x;
DCP dc;
if ( !p || NUM(p) || (VR(p) != v) )
*r = p;
else {
MKV(v,x);
for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
pwrp(vl,x,DEG(dc),&u); mulp(vl,COEF(dc),u,&a);
addp(vl,s,a,&b); s = b;
}
*r = s;
}
}
void getdeglist(P p,V v,NODE *d)
{
NODE n,n0,d0,d1,d2;
DCP dc;
if ( !p || NUM(p) ) {
MKNODE(n,0,0); *d = n;
} else if ( VR(p) == v ) {
for ( n0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
NEXTNODE(n0,n); BDY(n) = (pointer)DEG(dc);
}
NEXT(n) = 0; *d = n0;
} else {
for ( dc = DC(p), d0 = 0; dc; dc = NEXT(dc) ) {
getdeglist(COEF(dc),v,&d1); mergedeglist(d0,d1,&d2); d0 = d2;
}
*d = d0;
}
}
void Pmergelist(NODE arg,LIST *rp)
{
NODE n;
asir_assert(ARG0(arg),O_LIST,"mergelist");
asir_assert(ARG1(arg),O_LIST,"mergelist");
mergedeglist(BDY((LIST)ARG0(arg)),BDY((LIST)ARG1(arg)),&n);
MKLIST(*rp,n);
}
void mergedeglist(NODE d0,NODE d1,NODE *dr)
{
NODE t0,t,dt;
Z d;
int c;
if ( !d0 )
*dr = d1;
else {
while ( d1 ) {
dt = d1; d1 = NEXT(d1); d = (Z)BDY(dt);
for ( t0 = 0, t = d0; t; t0 = t, t = NEXT(t) ) {
c = cmpz(d,(Z)BDY(t));
if ( !c )
break;
else if ( c > 0 ) {
if ( !t0 ) {
NEXT(dt) = d0; d0 = dt;
} else {
NEXT(t0) = dt; NEXT(dt) = t;
}
break;
}
}
if ( !t ) {
NEXT(t0) = dt; *dr = d0; return;
}
}
*dr = d0;
}
}
void Pptomp(NODE arg,P *rp)
{
int mod;
if ( argc(arg) == 1 ) {
if ( !current_mod )
error("ptomp : current_mod is not set");
else
mod = current_mod;
} else
mod = ZTOS((Q)ARG1(arg));
ptomp(mod,(P)ARG0(arg),rp);
}
void Pmptop(NODE arg,P *rp)
{
mptop((P)ARG0(arg),rp);
}
void Pptolmp(NODE arg,P *rp)
{
ptolmp((P)ARG0(arg),rp);
}
void Plmptop(NODE arg,P *rp)
{
lmptop((P)ARG0(arg),rp);
}
void Psf_galois_action(NODE arg,P *rp)
{
sf_galois_action(ARG0(arg),ARG1(arg),rp);
}
/*
sf_embed(F,B,PM)
F : an element of GF(pn)
B : the image of the primitive root of GF(pn)
PM : order of GF(pm)
*/
void Psf_embed(NODE arg,P *rp)
{
int k,pm;
/* GF(pn)={0,1,a,a^2,...}->GF(pm)={0,1,b,b^2,...}; a->b^k */
k = CONT((GFS)ARG1(arg));
pm = ZTOS((Q)ARG2(arg));
sf_embed((P)ARG0(arg),k,pm,rp);
}
void Psf_log(NODE arg,Z *rp)
{
int k;
if ( !ARG0(arg) )
error("sf_log : invalid armument");
k = CONT((GFS)ARG0(arg));
STOZ(k,*rp);
}
void Psf_find_root(NODE arg,GFS *rp)
{
P p;
Obj t;
int d;
UM u;
int *root;
p = (P)ARG0(arg);
simp_ff((Obj)p,&t); p = (P)t;
d = getdeg(VR(p),p);
u = W_UMALLOC(d);
ptosfum(p,u);
root = (int *)ALLOCA(d*sizeof(int));
find_rootsf(u,root);
MKGFS(IFTOF(root[0]),*rp);
}
void Psf_minipoly(NODE arg,P *rp)
{
Obj t;
P p1,p2;
int d1,d2;
UM up1,up2,m;
p1 = (P)ARG0(arg); simp_ff((Obj)p1,&t); p1 = (P)t;
p2 = (P)ARG1(arg); simp_ff((Obj)p2,&t); p2 = (P)t;
d1 = getdeg(VR(p1),p1); up1 = W_UMALLOC(d1); ptosfum(p1,up1);
d2 = getdeg(VR(p2),p2); up2 = W_UMALLOC(d2); ptosfum(p2,up2);
m = W_UMALLOC(d2);
minipolysf(up1,up2,m);
sfumtop(VR(p2),m,&p1);
sfptop(p1,rp);
}
void Pptosfp(NODE arg,P *rp)
{
ptosfp(ARG0(arg),rp);
}
void Psfptop(NODE arg,P *rp)
{
sfptop((P)ARG0(arg),rp);
}
void Psfptopsfp(NODE arg,P *rp)
{
sfptopsfp((P)ARG0(arg),VR((P)ARG1(arg)),rp);
}
void Pptogf2n(NODE arg,GF2N *rp)
{
ptogf2n((Obj)ARG0(arg),rp);
}
void Pgf2ntop(NODE arg,P *rp)
{
if ( argc(arg) == 2 )
up2_var = VR((P)ARG1(arg));
gf2ntop((GF2N)ARG0(arg),rp);
}
void Pgf2ntovect(NODE arg,VECT *rp)
{
gf2ntovect((GF2N)ARG0(arg),rp);
}
void Pptogfpn(NODE arg,GFPN *rp)
{
ptogfpn((Obj)ARG0(arg),rp);
}
void Pgfpntop(NODE arg,P *rp)
{
if ( argc(arg) == 2 )
up_var = VR((P)ARG1(arg));
gfpntop((GFPN)ARG0(arg),rp);
}
void Pureverse(NODE arg,P *rp)
{
UP p,r;
ptoup((P)ARG0(arg),&p);
if ( argc(arg) == 1 )
reverseup(p,p->d,&r);
else
reverseup(p,ZTOS((Q)ARG1(arg)),&r);
uptop(r,rp);
}
void Putrunc(NODE arg,P *rp)
{
UP p,r;
ptoup((P)ARG0(arg),&p);
truncup(p,ZTOS((Q)ARG1(arg))+1,&r);
uptop(r,rp);
}
void Pudecomp(NODE arg,LIST *rp)
{
P u,l;
UP p,up,low;
NODE n0,n1;
ptoup((P)ARG0(arg),&p);
decompup(p,ZTOS((Q)ARG1(arg))+1,&low,&up);
uptop(low,&l);
uptop(up,&u);
MKNODE(n1,u,0); MKNODE(n0,l,n1);
MKLIST(*rp,n0);
}
void Purembymul(NODE arg,P *rp)
{
UP p1,p2,r;
if ( !ARG0(arg) || !ARG1(arg) )
*rp = 0;
else {
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
rembymulup(p1,p2,&r);
uptop(r,rp);
}
}
/*
* d1 = deg(p1), d2 = deg(p2)
* d1 <= 2*d2-1
* p2*inv = 1 mod x^d2
*/
void Purembymul_precomp(NODE arg,P *rp)
{
UP p1,p2,inv,r;
if ( !ARG0(arg) || !ARG1(arg) )
*rp = 0;
else {
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
ptoup((P)ARG2(arg),&inv);
if ( p1->d >= 2*p2->d ) {
error("urembymul_precomp : degree of 1st arg is too large");
/* fprintf(stderr,"urembymul_precomp : degree of 1st arg is too large"); */
remup(p1,p2,&r);
} else
hybrid_rembymulup_special(current_ff,p1,p2,inv,&r);
uptop(r,rp);
}
}
void Puinvmod(NODE arg,P *rp)
{
UP p,r;
ptoup((P)ARG0(arg),&p);
invmodup(p,ZTOS((Q)ARG1(arg)),&r);
uptop(r,rp);
}
void Purevinvmod(NODE arg,P *rp)
{
UP p,pr,r;
ptoup((P)ARG0(arg),&p);
reverseup(p,p->d,&pr);
invmodup(pr,ZTOS((Q)ARG1(arg)),&r);
uptop(r,rp);
}
void Ppwrmod_ff(NODE arg,P *rp)
{
UP p1,p2;
ptoup((P)ARG0(arg),&p1);
switch ( current_ff ) {
case FF_GFP:
/* XXX : hybrid version may not be useful ... */
#if 1
hybrid_powermodup(p1,&p2); break;
#else
powermodup(p1,&p2); break;
#endif
case FF_GF2N:
powermodup_gf2n(p1,&p2); break;
case FF_GFPN:
case FF_GFS:
case FF_GFSN:
powermodup(p1,&p2); break;
default:
error("pwrmod_ff : current_ff is not set");
}
uptop(p2,rp);
}
void Pgeneric_pwrmod_ff(NODE arg,P *rp)
{
UP g,f,r;
ptoup((P)ARG0(arg),&g);
ptoup((P)ARG1(arg),&f);
switch ( current_ff ) {
case FF_GFP:
hybrid_generic_powermodup(g,f,(Z)ARG2(arg),&r); break;
case FF_GF2N:
generic_powermodup_gf2n(g,f,(Z)ARG2(arg),&r); break;
case FF_GFPN:
case FF_GFS:
case FF_GFSN:
generic_powermodup(g,f,(Z)ARG2(arg),&r); break;
default:
error("generic_pwrmod_ff : current_ff is not set");
}
uptop(r,rp);
}
void Ppwrtab_ff(NODE arg,VECT *rp)
{
UP f,xp;
UP *tab;
VECT r;
int i,d;
ptoup((P)ARG0(arg),&f);
ptoup((P)ARG1(arg),&xp);
d = f->d;
tab = (UP *)ALLOCA(d*sizeof(UP));
switch ( current_ff ) {
case FF_GFP:
hybrid_powertabup(f,xp,tab); break;
case FF_GF2N:
powertabup_gf2n(f,xp,tab); break;
case FF_GFPN:
case FF_GFS:
case FF_GFSN:
powertabup(f,xp,tab); break;
default:
error("pwrtab_ff : current_ff is not set");
}
MKVECT(r,d); *rp = r;
for ( i = 0; i < d; i++ )
uptop(tab[i],(P *)&BDY(r)[i]);
}
void Pkpwrmod_lm(NODE arg,P *rp)
{
UP p1,p2;
ptoup((P)ARG0(arg),&p1);
powermodup(p1,&p2);
uptop(p2,rp);
}
void Pkgeneric_pwrmod_lm(NODE arg,P *rp)
{
UP g,f,r;
ptoup((P)ARG0(arg),&g);
ptoup((P)ARG1(arg),&f);
generic_powermodup(g,f,(Z)ARG2(arg),&r);
uptop(r,rp);
}
void Pkpwrtab_lm(NODE arg,VECT *rp)
{
UP f,xp;
UP *tab;
VECT r;
int i,d;
ptoup((P)ARG0(arg),&f);
ptoup((P)ARG1(arg),&xp);
d = f->d;
tab = (UP *)ALLOCA(d*sizeof(UP));
powertabup(f,xp,tab);
MKVECT(r,d); *rp = r;
for ( i = 0; i < d; i++ )
uptop(tab[i],(P *)&BDY(r)[i]);
}
void Plazy_lm(NODE arg,Q *rp)
{
lm_lazy = ZTOS((Q)ARG0(arg));
*rp = 0;
}
void Pkmul(NODE arg,P *rp)
{
P n1,n2;
n1 = (P)ARG0(arg); n2 = (P)ARG1(arg);
asir_assert(n1,O_P,"kmul");
asir_assert(n2,O_P,"kmul");
kmulp(CO,n1,n2,rp);
}
void Pksquare(NODE arg,P *rp)
{
P n1;
n1 = (P)ARG0(arg);
asir_assert(n1,O_P,"ksquare");
ksquarep(CO,n1,rp);
}
void Pktmul(NODE arg,P *rp)
{
UP p1,p2,r;
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
tkmulup(p1,p2,ZTOS((Q)ARG2(arg))+1,&r);
uptop(r,rp);
}
void Pumul(NODE arg,P *rp)
{
P a1,a2;
UP p1,p2,r;
a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
mulp(CO,a1,a2,rp);
else {
if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
error("umul : invalid argument");
ptoup(a1,&p1);
ptoup(a2,&p2);
hybrid_mulup(0,p1,p2,&r);
uptop(r,rp);
}
}
void Pusquare(NODE arg,P *rp)
{
UP p1,r;
ptoup((P)ARG0(arg),&p1);
hybrid_squareup(0,p1,&r);
uptop(r,rp);
}
void Putmul(NODE arg,P *rp)
{
UP p1,p2,r;
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
hybrid_tmulup(0,p1,p2,ZTOS((Q)ARG2(arg))+1,&r);
uptop(r,rp);
}
void Pumul_ff(NODE arg,Obj *rp)
{
P a1,a2;
UP p1,p2,r;
P p;
a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
ptoup(a1,&p1);
ptoup(a2,&p2);
hybrid_mulup(current_ff,p1,p2,&r);
uptop(r,&p);
simp_ff((Obj)p,rp);
}
void Pusquare_ff(NODE arg,Obj *rp)
{
UP p1,r;
P p;
ptoup((P)ARG0(arg),&p1);
hybrid_squareup(current_ff,p1,&r);
uptop(r,&p);
simp_ff((Obj)p,rp);
}
void Putmul_ff(NODE arg,Obj *rp)
{
UP p1,p2,r;
P p;
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
hybrid_tmulup(current_ff,p1,p2,ZTOS((Q)ARG2(arg))+1,&r);
uptop(r,&p);
simp_ff((Obj)p,rp);
}
void Phfmul_lm(NODE arg,P *rp)
{
UP p1,p2,r;
UP hi,lo,mid,t,s,p10,p11,p20,p21;
unsigned int m,d;
int i;
ptoup((P)ARG0(arg),&p1);
ptoup((P)ARG1(arg),&p2);
d = MAX(p1->d,p2->d);
for ( m = 1; m < d; m <<= 1 );
if ( m > d )
m >>= 1;
if ( d-m < 10000 ) {
decompup(p1,m,&p10,&p11); /* p1 = p11*x^m+p10 */
decompup(p2,m,&p20,&p21); /* p2 = p21*x^m+p20 */
fft_mulup_lm(p10,p20,&lo);
kmulup(p11,p21,&hi);
kmulup(p11,p20,&t); kmulup(p10,p21,&s); addup(t,s,&mid);
r = UPALLOC(2*d);
bzero((char *)COEF(r),(2*d+1)*sizeof(Num));
if ( lo )
bcopy((char *)COEF(lo),(char *)COEF(r),(DEG(lo)+1)*sizeof(Num));
if ( hi )
bcopy((char *)COEF(hi),(char *)(COEF(r)+2*m),(DEG(hi)+1)*sizeof(Num));
for ( i = 2*d; i >= 0 && !COEF(r)[i]; i-- );
if ( i < 0 )
r = 0;
else {
DEG(r) = i;
t = UPALLOC(DEG(mid)+m); DEG(t) = DEG(mid)+m;
bzero((char *)COEF(t),(DEG(mid)+m+1)*sizeof(Num));
bcopy((char *)COEF(mid),(char *)(COEF(t)+m),(DEG(mid)+1)*sizeof(Num));
addup(t,r,&s);
r = s;
}
} else
fft_mulup_lm(p1,p2,&r);
uptop(r,rp);
}
void Pfmultest(NODE arg,LIST *rp)
{
P p1,p2,r;
int d1,d2;
UM w1,w2,wr;
unsigned int *f1,*f2,*fr,*w;
int index,mod,root,d,maxint,i;
int cond;
Z prime;
NODE n0,n1;
p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); index = ZTOS((Q)ARG2(arg));
FFT_primes(index,&mod,&root,&d);
maxint = 1<<d;
d1 = UDEG(p1); d2 = UDEG(p2);
if ( maxint < d1+d2+1 )
*rp = 0;
else {
w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
wr = W_UMALLOC(d1+d2);
ptoum(mod,p1,w1); ptoum(mod,p2,w2);
f1 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
f2 = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
fr = (unsigned int *)ALLOCA(maxint*sizeof(unsigned int));
w = (unsigned int *)ALLOCA(12*maxint*sizeof(unsigned int));
for ( i = 0; i <= d1; i++ )
f1[i] = (unsigned int)w1->c[i];
for ( i = 0; i <= d2; i++ )
f2[i] = (unsigned int)w2->c[i];
cond = FFT_pol_product(d1,f1,d2,f2,fr,index,w);
if ( cond )
error("fmultest : ???");
wr->d = d1+d2;
for ( i = 0; i <= d1+d2; i++ )
wr->c[i] = (unsigned int)fr[i];
umtop(VR(p1),wr,&r);
STOZ(mod,prime);
MKNODE(n1,prime,0);
MKNODE(n0,r,n1);
MKLIST(*rp,n0);
}
}
void Pkmulum(NODE arg,P *rp)
{
P p1,p2;
int d1,d2,mod;
UM w1,w2,wr;
p1 = (P)ARG0(arg); p2 = (P)ARG1(arg); mod = ZTOS((Q)ARG2(arg));
d1 = UDEG(p1); d2 = UDEG(p2);
w1 = W_UMALLOC(d1); w2 = W_UMALLOC(d2);
wr = W_UMALLOC(d1+d2);
ptoum(mod,p1,w1); ptoum(mod,p2,w2);
kmulum(mod,w1,w2,wr);
umtop(VR(p1),wr,rp);
}
void Pksquareum(NODE arg,P *rp)
{
P p1;
int d1,mod;
UM w1,wr;
p1 = (P)ARG0(arg); mod = ZTOS((Q)ARG1(arg));
d1 = UDEG(p1);
w1 = W_UMALLOC(d1);
wr = W_UMALLOC(2*d1);
ptoum(mod,p1,w1);
kmulum(mod,w1,w1,wr);
umtop(VR(p1),wr,rp);
}
void Ptracemod_gf2n(NODE arg,P *rp)
{
UP g,f,r;
ptoup((P)ARG0(arg),&g);
ptoup((P)ARG1(arg),&f);
tracemodup_gf2n(g,f,(Z)ARG2(arg),&r);
uptop(r,rp);
}
void Pumul_specialmod(NODE arg,P *rp)
{
P a1,a2;
UP p1,p2,r;
int i,nmod;
int *modind;
NODE t,n;
a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
mulp(CO,a1,a2,rp);
else {
if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
error("umul_specialmod : invalid argument");
ptoup(a1,&p1);
ptoup(a2,&p2);
n = BDY((LIST)ARG2(arg));
nmod = length(n);
modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
modind[i] = ZTOS((Q)BDY(t));
fft_mulup_specialmod_main(p1,p2,0,modind,nmod,&r);
uptop(r,rp);
}
}
void Pusquare_specialmod(NODE arg,P *rp)
{
P a1;
UP p1,r;
int i,nmod;
int *modind;
NODE t,n;
a1 = (P)ARG0(arg);
if ( !a1 || NUM(a1) )
mulp(CO,a1,a1,rp);
else {
if ( !uzpcheck((Obj)a1) )
error("usquare_specialmod : invalid argument");
ptoup(a1,&p1);
n = BDY((LIST)ARG1(arg));
nmod = length(n);
modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
modind[i] = ZTOS((Q)BDY(t));
fft_mulup_specialmod_main(p1,p1,0,modind,nmod,&r);
uptop(r,rp);
}
}
void Putmul_specialmod(NODE arg,P *rp)
{
P a1,a2;
UP p1,p2,r;
int i,nmod;
int *modind;
NODE t,n;
a1 = (P)ARG0(arg); a2 = (P)ARG1(arg);
if ( !a1 || !a2 || NUM(a1) || NUM(a2) )
mulp(CO,a1,a2,rp);
else {
if ( !uzpcheck((Obj)a1) || !uzpcheck((Obj)a2) || VR(a1) != VR(a2) )
error("utmul_specialmod : invalid argument");
ptoup(a1,&p1);
ptoup(a2,&p2);
n = BDY((LIST)ARG3(arg));
nmod = length(n);
modind = (int *)MALLOC_ATOMIC(nmod*sizeof(int));
for ( i = 0, t = n; i < nmod; i++, t = NEXT(t) )
modind[i] = ZTOS((Q)BDY(t));
fft_mulup_specialmod_main(p1,p2,ZTOS((Q)ARG2(arg))+1,modind,nmod,&r);
uptop(r,rp);
}
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to poly.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin