===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/up.c,v
retrieving revision 1.1.1.1
retrieving revision 1.5
diff -u -p -r1.1.1.1 -r1.5
--- OpenXM_contrib2/asir2000/engine/up.c	1999/12/03 07:39:08	1.1.1.1
+++ OpenXM_contrib2/asir2000/engine/up.c	2001/10/09 01:36:13	1.5
@@ -1,4 +1,52 @@
-/* $OpenXM: OpenXM/src/asir99/engine/up.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/up.c,v 1.4 2000/08/22 05:04:06 noro Exp $ 
+*/
 #include "ca.h"
 #include <math.h>
 
@@ -28,9 +76,7 @@ extern int lm_lazy;
 extern int current_ff;
 extern int GC_dont_gc;
 
-void monicup(a,b)
-UP a;
-UP *b;
+void monicup(UP a,UP *b)
 {
 	UP w;
 
@@ -43,9 +89,7 @@ UP *b;
 	}
 }
 
-void simpup(a,b)
-UP a;
-UP *b;
+void simpup(UP a,UP *b)
 {
 	int i,d;
 	UP c;
@@ -68,9 +112,7 @@ UP *b;
 	}
 }
 
-void simpnum(a,b)
-Num a;
-Num *b;
+void simpnum(Num a,Num *b)
 {
 	LM lm;
 	GF2N gf;
@@ -91,9 +133,7 @@ Num *b;
 		}
 }
 
-void uremp(p1,p2,rp)
-P p1,p2;
-P *rp;
+void uremp(P p1,P p2,P *rp)
 {
 	UP n1,n2,r;
 
@@ -106,9 +146,7 @@ P *rp;
 	}
 }
 
-void ugcdp(p1,p2,rp)
-P p1,p2;
-P *rp;
+void ugcdp(P p1,P p2,P *rp)
 {
 	UP n1,n2,r;
 
@@ -117,10 +155,7 @@ P *rp;
 	uptop(r,rp);
 }
 
-void reversep(p1,d,rp)
-P p1;
-Q d;
-P *rp;
+void reversep(P p1,Q d,P *rp)
 {
 	UP n1,r;
 
@@ -133,10 +168,7 @@ P *rp;
 	}
 }
 
-void invmodp(p1,d,rp)
-P p1;
-Q d;
-P *rp;
+void invmodp(P p1,Q d,P *rp)
 {
 	UP n1,r;
 
@@ -149,9 +181,7 @@ P *rp;
 	}
 }
 
-void addup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void addup(UP n1,UP n2,UP *nr)
 {
 	UP r,t;
 	int i,d1,d2;
@@ -180,9 +210,7 @@ UP *nr;
 	}
 }
 
-void subup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void subup(UP n1,UP n2,UP *nr)
 {
 	UP r;
 	int i,d1,d2,d;
@@ -215,9 +243,7 @@ UP *nr;
 	}
 }
 
-void chsgnup(n1,nr)
-UP n1;
-UP *nr;
+void chsgnup(UP n1,UP *nr)
 {
 	UP r;
 	int d1,i;
@@ -234,10 +260,7 @@ UP *nr;
 	}
 }
 
-void hybrid_mulup(ff,n1,n2,nr)
-int ff;
-UP n1,n2;
-UP *nr;
+void hybrid_mulup(int ff,UP n1,UP n2,UP *nr)
 {
 	if ( !n1 || !n2 )
 		*nr = 0;
@@ -254,10 +277,7 @@ UP *nr;
 		}
 }
 
-void hybrid_squareup(ff,n1,nr)
-int ff;
-UP n1;
-UP *nr;
+void hybrid_squareup(int ff,UP n1,UP *nr)
 {
 	if ( !n1 )
 		*nr = 0;
@@ -274,11 +294,7 @@ UP *nr;
 		}
 }
 
-void hybrid_tmulup(ff,n1,n2,d,nr)
-int ff;
-UP n1,n2;
-int d;
-UP *nr;
+void hybrid_tmulup(int ff,UP n1,UP n2,int d,UP *nr)
 {
 	if ( !n1 || !n2 )
 		*nr = 0;
@@ -295,9 +311,7 @@ UP *nr;
 		}
 }
 
