Return to dmul102 CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/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)
        !            16: hostlist : [\"iyokan-0\",\"iyokan-1\",...]
        !            17: nserver : number of servers to be used in hostlist
        !            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 urandompoly(N,D)
        !            44: "urandompoly(N,D)
        !            45:  generate a univariate random polynomial of degree N,
        !            46:  with D bit random coefficients."
        !            47: {
        !            48:        for(I=0,R=0;I<=N;I++)R+= lrandom(D)*x^I; return R;
        !            49: }
        !            50:
        !            51: /*
        !            52:        return: F1*F2
        !            53:        if option 'proc' is supplied as a list of server id's,
        !            54:        F1*F2 is calculated by distributed computation.
        !            55: */
        !            56:
        !            57: def d_mul(F1,F2)
        !            58: "d_mul(F1,F2|proc=ProcList)
        !            59:  computes the product of F1 and F2.
        !            60:  If ProcList is specified, the product is computed in parallel."
        !            61: {
        !            62:        Procs = getopt(proc);
        !            63:        if ( type(Procs) == -1 ) return umul(F1,F2);
        !            64:        if ( !var(F1) || !var(F2) ) return F1*F2;
        !            65:        NP = length(Procs);
        !            66:        ox_push_cmo(Procs[0],[F1,F2]);
        !            67:        for ( I = 0; I < NP; I++ )
        !            68:                ox_cmo_rpc(Procs[I],"d_mul_main",0);
        !            69:        R = ox_pop_cmo(Procs[0]);
        !            70:        return R;
        !            71: }
        !            72:
        !            73: def d_mul_main(Root)
        !            74: {
        !            75:        Id = ox_get_rank_102();
        !            76:        NP = Id[0]; Rank = Id[1];
        !            77:        Arg = ox_bcast_102(Root);
        !            78:        F1 = Arg[0]; F2 = Arg[1];
        !            79:        L = setup_modarrays(F1,F2,NP);
        !            80:        Marray = L[0]; MIarray = L[1]; M = L[2];
        !            81:        R = umul_chrem(F1,F2,MIarray[Rank],Marray[Rank],M);
        !            82:        Arg = 0; F1 = 0; F2 = 0;
        !            83:        R = ox_reduce_102(Root,"+",R);
        !            84:        if ( Rank == Root )
        !            85:                R = uadj_coef(R%M,M,ishift(M,1));
        !            86:        return R;
        !            87: }
        !            88:
        !            89: /*
        !            90:  * Marray[J] = FFTprime[Marray[J][0]]*...*FFTprime[Marray[J][...]]
        !            91:  * M = Marray[0]*...*Marray[NP-1]
        !            92:  */
        !            93:
        !            94: def setup_modarrays(F1,F2,NP)
        !            95: {
        !            96:        V = var(F1);
        !            97:        D1 = deg(F1,V); D2 = deg(F2,V);
        !            98:        Dmin = MIN(D1,D2);
        !            99:        Dfft = p_mag(D1+D2+1)+1;
        !           100:        Bound = maxblen(F1)+maxblen(F2)+p_mag(Dmin)+1;
        !           101:        if ( Bound < 32 ) Bound = 32;
        !           102:        Marray = newvect(NP); MIarray = newvect(NP);
        !           103:        for ( I = 0; I < NP; I++ ) {
        !           104:                Marray[I] = 1; MIarray[I] = [];
        !           105:        }
        !           106:        for ( M = 1, I = 0, J = 0; p_mag(M) <= Bound; J = (J+1)%NP ) {
        !           107:                T = get_next_fft_prime(I,Dfft);
        !           108:                if ( !T )
        !           109:                        error("fft_mul_d : fft_prime exhausted.");
        !           110:                Marray[J] *= T[1];
        !           111:                MIarray[J] = cons(T[0],MIarray[J]);
        !           112:                M *= T[1];
        !           113:                I = T[0]+1;
        !           114:        }
        !           115:        return [Marray,MIarray,M];
        !           116: }
        !           117:
        !           118: def umul_chrem(F1,F2,Ind,M1,M)
        !           119: {
        !           120:        T0 = time();
        !           121:        C = umul_specialmod(F1,F2,Ind);
        !           122:        Mhat = idiv(M,M1);
        !           123:        MhatInv = inv(Mhat,M1);
        !           124:        R = Mhat*((MhatInv*C)%M1);
        !           125:        return R;
        !           126: }
        !           127: end$