Return to dmul102 CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / lib

Annotation of OpenXM_contrib2/asir2018/lib/dmul102, Revision 1.1

1.1     ! noro        1: /* $OpenXM$ */
        !             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)
        !            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\
        !            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);
        !            39:                }
        !            40:        return Procs;
        !            41: }
        !            42:
        !            43: def spawn102_local(N)
        !            44: "spawn102_local(nserver|debug=1)\n\
        !            45: nserver : number of servers to be used in hostlist\n\
        !            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);
        !            67:                }
        !            68:        return Procs;
        !            69: }
        !            70:
        !            71: def urandompoly(N,D)
        !            72: "urandompoly(N,D)\n\
        !            73:  generate a univariate random polynomial of degree N,\n\
        !            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)
        !            86: "d_mul(F1,F2|proc=ProcList)\n\
        !            87:  computes the product of F1 and F2.\n\
        !            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$