-void mulup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void mulup(UP n1,UP n2,UP *nr)
 {
 	UP r;
 	Num *pc1,*pc,*c1,*c2,*c;
@@ -326,10 +340,7 @@ UP *nr;
 
 /* nr = c*n1 */
 
-void mulcup(c,n1,nr)
-Num c;
-UP n1;
-UP *nr;
+void mulcup(Num c,UP n1,UP *nr)
 {
 	int d;
 	UP r;
@@ -354,10 +365,7 @@ UP *nr;
 	}
 }
 
-void tmulup(n1,n2,d,nr)
-UP n1,n2;
-int d;
-UP *nr;
+void tmulup(UP n1,UP n2,int d,UP *nr)
 {
 	UP r;
 	Num *pc1,*pc,*c1,*c2,*c;
@@ -391,9 +399,7 @@ UP *nr;
 	}
 }
 
-void squareup(n1,nr)
-UP n1;
-UP *nr;
+void squareup(UP n1,UP *nr)
 {
 	UP r;
 	Num *c1,*c;
@@ -428,9 +434,7 @@ UP *nr;
 	}
 }
 
-void remup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void remup(UP n1,UP n2,UP *nr)
 {
 	UP w,r;
 
@@ -451,8 +455,7 @@ UP *nr;
 	}
 }
 
-void remup_destructive(n1,n2)
-UP n1,n2;
+void remup_destructive(UP n1,UP n2)
 {
 	Num *c1,*c2;
 	int d1,d2,i,j;
@@ -480,9 +483,7 @@ UP n1,n2;
 	n1->d = i;
 }
 
-void qrup(n1,n2,nq,nr)
-UP n1,n2;
-UP *nq,*nr;
+void qrup(UP n1,UP n2,UP *nq,UP *nr)
 {
 	UP w,r,q;
 	struct oUP t;
@@ -517,8 +518,7 @@ UP *nq,*nr;
 	}
 }
 
-void qrup_destructive(n1,n2)
-UP n1,n2;
+void qrup_destructive(UP n1,UP n2)
 {
 	Num *c1,*c2;
 	int d1,d2,i,j;
@@ -548,9 +548,7 @@ UP n1,n2;
 	n1->d = i;
 }
 
-void gcdup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void gcdup(UP n1,UP n2,UP *nr)
 {
 	UP w1,w2,t,r;
 	int d1,d2;
@@ -588,9 +586,7 @@ UP *nr;
 
 /* compute r s.t. a*r = 1 mod m */
 
-void extended_gcdup(a,m,r)
-UP a,m;
-UP *r;
+void extended_gcdup(UP a,UP m,UP *r)
 {
 	UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;
 	Num i;
@@ -611,10 +607,7 @@ UP *r;
 	mulup(b2,inv,r);
 }
 
-void reverseup(n1,d,nr)
-UP n1;
-int d;
-UP *nr;
+void reverseup(UP n1,int d,UP *nr)
 {
 	Num *c1,*c;
 	int i,d1;
@@ -638,10 +631,7 @@ UP *nr;
 	}
 }
 
-void invmodup(n1,d,nr)
-UP n1;
-int d;
-UP *nr;
+void invmodup(UP n1,int d,UP *nr)
 {
 	UP r;
 	Num s,t,u,hinv;
@@ -679,10 +669,7 @@ UP *nr;
 	}
 }
 
-void pwrup(n,e,nr)
-UP n;
-Q e;
-UP *nr;
+void pwrup(UP n,Q e,UP *nr)
 {
 	UP y,z,t;
 	N en,qn;
@@ -710,8 +697,7 @@ UP *nr;
 	}
 }
 
-int compup(n1,n2)
-UP n1,n2;
+int compup(UP n1,UP n2)
 {
 	int i,r;
 
@@ -733,10 +719,7 @@ UP n1,n2;
 	}
 }
 
