Annotation of OpenXM_contrib2/asir2018/lib/dmul, Revision 1.1
1.1 ! 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@sec.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$
! 49: */
! 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 = [];
! 66: Mod = getopt(mod);
! 67: if ( type(Mod) == -1 )
! 68: Mod = 0;
! 69: NP = length(Procs)+1;
! 70: V =var(F1);
! 71: if ( !V ) {
! 72: T = F1*F2;
! 73: if ( Mod )
! 74: return T % Mod;
! 75: else
! 76: return T;
! 77: }
! 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;
! 83: if ( Bound < 32 )
! 84: Bound = 32;
! 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();
! 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));
! 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 = [];
! 133: Mod = getopt(mod);
! 134: if ( type(Mod) == -1 )
! 135: Mod = 0;
! 136: NP = length(Procs)+1;
! 137: V =var(F1);
! 138: if ( !V ) {
! 139: T = F1^2;
! 140: if ( Mod )
! 141: return T % Mod;
! 142: else
! 143: return T;
! 144: }
! 145: D1 = deg(F1,V);
! 146: Dfft = p_mag(2*D1+1)+1;
! 147: Bound = 2*maxblen(F1)+p_mag(D1)+1;
! 148: if ( Bound < 32 )
! 149: Bound = 32;
! 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();
! 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));
! 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 = [];
! 198: Mod = getopt(mod);
! 199: if ( type(Mod) == -1 )
! 200: Mod = 0;
! 201: NP = length(Procs)+1;
! 202: V =var(F1);
! 203: if ( !V ) {
! 204: T = utrunc(F1*F2,D);
! 205: if ( Mod )
! 206: return T % Mod;
! 207: else
! 208: return T;
! 209: }
! 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;
! 215: if ( Bound < 32 )
! 216: Bound = 32;
! 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();
! 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);
! 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);
! 293: return Mhat*((MhatInv*C)%M1);
! 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);
! 301: return Mhat*((MhatInv*C)%M1);
! 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);
! 309: return Mhat*((MhatInv*C)%M1);
! 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 / asir2018 / lib