Annotation of OpenXM_contrib2/asir2000/lib/dmul, Revision 1.3
1.3 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/dmul,v 1.2 2000/01/11 06:43:37 noro Exp $
! 49: */
1.1 noro 50: #define MAX(a,b) ((a)>(b)?(a):(b))
51: #define MIN(a,b) ((a)>(b)?(b):(a))
52:
53: /* CAUTION: functions in this file are experimental. */
54:
55: /*
56: return: F1*F2
57: if option 'proc' is supplied as a list of server id's,
58: F1*F2 is calculated by distributed computation.
59: */
60:
61: def d_mul(F1,F2)
62: {
63: Procs = getopt(proc);
64: if ( type(Procs) == -1 )
65: Procs = [];
1.2 noro 66: Mod = getopt(mod);
67: if ( type(Mod) == -1 )
68: Mod = 0;
1.1 noro 69: NP = length(Procs)+1;
70: V =var(F1);
1.2 noro 71: if ( !V ) {
72: T = F1*F2;
73: if ( Mod )
74: return T % Mod;
75: else
76: return T;
77: }
1.1 noro 78: D1 = deg(F1,V);
79: D2 = deg(F2,V);
80: Dmin = MIN(D1,D2);
81: Dfft = p_mag(D1+D2+1)+1;
82: Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
1.2 noro 83: if ( Bound < 32 )
84: Bound = 32;
1.1 noro 85: Marray = newvect(NP);
86: MIarray = newvect(NP);
87: for ( I = 0; I < NP; I++ ) {
88: Marray[I] = 1;
89: MIarray[I] = [];
90: }
91:
92: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
93: T = get_next_fft_prime(I,Dfft);
94: if ( !T )
95: error("fft_mul_d : fft_prime exhausted.");
96: Marray[J] *= T[1];
97: MIarray[J] = cons(T[0],MIarray[J]);
98: M *= T[1];
99: I = T[0]+1;
100: }
101: /* Now,
102: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
103: M = Marray[0]*...*Marray[NP-1]
104: */
105: C = newvect(NP);
106: T0 = time();
107: for ( J = 0; J < NP-1; J++ )
108: ox_cmo_rpc(Procs[J],"call_umul",F1,F2,MIarray[J],Marray[J],M);
109: T1 = time();
110: R = call_umul(F1,F2,MIarray[NP-1],Marray[NP-1],M);
111: T2 = time();
112: for ( J = 0; J < NP-1; J++ )
113: R += ox_pop_cmo(Procs[J]);
114: T3 = time();
1.2 noro 115: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
116: if ( Mod )
117: return (R%M)%Mod;
118: else
119: return uadj_coef(R%M,M,ishift(M,1));
1.1 noro 120: }
121:
122: /*
123: return: F1^2
124: if option 'proc' is supplied as a list of server id's,
125: F1^2 is calculated by distributed computation.
126: */
127:
128: def d_square(F1)
129: {
130: Procs = getopt(proc);
131: if ( type(Procs) == -1 )
132: Procs = [];
1.2 noro 133: Mod = getopt(mod);
134: if ( type(Mod) == -1 )
135: Mod = 0;
1.1 noro 136: NP = length(Procs)+1;
137: V =var(F1);
1.2 noro 138: if ( !V ) {
139: T = F1^2;
140: if ( Mod )
141: return T % Mod;
142: else
143: return T;
144: }
1.1 noro 145: D1 = deg(F1,V);
146: Dfft = p_mag(2*D1+1)+1;
147: Bound = 2*maxblen(F1)+p_mag(D1)+1;
1.2 noro 148: if ( Bound < 32 )
149: Bound = 32;
1.1 noro 150: Marray = newvect(NP);
151: MIarray = newvect(NP);
152: for ( I = 0; I < NP; I++ ) {
153: Marray[I] = 1;
154: MIarray[I] = [];
155: }
156:
157: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
158: T = get_next_fft_prime(I,Dfft);
159: if ( !T )
160: error("fft_mul_d : fft_prime exhausted.");
161: Marray[J] *= T[1];
162: MIarray[J] = cons(T[0],MIarray[J]);
163: M *= T[1];
164: I = T[0]+1;
165: }
166: /* Now,
167: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
168: M = Marray[0]*...*Marray[NP-1]
169: */
170: C = newvect(NP);
171: T0 = time();
172: for ( J = 0; J < NP-1; J++ )
173: ox_cmo_rpc(Procs[J],"call_usquare",F1,MIarray[J],Marray[J],M);
174: T1 = time();
175: R = call_usquare(F1,MIarray[NP-1],Marray[NP-1],M);
176: T2 = time();
177: for ( J = 0; J < NP-1; J++ )
178: R += ox_pop_cmo(Procs[J]);
179: T3 = time();
1.2 noro 180: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
181: if ( Mod )
182: return (R%M)%Mod;
183: else
184: return uadj_coef(R%M,M,ishift(M,1));
1.1 noro 185: }
186:
187: /*
188: return: F1^2 mod V^(D+1)
189: if option 'proc' is supplied as a list of server id's,
190: F1*F2 mod V^(D+1) is calculated by distributed computation.
191: */
192:
193: def d_tmul(F1,F2,D)
194: {
195: Procs = getopt(proc);
196: if ( type(Procs) == -1 )
197: Procs = [];
1.2 noro 198: Mod = getopt(mod);
199: if ( type(Mod) == -1 )
200: Mod = 0;
1.1 noro 201: NP = length(Procs)+1;
202: V =var(F1);
1.2 noro 203: if ( !V ) {
204: T = utrunc(F1*F2,D);
205: if ( Mod )
206: return T % Mod;
207: else
208: return T;
209: }
1.1 noro 210: D1 = deg(F1,V);
211: D2 = deg(F2,V);
212: Dmin = MIN(D1,D2);
213: Dfft = p_mag(D1+D2+1)+1;
214: Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
1.2 noro 215: if ( Bound < 32 )
216: Bound = 32;
1.1 noro 217: Marray = newvect(NP);
218: MIarray = newvect(NP);
219: for ( I = 0; I < NP; I++ ) {
220: Marray[I] = 1;
221: MIarray[I] = [];
222: }
223:
224: for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
225: T = get_next_fft_prime(I,Dfft);
226: if ( !T )
227: error("fft_mul_d : fft_prime exhausted.");
228: Marray[J] *= T[1];
229: MIarray[J] = cons(T[0],MIarray[J]);
230: M *= T[1];
231: I = T[0]+1;
232: }
233: /* Now,
234: Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
235: M = Marray[0]*...*Marray[NP-1]
236: */
237: C = newvect(NP);
238: T0 = time();
239: for ( J = 0; J < NP-1; J++ )
240: ox_cmo_rpc(Procs[J],"call_utmul",F1,F2,D,MIarray[J],Marray[J],M);
241: T1 = time();
242: R = call_utmul(F1,F2,D,MIarray[NP-1],Marray[NP-1],M);
243: T2 = time();
244: for ( J = 0; J < NP-1; J++ )
245: R += ox_pop_cmo(Procs[J]);
246: T3 = time();
1.2 noro 247: /* print(["send",T1[3]-T0[3],"self",T2[3]-T1[3],"recv",T3[3]-T2[3]]); */
248: if ( Mod )
249: return (R%M)%Mod;
250: else
251: return uadj_coef(R%M,M,ishift(M,1));
252: }
253:
254: def d_rembymul(F1,F2,INVF2)
255: {
256: Procs = getopt(proc);
257: if ( type(Procs) == -1 )
258: Procs = [];
259: Mod = getopt(mod);
260: if ( type(Mod) == -1 )
261: Mod = 0;
262: NP = length(Procs)+1;
263: if ( !F2 )
264: error("d_rembymul : division by 0");
265: V =var(F1);
266: if ( !V ) {
267: T = srem(F1,F2);
268: if ( Mod )
269: return T % Mod;
270: else
271: return T;
272: }
273: D1 = deg(F1,V);
274: D2 = deg(F2,V);
275: if ( !F1 || !D2 )
276: return 0;
277: if ( D1 < D2 )
278: return F1;
279: D = D1-D2;
280: R1 = utrunc(ureverse(F1),D);
281: Q = ureverse(utrunc(d_tmul(R1,INVF2,D|proc=Procs,mod=Mod),D));
282: if ( Mod )
283: return (utrunc(F1,D2-1)-d_tmul(Q,F2,D2-1|proc=Procs,mod=Mod))%Mod;
284: else
285: return utrunc(F1,D2-1)-d_tmul(Q,F2,D2-1|proc=Procs);
1.1 noro 286: }
287:
288: def call_umul(F1,F2,Ind,M1,M)
289: {
290: C = umul_specialmod(F1,F2,Ind);
291: Mhat = idiv(M,M1);
292: MhatInv = inv(Mhat,M1);
1.2 noro 293: return Mhat*((MhatInv*C)%M1);
1.1 noro 294: }
295:
296: def call_usquare(F1,Ind,M1,M)
297: {
298: C = usquare_specialmod(F1,Ind);
299: Mhat = idiv(M,M1);
300: MhatInv = inv(Mhat,M1);
1.2 noro 301: return Mhat*((MhatInv*C)%M1);
1.1 noro 302: }
303:
304: def call_utmul(F1,F2,D,Ind,M1,M)
305: {
306: C = utmul_specialmod(F1,F2,D,Ind);
307: Mhat = idiv(M,M1);
308: MhatInv = inv(Mhat,M1);
1.2 noro 309: return Mhat*((MhatInv*C)%M1);
1.1 noro 310: }
311: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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