-void kmulp(vl,n1,n2,nr)
-VL vl;
-P n1,n2;
-P *nr;
+void kmulp(VL vl,P n1,P n2,P *nr)
 {
 	UP b1,b2,br;
 
@@ -751,10 +734,7 @@ P *nr;
 	}
 }
 
-void ksquarep(vl,n1,nr)
-VL vl;
-P n1;
-P *nr;
+void ksquarep(VL vl,P n1,P *nr)
 {
 	UP b1,br;
 
@@ -769,12 +749,9 @@ P *nr;
 	}
 }
 
-void kmulup(n1,n2,nr)
-UP n1,n2,*nr;
+void kmulup(UP n1,UP n2,UP *nr)
 {
-	UP n,t,s,m,carry;
-	int d,d1,d2,len,i,l;
-	pointer *r,*r0;
+	int d1,d2;
 
 	if ( !n1 || !n2 ) {
 		*nr = 0; return;
@@ -786,8 +763,7 @@ UP n1,n2,*nr;
 		kmulupmain(n1,n2,nr);
 }
 
-void ksquareup(n1,nr)
-UP n1,*nr;
+void ksquareup(UP n1,UP *nr)
 {
 	int d1;
 	extern Q TWO;
@@ -802,26 +778,21 @@ UP n1,*nr;
 		ksquareupmain(n1,nr);
 }
 
-void copyup(n1,n2)
-UP n1,n2;
+void copyup(UP n1,UP n2)
 {
 	n2->d = n1->d;
 	bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));
 }
 
-void c_copyup(n,len,p)
-UP n;
-int len;
-pointer *p;
+void c_copyup(UP n,int len,pointer *p)
 {
 	if ( n )
 		bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));
 }
 
-void kmulupmain(n1,n2,nr)
-UP n1,n2,*nr;
+void kmulupmain(UP n1,UP n2,UP *nr)
 {
-	int d1,d2,h,len;
+	int d1,d2,h;
 	UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
 
 	d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;
@@ -838,10 +809,9 @@ UP n1,n2,*nr;
 	addup(t1,t2,nr);
 }
 
-void ksquareupmain(n1,nr)
-UP n1,*nr;
+void ksquareupmain(UP n1,UP *nr)
 {
-	int d1,h,len;
+	int d1,h;
 	UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
 
 	d1 = DEG(n1)+1; h = (d1+1)/2;
@@ -856,9 +826,7 @@ UP n1,*nr;
 	addup(t1,t2,nr);
 }
 
-void rembymulup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void rembymulup(UP n1,UP n2,UP *nr)
 {
 	int d1,d2,d;
 	UP r1,r2,inv2,t,s,q;
@@ -887,13 +855,10 @@ UP *nr;
 	inv2 = inverse of reversep(n2) mod x^d
 */
 
-void hybrid_rembymulup_special(ff,n1,n2,inv2,nr)
-int ff;
-UP n1,n2,inv2;
-UP *nr;
+void hybrid_rembymulup_special(int ff,UP n1,UP n2,UP inv2,UP *nr)
 {
 	int d1,d2,d;
-	UP r1,t,s,q,u;
+	UP r1,t,s,q;
 
 	if ( !n2 )
 		error("hybrid_rembymulup : division by 0");
@@ -914,12 +879,10 @@ UP *nr;
 	}
 }
 
