Return to dmul102 CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/lib/dmul102, Revision 1.3

1.3     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2000/lib/dmul102,v 1.2 2003/12/12 09:17:31 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: Hosts = [
                      8: "iyokan-0-g","iyokan-0-g",
                      9: "iyokan-1-g","iyokan-1-g",
                     10: "iyokan-2-g","iyokan-2-g",
                     11: "iyokan-3-g","iyokan-3-g"
                     12: ]$
                     13:
                     14: def spawn102(Hosts,N)
1.3     ! noro       15: "spawn102(hostlist,nserver|debug=1)\n\
        !            16: hostlist : [\"iyokan-0\",\"iyokan-1\",...]\n\
        !            17: nserver : number of servers to be used in hostlist\n\
1.1       noro       18: If debug is specified, the debug windows will appear."
                     19: {
                     20:        Debug = getopt(debug);
                     21:        if ( type(Debug) == -1 )
                     22:                Debug = 0;
                     23:        else if ( Debug )
                     24:                Debug = 1;
                     25:        Procs = newvect(N);
                     26:        for ( I = 0; I < N; I++ ) {
                     27:                if ( Debug )
                     28:                        Procs[I] = ox_launch(Hosts[I],get_rootdir(),"ox_asir");
                     29:                else
                     30:                        Procs[I] = ox_launch_nox(Hosts[I],get_rootdir(),"ox_asir");
                     31:                ox_set_rank_102(Procs[I],N,I);
                     32:                sleep(1000);
                     33:        }
                     34:        for ( I = 0; I < N; I++ )
                     35:                for ( J = I+1; J < N; J++ ) {
                     36:                        P = generate_port();
                     37:                        ox_tcp_accept_102(Procs[I],P,J);
                     38:                        ox_tcp_connect_102(Procs[J],Hosts[I],P,I);
1.2       noro       39:                }
                     40:        return Procs;
                     41: }
                     42:
                     43: def spawn102_local(N)
1.3     ! noro       44: "spawn102_local(nserver|debug=1)\n\
        !            45: nserver : number of servers to be used in hostlist\n\
1.2       noro       46: If debug is specified, the debug windows will appear."
                     47: {
                     48:        Debug = getopt(debug);
                     49:        if ( type(Debug) == -1 )
                     50:                Debug = 0;
                     51:        else if ( Debug )
                     52:                Debug = 1;
                     53:        Procs = newvect(N);
                     54:        for ( I = 0; I < N; I++ ) {
                     55:                if ( Debug )
                     56:                        Procs[I] = ox_launch();
                     57:                else
                     58:                        Procs[I] = ox_launch_nox();
                     59:                ox_set_rank_102(Procs[I],N,I);
                     60:                sleep(1000);
                     61:        }
                     62:        for ( I = 0; I < N; I++ )
                     63:                for ( J = I+1; J < N; J++ ) {
                     64:                        P = generate_port(1);
                     65:                        ox_tcp_accept_102(Procs[I],P,J);
                     66:                        ox_tcp_connect_102(Procs[J],0,P,I);
1.1       noro       67:                }
                     68:        return Procs;
                     69: }
                     70:
                     71: def urandompoly(N,D)
1.3     ! noro       72: "urandompoly(N,D)\n\
        !            73:  generate a univariate random polynomial of degree N,\n\
1.1       noro       74:  with D bit random coefficients."
                     75: {
                     76:        for(I=0,R=0;I<=N;I++)R+= lrandom(D)*x^I; return R;
                     77: }
                     78:
                     79: /*
                     80:        return: F1*F2
                     81:        if option 'proc' is supplied as a list of server id's,
                     82:        F1*F2 is calculated by distributed computation.
                     83: */
                     84:
                     85: def d_mul(F1,F2)
1.3     ! noro       86: "d_mul(F1,F2|proc=ProcList)\n\
        !            87:  computes the product of F1 and F2.\n\
1.1       noro       88:  If ProcList is specified, the product is computed in parallel."
                     89: {
                     90:        Procs = getopt(proc);
                     91:        if ( type(Procs) == -1 ) return umul(F1,F2);
                     92:        if ( !var(F1) || !var(F2) ) return F1*F2;
                     93:        NP = length(Procs);
                     94:        ox_push_cmo(Procs[0],[F1,F2]);
                     95:        for ( I = 0; I < NP; I++ )
                     96:                ox_cmo_rpc(Procs[I],"d_mul_main",0);
                     97:        R = ox_pop_cmo(Procs[0]);
                     98:        return R;
                     99: }
                    100:
                    101: def d_mul_main(Root)
                    102: {
                    103:        Id = ox_get_rank_102();
                    104:        NP = Id[0]; Rank = Id[1];
                    105:        Arg = ox_bcast_102(Root);
                    106:        F1 = Arg[0]; F2 = Arg[1];
                    107:        L = setup_modarrays(F1,F2,NP);
                    108:        Marray = L[0]; MIarray = L[1]; M = L[2];
                    109:        R = umul_chrem(F1,F2,MIarray[Rank],Marray[Rank],M);
                    110:        Arg = 0; F1 = 0; F2 = 0;
                    111:        R = ox_reduce_102(Root,"+",R);
                    112:        if ( Rank == Root )
                    113:                R = uadj_coef(R%M,M,ishift(M,1));
                    114:        return R;
                    115: }
                    116:
                    117: /*
                    118:  * Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
                    119:  * M = Marray[0]*...*Marray[NP-1]
                    120:  */
                    121:
                    122: def setup_modarrays(F1,F2,NP)
                    123: {
                    124:        V = var(F1);
                    125:        D1 = deg(F1,V); D2 = deg(F2,V);
                    126:        Dmin = MIN(D1,D2);
                    127:        Dfft = p_mag(D1+D2+1)+1;
                    128:        Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
                    129:        if ( Bound < 32 ) Bound = 32;
                    130:        Marray = newvect(NP); MIarray = newvect(NP);
                    131:        for ( I = 0; I < NP; I++ ) {
                    132:                Marray[I] = 1; MIarray[I] = [];
                    133:        }
                    134:        for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
                    135:                T = get_next_fft_prime(I,Dfft);
                    136:                if ( !T )
                    137:                        error("fft_mul_d : fft_prime exhausted.");
                    138:                Marray[J] *= T[1];
                    139:                MIarray[J] = cons(T[0],MIarray[J]);
                    140:                M *= T[1];
                    141:                I = T[0]+1;
                    142:        }
                    143:        return [Marray,MIarray,M];
                    144: }
                    145:
                    146: def umul_chrem(F1,F2,Ind,M1,M)
                    147: {
                    148:        T0 = time();
                    149:        C = umul_specialmod(F1,F2,Ind);
                    150:        Mhat = idiv(M,M1);
                    151:        MhatInv = inv(Mhat,M1);
                    152:        R = Mhat*((MhatInv*C)%M1);
                    153:        return R;
                    154: }
                    155: end$