version 1.1, 1999/12/03 07:39:06 |
version 1.9, 2018/03/29 01:32:50 |
|
|
/* $OpenXM: OpenXM/src/asir99/asm/ddM.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/asm/ddM.c,v 1.8 2009/03/02 19:01:43 ohara Exp $ |
|
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
#include "inline.h" |
#include "inline.h" |
|
|
void ksquareummain(int,UM,UM); |
|
void kmulummain(int,UM,UM,UM); |
|
void c_copyum(UM,int,int *); |
|
void copyum(UM,UM); |
|
void extractum(UM,int,int,UM); |
|
void ksquareum(int,UM,UM); |
|
void kmulum(int,UM,UM,UM); |
|
|
|
/* |
/* |
* mod is declared as 'int', because several xxxum functions contains signed |
* mod is declared as 'int', because several xxxum functions contains signed |
* integer addition/subtraction. So mod should be less than 2^31. |
* integer addition/subtraction. So mod should be less than 2^31. |
*/ |
*/ |
|
|
void mulum(mod,p1,p2,pr) |
void mulum(int mod,UM p1,UM p2,UM pr) |
int mod; |
|
UM p1,p2,pr; |
|
{ |
{ |
int *pc1,*pcr; |
int *pc1,*pcr; |
int *c1,*c2,*cr; |
int *c1,*c2,*cr; |
unsigned int mul; |
unsigned int mul; |
int i,j,d1,d2; |
int i,j,d1,d2; |
|
|
if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) { |
if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) { |
DEG(pr) = -1; |
DEG(pr) = -1; |
return; |
return; |
} |
} |
c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr); |
c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr); |
bzero((char *)cr,(int)((d1+d2+1)*sizeof(int))); |
bzero((char *)cr,(int)((d1+d2+1)*sizeof(int))); |
for ( i = 0; i <= d2; i++, cr++ ) |
for ( i = 0; i <= d2; i++, cr++ ) |
if ( mul = *c2++ ) |
if ( mul = *c2++ ) |
for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ ) { |
for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ ) { |
DMAR(*pc1,mul,*pcr,mod,*pcr) |
DMAR(*pc1,mul,*pcr,mod,*pcr) |
} |
} |
DEG(pr) = d1 + d2; |
DEG(pr) = d1 + d2; |
} |
} |
|
|
void mulsum(mod,p,n,pr) |
void mulsum(int mod,UM p,int n,UM pr) |
int mod,n; |
|
UM p,pr; |
|
{ |
{ |
int *sp,*dp; |
int *sp,*dp; |
int i; |
int i; |
|
|
for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i; |
for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i; |
i >= 0; i--, dp--, sp-- ) { |
i >= 0; i--, dp--, sp-- ) { |
DMAR(*sp,n,0,mod,*dp) |
DMAR(*sp,n,0,mod,*dp) |
} |
} |
} |
} |
|
|
int divum(mod,p1,p2,pq) |
int divum(int mod,UM p1,UM p2,UM pq) |
int mod; |
|
UM p1,p2,pq; |
|
{ |
{ |
int *pc1,*pct; |
int *pc1,*pct; |
int *c1,*c2,*ct; |
int *c1,*c2,*ct; |
unsigned int inv,hd,tmp; |
unsigned int inv,hd,tmp; |
int i,j, d1,d2,dd; |
int i,j, d1,d2,dd; |
|
|
if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) { |
if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) { |
DEG(pq) = -1; |
DEG(pq) = -1; |
return d1; |
return d1; |
} |
} |
c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2; |
c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2; |
if ( ( hd = c2[d2] ) != 1 ) { |
if ( ( hd = c2[d2] ) != 1 ) { |
inv = invm(hd,mod); |
inv = invm(hd,mod); |
for ( pc1 = c2 + d2; pc1 >= c2; pc1-- ) { |
for ( pc1 = c2 + d2; pc1 >= c2; pc1-- ) { |
DMAR(*pc1,inv,0,mod,*pc1) |
DMAR(*pc1,inv,0,mod,*pc1) |
} |
} |
} else |
} else |
inv = 1; |
inv = 1; |
for ( i = dd, ct = c1+d1; i >= 0; i-- ) |
for ( i = dd, ct = c1+d1; i >= 0; i-- ) |
if ( tmp = *ct-- ) { |
if ( tmp = *ct-- ) { |
tmp = mod - tmp; |
tmp = mod - tmp; |
for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- ) { |
for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- ) { |
DMAR(*pc1,tmp,*pct,mod,*pct) |
DMAR(*pc1,tmp,*pct,mod,*pct) |
} |
} |
} |
} |
if ( inv != 1 ) { |
if ( inv != 1 ) { |
for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ ) { |
for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ ) { |
DMAR(*pc1,inv,0,mod,*pc1) |
DMAR(*pc1,inv,0,mod,*pc1) |
} |
} |
for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ ) { |
for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ ) { |
DMAR(*pc1,inv,0,mod,*pc1) |
DMAR(*pc1,inv,0,mod,*pc1) |
} |
} |
} |
} |
for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- ); |
for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- ); |
for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- ) |
for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- ) |
*pct-- = *pc1--; |
*pct-- = *pc1--; |
return i; |
return i; |
} |
} |
|
|
void diffum(mod,f,fd) |
void diffum(int mod,UM f,UM fd) |
int mod; |
|
UM f,fd; |
|
{ |
{ |
int *dp,*sp; |
int *dp,*sp; |
int i; |
int i; |
UL ltmp; |
|
|
|
for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i; |
for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i; |
i >= 1; i--, dp--, sp-- ) { |
i >= 1; i--, dp--, sp-- ) { |
DMAR(*sp,i,0,mod,*dp) |
DMAR(*sp,i,0,mod,*dp) |
} |
} |
degum(fd,DEG(f) - 1); |
degum(fd,DEG(f) - 1); |
} |
} |
|
|
unsigned int pwrm(mod,a,n) |
unsigned int pwrm(int mod,int a,int n) |
int mod,a; |
|
int n; |
|
{ |
{ |
unsigned int s,t; |
unsigned int s,t; |
|
|
if ( !n ) |
if ( !n ) |
return 1; |
return 1; |
else if ( n == 1 ) |
else if ( n == 1 ) |
return a; |
return a; |
else { |
else { |
t = pwrm(mod,a,n/2); |
t = pwrm(mod,a,n/2); |
DMAR(t,t,0,mod,s) |
DMAR(t,t,0,mod,s) |
if ( n % 2 ) { |
if ( n % 2 ) { |
DMAR(s,a,0,mod,t) |
DMAR(s,a,0,mod,t) |
return t; |
return t; |
} else |
} else |
return s; |
return s; |
} |
} |
} |
} |
|
|
unsigned int invm(s,mod) |
unsigned int invm(unsigned int s,int mod) |
unsigned int s; |
|
int mod; |
|
{ |
{ |
unsigned int r,a2,q; |
unsigned int r,a2,q; |
unsigned int f1,f2,a1; |
unsigned int f1,f2,a1; |
|
|
for ( f1 = s, f2 = mod, a1 = 1, a2 = 0; ; ) { |
for ( f1 = s, f2 = mod, a1 = 1, a2 = 0; ; ) { |
q = f1/f2; r = f1 - f2*q; f1 = f2; |
q = f1/f2; r = f1 - f2*q; f1 = f2; |
if ( !(f2 = r) ) |
if ( !(f2 = r) ) |
break; |
break; |
DMAR(a2,q,0,mod,r) |
DMAR(a2,q,0,mod,r) |
/* r = ( a1 - r + mod ) % mod; */ |
/* r = ( a1 - r + mod ) % mod; */ |
if ( a1 >= r ) |
if ( a1 >= r ) |
r = a1 - r; |
r = a1 - r; |
else { |
else { |
r = (mod - r) + a1; |
r = (mod - r) + a1; |
} |
} |
a1 = a2; a2 = r; |
a1 = a2; a2 = r; |
} |
} |
/* return( ( a2 >= 0 ? a2 : a2 + mod ) ); */ |
/* return( ( a2 >= 0 ? a2 : a2 + mod ) ); */ |
return a2; |
return a2; |
} |
} |
|
|
unsigned int rem(n,m) |
unsigned int rem(N n,int m) |
N n; |
|
unsigned int m; |
|
{ |
{ |
unsigned int *x; |
unsigned int *x; |
unsigned int t,r; |
unsigned int t,r; |
int i; |
int i; |
|
|
if ( !n ) |
if ( !n ) |
return 0; |
return 0; |
for ( i = PL(n)-1, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) { |
for ( i = PL(n)-1, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) { |
#if defined(sparc) |
#if defined(sparc) && !defined(__sparcv9) |
r = dsar(m,r,*x); |
r = dsar(m,r,*x); |
#else |
#else |
DSAB(m,r,*x,t,r) |
DSAB(m,r,*x,t,r) |
#endif |
#endif |
} |
} |
return r; |
return r; |
} |
} |
|
|
#ifndef sparc |
#if !defined(sparc) || defined(__sparcv9) |
void addpadic(mod,n,n1,n2) |
void addpadic(int mod,int n,unsigned int *n1,unsigned int *n2) |
int mod; |
|
int n; |
|
unsigned int *n1,*n2; |
|
{ |
{ |
unsigned int carry,tmp; |
unsigned int carry,tmp; |
int i; |
int i; |
|
|
for ( i = 0, carry = 0; i < n; i++ ) { |
for ( i = 0, carry = 0; i < n; i++ ) { |
tmp = *n1++ + *n2 + carry; |
tmp = *n1++ + *n2 + carry; |
DQR(tmp,mod,carry,*n2++) |
DQR(tmp,mod,carry,*n2++) |
/* |
/* |
carry = tmp / mod; |
carry = tmp / mod; |
*n2++ = tmp - ( carry * mod ); |
*n2++ = tmp - ( carry * mod ); |
*/ |
*/ |
} |
} |
} |
} |
#endif |
#endif |
|
|
void mulpadic(mod,n,n1,n2,nr) |
void mulpadic(int mod,int n,unsigned int *n1,unsigned int *n2,unsigned int *nr) |
int mod; |
|
int n; |
|
unsigned int *n1; |
|
unsigned int *n2,*nr; |
|
{ |
{ |
unsigned int *pn1,*pnr; |
unsigned int *pn1,*pnr; |
unsigned int carry,mul; |
unsigned int carry,mul; |
int i,j; |
int i,j; |
|
|
bzero((char *)nr,(int)(n*sizeof(int))); |
bzero((char *)nr,(int)(n*sizeof(int))); |
for ( j = 0; j < n; j++, n2++, nr++ ) |
for ( j = 0; j < n; j++, n2++, nr++ ) |
if ( mul = *n2 ) |
if ( mul = *n2 ) |
for ( i = j, carry = 0, pn1 = n1, pnr = nr; |
for ( i = j, carry = 0, pn1 = n1, pnr = nr; |
i < n; i++, pn1++, pnr++ ) { |
i < n; i++, pn1++, pnr++ ) { |
carry += *pnr; |
carry += *pnr; |
DMAB(mod,*pn1,mul,carry,carry,*pnr) |
DMAB(mod,*pn1,mul,carry,carry,*pnr) |
} |
} |
} |
} |
|
|
extern up_kara_mag; |
extern up_kara_mag; |
|
|
void kmulum(mod,n1,n2,nr) |
void kmulum(int mod,UM n1,UM n2,UM nr) |
UM n1,n2,nr; |
|
{ |
{ |
UM n,t,s,m,carry; |
UM n,t,s,m,carry; |
int d,d1,d2,len,i,l; |
int d,d1,d2,len,i,l; |
unsigned int *r,*r0; |
unsigned int *r,*r0; |
|
|
if ( !n1 || !n2 ) { |
if ( !n1 || !n2 ) { |
nr->d = -1; return; |
nr->d = -1; return; |
} |
} |
d1 = DEG(n1)+1; d2 = DEG(n2)+1; |
d1 = DEG(n1)+1; d2 = DEG(n2)+1; |
if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) { |
if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) { |
mulum(mod,n1,n2,nr); return; |
mulum(mod,n1,n2,nr); return; |
} |
} |
if ( d1 < d2 ) { |
if ( d1 < d2 ) { |
n = n1; n1 = n2; n2 = n; |
n = n1; n1 = n2; n2 = n; |
d = d1; d1 = d2; d2 = d; |
d = d1; d1 = d2; d2 = d; |
} |
} |
if ( d2 > (d1+1)/2 ) { |
if ( d2 > (d1+1)/2 ) { |
kmulummain(mod,n1,n2,nr); return; |
kmulummain(mod,n1,n2,nr); return; |
} |
} |
d = (d1/d2)+((d1%d2)!=0); |
d = (d1/d2)+((d1%d2)!=0); |
len = (d+1)*d2; |
len = (d+1)*d2; |
r0 = (unsigned int *)ALLOCA(len*sizeof(int)); |
r0 = (unsigned int *)ALLOCA(len*sizeof(int)); |
bzero((char *)r0,len*sizeof(int)); |
bzero((char *)r0,len*sizeof(int)); |
m = W_UMALLOC(d2+1); |
m = W_UMALLOC(d2+1); |
carry = W_UMALLOC(d2+1); |
carry = W_UMALLOC(d2+1); |
t = W_UMALLOC(d1+d2+1); |
t = W_UMALLOC(d1+d2+1); |
s = W_UMALLOC(d1+d2+1); |
s = W_UMALLOC(d1+d2+1); |
for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) { |
for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) { |
extractum(n1,i*d2,d2,m); |
extractum(n1,i*d2,d2,m); |
if ( m ) { |
if ( m ) { |
kmulum(mod,m,n2,t); |
kmulum(mod,m,n2,t); |
addum(mod,t,carry,s); |
addum(mod,t,carry,s); |
c_copyum(s,d2,r); |
c_copyum(s,d2,r); |
extractum(s,d2,d2,carry); |
extractum(s,d2,d2,carry); |
} else { |
} else { |
c_copyum(carry,d2,r); |
c_copyum(carry,d2,r); |
carry = 0; |
carry = 0; |
} |
} |
} |
} |
c_copyum(carry,d2,r); |
c_copyum(carry,d2,r); |
for ( l = len - 1; !r0[l]; l-- ); |
for ( l = len - 1; !r0[l]; l-- ); |
l++; |
l++; |
DEG(nr) = l-1; |
DEG(nr) = l-1; |
bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int)); |
bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int)); |
} |
} |
|
|
void ksquareum(mod,n1,nr) |
void ksquareum(int mod,UM n1,UM nr) |
int mod; |
|
UM n1,nr; |
|
{ |
{ |
int d1; |
int d1; |
|
|
if ( !n1 ) { |
if ( !n1 ) { |
nr->d = -1; return; |
nr->d = -1; return; |
} |
} |
d1 = DEG(n1)+1; |
d1 = DEG(n1)+1; |
if ( (d1 < up_kara_mag) ) { |
if ( (d1 < up_kara_mag) ) { |
pwrum(mod,n1,2,nr); return; |
pwrum(mod,n1,2,nr); return; |
} |
} |
ksquareummain(mod,n1,nr); |
ksquareummain(mod,n1,nr); |
} |
} |
|
|
void extractum(n,index,len,nr) |
void extractum(UM n,int index,int len,UM nr) |
UM n; |
|
int index,len; |
|
UM nr; |
|
{ |
{ |
int *m; |
int *m; |
int l; |
int l; |
|
|
if ( !n ) { |
if ( !n ) { |
nr->d = -1; return; |
nr->d = -1; return; |
} |
} |
m = COEF(n)+index; |
m = COEF(n)+index; |
if ( (l = (DEG(n)+1)-index) >= len ) { |
if ( (l = (DEG(n)+1)-index) >= len ) { |
for ( l = len - 1; (l >= 0) && !m[l]; l-- ); |
for ( l = len - 1; (l >= 0) && !m[l]; l-- ); |
l++; |
l++; |
} |
} |
if ( l <= 0 ) { |
if ( l <= 0 ) { |
nr->d = -1; return; |
nr->d = -1; return; |
} else { |
} else { |
DEG(nr) = l-1; |
DEG(nr) = l-1; |
bcopy((char *)m,(char *)COEF(nr),l*sizeof(Q)); |
bcopy((char *)m,(char *)COEF(nr),l*sizeof(Q)); |
} |
} |
} |
} |
|
|
void copyum(n1,n2) |
void copyum(UM n1,UM n2) |
UM n1,n2; |
|
{ |
{ |
n2->d = n1->d; |
n2->d = n1->d; |
bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(int)); |
bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(int)); |
} |
} |
|
|
void c_copyum(n,len,p) |
void c_copyum(UM n,int len,int *p) |
UM n; |
|
int len; |
|
int *p; |
|
{ |
{ |
if ( n ) |
if ( n ) |
bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(int)); |
bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(int)); |
} |
} |
|
|
void kmulummain(mod,n1,n2,nr) |
void kmulummain(int mod,UM n1,UM n2,UM nr) |
int mod; |
|
UM n1,n2,nr; |
|
{ |
{ |
int d1,d2,h,len; |
int d1,d2,h,len; |
UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; |
UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2; |
|
|
d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2; |
d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2; |
n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1); |
n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1); |
n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1); |
n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1); |
lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1); |
lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1); |
mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1); |
mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1); |
mid = W_UMALLOC(d1+d2+1); |
mid = W_UMALLOC(d1+d2+1); |
s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1); |
s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1); |
extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi); |
extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi); |
extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi); |
extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi); |
kmulum(mod,n1hi,n2hi,hi); kmulum(mod,n1lo,n2lo,lo); |
kmulum(mod,n1hi,n2hi,hi); kmulum(mod,n1lo,n2lo,lo); |
len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1; |
len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1; |
bzero((char *)COEF(t1),len*sizeof(int)); |
bzero((char *)COEF(t1),len*sizeof(int)); |
if ( lo ) |
if ( lo ) |
bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int)); |
bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int)); |
if ( hi ) |
if ( hi ) |
bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int)); |
bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int)); |
|
|
addum(mod,hi,lo,mid1); |
addum(mod,hi,lo,mid1); |
subum(mod,n1hi,n1lo,s1); subum(mod,n2lo,n2hi,s2); kmulum(mod,s1,s2,mid2); |
subum(mod,n1hi,n1lo,s1); subum(mod,n2lo,n2hi,s2); kmulum(mod,s1,s2,mid2); |
addum(mod,mid1,mid2,mid); |
addum(mod,mid1,mid2,mid); |
if ( mid ) { |
if ( mid ) { |
len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1; |
len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1; |
bzero((char *)COEF(t2),len*sizeof(int)); |
bzero((char *)COEF(t2),len*sizeof(int)); |
bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int)); |
bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int)); |
addum(mod,t1,t2,nr); |
addum(mod,t1,t2,nr); |
} else |
} else |
copyum(t1,nr); |
copyum(t1,nr); |
} |
} |
|
|
void ksquareummain(mod,n1,nr) |
void ksquareummain(int mod,UM n1,UM nr) |
int mod; |
|
UM n1,nr; |
|
{ |
{ |
int d1,h,len; |
int d1,h,len; |
UM n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
UM n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2; |
|
|
d1 = DEG(n1)+1; h = (d1+1)/2; |
d1 = DEG(n1)+1; h = (d1+1)/2; |
n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1); |
n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1); |
lo = W_UMALLOC(2*d1+1); hi = W_UMALLOC(2*d1+1); |
lo = W_UMALLOC(2*d1+1); hi = W_UMALLOC(2*d1+1); |
mid1 = W_UMALLOC(2*d1+1); mid2 = W_UMALLOC(2*d1+1); |
mid1 = W_UMALLOC(2*d1+1); mid2 = W_UMALLOC(2*d1+1); |
mid = W_UMALLOC(2*d1+1); |
mid = W_UMALLOC(2*d1+1); |
s1 = W_UMALLOC(2*d1+1); |
s1 = W_UMALLOC(2*d1+1); |
extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi); |
extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi); |
ksquareum(mod,n1hi,hi); ksquareum(mod,n1lo,lo); |
ksquareum(mod,n1hi,hi); ksquareum(mod,n1lo,lo); |
len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1; |
len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1; |
bzero((char *)COEF(t1),len*sizeof(int)); |
bzero((char *)COEF(t1),len*sizeof(int)); |
if ( lo ) |
if ( lo ) |
bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int)); |
bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int)); |
if ( hi ) |
if ( hi ) |
bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int)); |
bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int)); |
|
|
addum(mod,hi,lo,mid1); |
addum(mod,hi,lo,mid1); |
subum(mod,n1hi,n1lo,s1); ksquareum(mod,s1,mid2); |
subum(mod,n1hi,n1lo,s1); ksquareum(mod,s1,mid2); |
subum(mod,mid1,mid2,mid); |
subum(mod,mid1,mid2,mid); |
if ( mid ) { |
if ( mid ) { |
len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1; |
len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1; |
bzero((char *)COEF(t2),len*sizeof(int)); |
bzero((char *)COEF(t2),len*sizeof(int)); |
bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int)); |
bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int)); |
addum(mod,t1,t2,nr); |
addum(mod,t1,t2,nr); |
} else |
} else |
copyum(t1,nr); |
copyum(t1,nr); |
} |
} |