-void rembymulup_special(n1,n2,inv2,nr)
-UP n1,n2,inv2;
-UP *nr;
+void rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
 {
 	int d1,d2,d;
-	UP r1,t,s,q,u;
+	UP r1,t,s,q;
 
 	if ( !n2 )
 		error("rembymulup : division by 0");
@@ -942,13 +905,10 @@ UP *nr;
 
 /* *nr = n1*n2 mod x^d */
 
-void tkmulup(n1,n2,d,nr)
-UP n1,n2;
-int d;
-UP *nr;
+void tkmulup(UP n1,UP n2,int d,UP *nr)
 {
 	int m;
-	UP n1l,n1h,n2l,n2h,l,h,t,s,u,afo;
+	UP n1l,n1h,n2l,n2h,l,h,t,s,u;
 
 	if ( d < 0 )
 		error("tkmulup : invalid argument");
@@ -978,10 +938,7 @@ UP *nr;
 
 /* n->n*x^d */
 
-void shiftup(n,d,nr)
-UP n;
-int d;
-UP *nr;
+void shiftup(UP n,int d,UP *nr)
 {
 	int dr;
 	UP r;
@@ -996,9 +953,7 @@ UP *nr;
 	}
 }
 
-void fft_rembymulup_special(n1,n2,inv2,nr)
-UP n1,n2,inv2;
-UP *nr;
+void fft_rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
 {
 	int d1,d2,d;
 	UP r1,t,s,q,u;
@@ -1021,9 +976,7 @@ UP *nr;
 	}
 }
 
