![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dmul CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to dmul CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/lib/dmul,v 1.1 1999/12/27 04:16:32 noro Exp $ */
1.1 noro 2: #define MAX(a,b) ((a)>(b)?(a):(b))
3: #define MIN(a,b) ((a)>(b)?(b):(a))
4:
5: /* CAUTION: functions in this file are experimental. */
6:
7: /*
8: return: F1*F2
9: if option 'proc' is supplied as a list of server id's,
10: F1*F2 is calculated by distributed computation.
11: */
12:
13: def d_mul(F1,F2)
14: {
15: Procs = getopt(proc);
16: if ( type(Procs) == -1 )
17: Procs = [];
1.2 ! noro 18: Mod = getopt(mod);
! 19: if ( type(Mod) == -1 )
! 20: Mod = 0;
1.1 noro 21: NP = length(Procs)+1;
22: V =var(F1);
1.2 ! noro 23: if ( !V ) {
! 24: T = F1*F2;
! 25: if ( Mod )
! 26: return T % Mod;
! 27: else
! 28: return T;
! 29: }
1.1 noro 30: D1 = deg(F1,V);
31: D2 = deg(F2,V);
32: Dmin = MIN(D1,D2);
33: Dfft = p_mag(D1+D2+1)+1;
34: Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
1.2 ! noro 35: if ( Bound < 32 )
! 36: Bound = 32;
1.1 noro 37: Marray = newvect(NP);
38: MIarray = newvect(NP);
39: for ( I = 0; I < NP; I++ ) {
40: Marray[I] = 1;
41: MIarray[I] = [];
42: }
43:
44: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
45: T = get_next_fft_prime(I,Dfft);
46: if ( !T )
47: error("fft_mul_d : fft_prime exhausted.");
48: Marray[J] *= T[1];
49: MIarray[J] = cons(T[0],MIarray[J]);
50: M *= T[1];
51: I = T[0]+1;
52: }
53: /* Now,
54: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
55: M = Marray[0]*...*Marray[NP-1]
56: */
57: C = newvect(NP);
58: T0 = time();
59: for ( J = 0; J < NP-1; J++ )
60: ox_cmo_rpc(Procs[J],"call_umul",F1,F2,MIarray[J],Marray[J],M);
61: T1 = time();
62: R = call_umul(F1,F2,MIarray[NP-1],Marray[NP-1],M);
63: T2 = time();
64: for ( J = 0; J < NP-1; J++ )
65: R += ox_pop_cmo(Procs[J]);
66: T3 = time();
1.2 ! noro 67: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
! 68: if ( Mod )
! 69: return (R%M)%Mod;
! 70: else
! 71: return uadj_coef(R%M,M,ishift(M,1));
1.1 noro 72: }
73:
74: /*
75: return: F1^2
76: if option 'proc' is supplied as a list of server id's,
77: F1^2 is calculated by distributed computation.
78: */
79:
80: def d_square(F1)
81: {
82: Procs = getopt(proc);
83: if ( type(Procs) == -1 )
84: Procs = [];
1.2 ! noro 85: Mod = getopt(mod);
! 86: if ( type(Mod) == -1 )
! 87: Mod = 0;
1.1 noro 88: NP = length(Procs)+1;
89: V =var(F1);
1.2 ! noro 90: if ( !V ) {
! 91: T = F1^2;
! 92: if ( Mod )
! 93: return T % Mod;
! 94: else
! 95: return T;
! 96: }
1.1 noro 97: D1 = deg(F1,V);
98: Dfft = p_mag(2*D1+1)+1;
99: Bound = 2*maxblen(F1)+p_mag(D1)+1;
1.2 ! noro 100: if ( Bound < 32 )
! 101: Bound = 32;
1.1 noro 102: Marray = newvect(NP);
103: MIarray = newvect(NP);
104: for ( I = 0; I < NP; I++ ) {
105: Marray[I] = 1;
106: MIarray[I] = [];
107: }
108:
109: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
110: T = get_next_fft_prime(I,Dfft);
111: if ( !T )
112: error("fft_mul_d : fft_prime exhausted.");
113: Marray[J] *= T[1];
114: MIarray[J] = cons(T[0],MIarray[J]);
115: M *= T[1];
116: I = T[0]+1;
117: }
118: /* Now,
119: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
120: M = Marray[0]*...*Marray[NP-1]
121: */
122: C = newvect(NP);
123: T0 = time();
124: for ( J = 0; J < NP-1; J++ )
125: ox_cmo_rpc(Procs[J],"call_usquare",F1,MIarray[J],Marray[J],M);
126: T1 = time();
127: R = call_usquare(F1,MIarray[NP-1],Marray[NP-1],M);
128: T2 = time();
129: for ( J = 0; J < NP-1; J++ )
130: R += ox_pop_cmo(Procs[J]);
131: T3 = time();
1.2 ! noro 132: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
! 133: if ( Mod )
! 134: return (R%M)%Mod;
! 135: else
! 136: return uadj_coef(R%M,M,ishift(M,1));
1.1 noro 137: }
138:
139: /*
140: return: F1^2 mod V^(D+1)
141: if option 'proc' is supplied as a list of server id's,
142: F1*F2 mod V^(D+1) is calculated by distributed computation.
143: */
144:
145: def d_tmul(F1,F2,D)
146: {
147: Procs = getopt(proc);
148: if ( type(Procs) == -1 )
149: Procs = [];
1.2 ! noro 150: Mod = getopt(mod);
! 151: if ( type(Mod) == -1 )
! 152: Mod = 0;
1.1 noro 153: NP = length(Procs)+1;
154: V =var(F1);
1.2 ! noro 155: if ( !V ) {
! 156: T = utrunc(F1*F2,D);
! 157: if ( Mod )
! 158: return T % Mod;
! 159: else
! 160: return T;
! 161: }
1.1 noro 162: D1 = deg(F1,V);
163: D2 = deg(F2,V);
164: Dmin = MIN(D1,D2);
165: Dfft = p_mag(D1+D2+1)+1;
166: Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
1.2 ! noro 167: if ( Bound < 32 )
! 168: Bound = 32;
1.1 noro 169: Marray = newvect(NP);
170: MIarray = newvect(NP);
171: for ( I = 0; I < NP; I++ ) {
172: Marray[I] = 1;
173: MIarray[I] = [];
174: }
175:
176: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
177: T = get_next_fft_prime(I,Dfft);
178: if ( !T )
179: error("fft_mul_d : fft_prime exhausted.");
180: Marray[J] *= T[1];
181: MIarray[J] = cons(T[0],MIarray[J]);
182: M *= T[1];
183: I = T[0]+1;
184: }
185: /* Now,
186: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
187: M = Marray[0]*...*Marray[NP-1]
188: */
189: C = newvect(NP);
190: T0 = time();
191: for ( J = 0; J < NP-1; J++ )
192: ox_cmo_rpc(Procs[J],"call_utmul",F1,F2,D,MIarray[J],Marray[J],M);
193: T1 = time();
194: R = call_utmul(F1,F2,D,MIarray[NP-1],Marray[NP-1],M);
195: T2 = time();
196: for ( J = 0; J < NP-1; J++ )
197: R += ox_pop_cmo(Procs[J]);
198: T3 = time();
1.2 ! noro 199: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
! 200: if ( Mod )
! 201: return (R%M)%Mod;
! 202: else
! 203: return uadj_coef(R%M,M,ishift(M,1));
! 204: }
! 205:
! 206: def d_rembymul(F1,F2,INVF2)
! 207: {
! 208: Procs = getopt(proc);
! 209: if ( type(Procs) == -1 )
! 210: Procs = [];
! 211: Mod = getopt(mod);
! 212: if ( type(Mod) == -1 )
! 213: Mod = 0;
! 214: NP = length(Procs)+1;
! 215: if ( !F2 )
! 216: error("d_rembymul : division by 0");
! 217: V =var(F1);
! 218: if ( !V ) {
! 219: T = srem(F1,F2);
! 220: if ( Mod )
! 221: return T % Mod;
! 222: else
! 223: return T;
! 224: }
! 225: D1 = deg(F1,V);
! 226: D2 = deg(F2,V);
! 227: if ( !F1 || !D2 )
! 228: return 0;
! 229: if ( D1 < D2 )
! 230: return F1;
! 231: D = D1-D2;
! 232: R1 = utrunc(ureverse(F1),D);
! 233: Q = ureverse(utrunc(d_tmul(R1,INVF2,D|proc=Procs,mod=Mod),D));
! 234: if ( Mod )
! 235: return (utrunc(F1,D2-1)-d_tmul(Q,F2,D2-1|proc=Procs,mod=Mod))%Mod;
! 236: else
! 237: return utrunc(F1,D2-1)-d_tmul(Q,F2,D2-1|proc=Procs);
1.1 noro 238: }
239:
240: def call_umul(F1,F2,Ind,M1,M)
241: {
242: C = umul_specialmod(F1,F2,Ind);
243: Mhat = idiv(M,M1);
244: MhatInv = inv(Mhat,M1);
1.2 ! noro 245: return Mhat*((MhatInv*C)%M1);
1.1 noro 246: }
247:
248: def call_usquare(F1,Ind,M1,M)
249: {
250: C = usquare_specialmod(F1,Ind);
251: Mhat = idiv(M,M1);
252: MhatInv = inv(Mhat,M1);
1.2 ! noro 253: return Mhat*((MhatInv*C)%M1);
1.1 noro 254: }
255:
256: def call_utmul(F1,F2,D,Ind,M1,M)
257: {
258: C = utmul_specialmod(F1,F2,D,Ind);
259: Mhat = idiv(M,M1);
260: MhatInv = inv(Mhat,M1);
1.2 ! noro 261: return Mhat*((MhatInv*C)%M1);
1.1 noro 262: }
263: end$