-void set_degreeup(n,d)
-UP n;
-int d;
+void set_degreeup(UP n,int d)
 {
 	int i;
 
@@ -1033,10 +986,7 @@ int d;
 
 /* n -> n0 + x^d n1 */
 
-void decompup(n,d,n0,n1)
-UP n;
-int d;
-UP *n0,*n1;
+void decompup(UP n,int d,UP *n0,UP *n1)
 {
 	int dn;
 	UP r0,r1;
@@ -1064,10 +1014,7 @@ UP *n0,*n1;
 
 /* n -> n mod x^d */
 
-void truncup(n1,d,nr)
-UP n1;
-int d;
-UP *nr;
+void truncup(UP n1,int d,UP *nr)
 {
 	int i;
 	UP r;
@@ -1089,8 +1036,7 @@ UP *nr;
 	}
 }
 
-int int_bits(t)
-int t;
+int int_bits(int t)
 {
 	int k;
 
@@ -1100,8 +1046,7 @@ int t;
 
 /* n is assumed to be LM or integer coefficient */
 
-int maxblenup(n)
-UP n;
+int maxblenup(UP n)
 {
 	int m,r,i,d;
 	Num *c;
@@ -1126,10 +1071,7 @@ UP n;
 	return r;
 }
 
-void uptofmarray(mod,n,f)
-int mod;
-UP n;
-ModNum *f;
+void uptofmarray(int mod,UP n,ModNum *f)
 {
 	int d,i;
 	unsigned int r;
@@ -1160,10 +1102,7 @@ ModNum *f;
 	}
 }
 
-void fmarraytoup(f,d,nr)
-ModNum *f;
-int d;
-UP *nr;
+void fmarraytoup(ModNum *f,int d,UP *nr)
 {
 	int i;
 	Q *c;
@@ -1186,10 +1125,7 @@ UP *nr;
 
 /* f[i]: an array of length n */
 
-void uiarraytoup(f,n,d,nr)
-unsigned int **f;
-int n,d;
-UP *nr;
+void uiarraytoup(unsigned int **f,int n,int d,UP *nr)
 {
 	int i,j;
 	unsigned int *fi;
@@ -1221,10 +1157,7 @@ UP *nr;
 	}
 }
 
-void adj_coefup(n,m,m2,nr)
-UP n;
-N m,m2;
-UP *nr;
+void adj_coefup(UP n,N m,N m2,UP *nr)
 {
 	int d;
 	Q *c,*cr;
@@ -1256,10 +1189,7 @@ UP *nr;
 
 /* n is assumed to have positive integer coefficients. */
 
-void remcup(n,mod,nr)
-UP n;
-N mod;
-UP *nr;
+void remcup(UP n,N mod,UP *nr)
 {
 	int i,d;
 	Q *c,*cr;
@@ -1290,9 +1220,7 @@ UP *nr;
 
 void fft_mulup_main(UP,UP,int,UP *);
 
-void fft_mulup(n1,n2,nr)
-UP n1,n2;
-UP *nr;
+void fft_mulup(UP n1,UP n2,UP *nr)
 {
 	int d1,d2,d,b1,b2,h;
 	UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;
@@ -1326,10 +1254,7 @@ UP *nr;
 	}
 }
 
-void trunc_fft_mulup(n1,n2,dbd,nr)
-UP n1,n2;
-int dbd;
-UP *nr;
+void trunc_fft_mulup(UP n1,UP n2,int dbd,UP *nr)
 {
 	int d1,d2,b1,b2,m;
 	UP n1l,n1h,n2l,n2h,l,h,t,s,u;
@@ -1366,11 +1291,9 @@ UP *nr;
 	}
 }
 
-void fft_squareup(n1,nr)
-UP n1;
-UP *nr;
+void fft_squareup(UP n1,UP *nr)
 {
-	int d1,d,h,len,b1;
+	int d1,d,h,b1;
 	UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
 
 	if ( !n1 )
@@ -1407,9 +1330,7 @@ UP *nr;
  * n1 == n2 => squaring
  */
 
-void fft_mulup_main(n1,n2,dbd,nr)
-UP n1,n2;
-UP *nr;
+void fft_mulup_main(UP n1,UP n2,int dbd,UP *nr)
 {
 	ModNum *f1,*f2,*w,*fr;
 	ModNum **frarray,**fa;
@@ -1478,13 +1399,7 @@ UP *nr;
 #if 0
 /* inefficient version */
 
-void crup(f,d,mod,index,m,r)
-ModNum **f;
-int d;
-int *mod;
-int index;
-N m;
-UP *r;
+void crup(ModNum **f,int d,int *mod,int index,N m,UP *r)
 {
 	N *cof,*c;
 	int *inv;
@@ -1523,13 +1438,7 @@ UP *r;
 #else
 /* improved version */
 
-void crup(f,d,mod,index,m,r)
-ModNum **f;
-int d;
-int *mod;
-int index;
-N m;
-UP *r;
+void crup(ModNum **f,int d,int *mod,int index,N m,UP *r)
 {
 	N cof,c,t,w,w1;
 	struct oN fc;
@@ -1582,3 +1491,57 @@ UP *r;
 }
 #endif
 
+/* 
+ * dbd == 0 => n1 * n2
+ * dbd > 0  => n1 * n2 mod x^dbd
+ * n1 == n2 => squaring
+ * return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime
+ */
+
+void fft_mulup_specialmod_main(UP n1,UP n2,int dbd,int *modind,int nmod,UP *nr)
+{
+	ModNum *f1,*f2,*w,*fr;
+	ModNum **frarray;
+	N m,m1,m2;
+	unsigned int *modarray;
+	int d1,d2,dmin,i,root,d,cond,bound;
+
+	if ( !n1 || !n2 ) {
+		*nr = 0; return;
+	}
+	d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
+	if ( !d1 || !d2 ) {
+		mulup(n1,n2,nr); return;
+	}
+	m = ONEN;
+	bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
+	f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
+	if ( n1 == n2 )
+		f2 = 0;
+	else
+		f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
+	w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
+	frarray = (ModNum **)MALLOC(nmod*sizeof(ModNum *));
+	modarray = (unsigned int *)MALLOC_ATOMIC(nmod*sizeof(unsigned int *));
+
+	for ( i = 0; i < nmod; i++ ) {
+		FFT_primes(modind[i],&modarray[i],&root,&d);
+			if ( (1<<d) < d1+d2+1 )
+				error("fft_mulup_specialmod_main : invalid modulus");
+		frarray[i] = fr 
+			= (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
+		uptofmarray(modarray[i],n1,f1); 
+		if ( !f2 )
+			cond = FFT_pol_square(d1,f1,fr,modind[i],w);
+		else {
+			uptofmarray(modarray[i],n2,f2); 
+			cond = FFT_pol_product(d1,f1,d2,f2,fr,modind[i],w);
+		}
+		if ( cond )
+			error("fft_mulup_specialmod_main : error in FFT_pol_product");
+		STON(modarray[i],m1); muln(m,m1,&m2); m = m2;
+	}
+	if ( !dbd )
+		dbd = d1+d2+1;
+	crup(frarray,MIN(d1+d2,dbd-1),modarray,nmod,m,nr);
+}