[BACK]Return to nf CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / lib

Annotation of OpenXM_contrib2/asir2000/lib/nf, Revision 1.1.1.1

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/lib/nf,v 1.1.1.1 1999/11/10 08:12:31 noro Exp $ */
                      2: extern IDIVV_REM,IDIVV_HIST;
                      3:
                      4: def nf_mod(B,G,MOD,DIR,PS)
                      5: {
                      6:        setmod(MOD);
                      7:        for ( D = 0, H = []; G; ) {
                      8:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                      9:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                     10:                                R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
                     11:                                CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
                     12:                                U = G-CR*R;
                     13:                                print(".",2);
                     14:                                if ( !U )
                     15:                                        return D;
                     16:                                break;
                     17:                        }
                     18:                }
                     19:                if ( U )
                     20:                        G = U;
                     21:                else {
                     22:                        D += dp_hm(G); G = dp_rest(G); print(!D,2);
                     23:                }
                     24:        }
                     25:        return D;
                     26: }
                     27:
                     28: def nf_mod_ext(B,G,MOD,DIR,PS)
                     29: {
                     30:        setmod(MOD);
                     31:        for ( D = 0, H = []; G; ) {
                     32:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                     33:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                     34:                                R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
                     35:                                CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
                     36:                                U = G-CR*R;
                     37:                                H = cons([CR,strtov("t"+rtostr(car(L)))],H);
                     38:                                print([length(H)]);
                     39:                                if ( !U )
                     40:                                        return [D,reverse(H)];
                     41:                                break;
                     42:                        }
                     43:                }
                     44:                if ( U )
                     45:                        G = U;
                     46:                else {
                     47:                        D += dp_hm(G); G = dp_rest(G);
                     48:                }
                     49:        }
                     50:        return [D,reverse(H)];
                     51: }
                     52:
                     53: def nf(B,G,M,PS)
                     54: {
                     55:        for ( D = 0, H = []; G; ) {
                     56:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                     57:                        if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
                     58:                                GCD = igcd(dp_hc(G),dp_hc(R));
                     59:                                CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
                     60:                                U = CG*G-dp_subd(G,R)*CR*R;
                     61:                                H = cons(car(L),H);
                     62:                                if ( !U )
                     63:                                        return [D,M,reverse(H)];
                     64:                                D *= CG; M *= CG;
                     65:                                break;
                     66:                        }
                     67:                }
                     68:                if ( U )
                     69:                        G = U;
                     70:                else {
                     71:                        D += dp_hm(G); G = dp_rest(G);
                     72:                }
                     73:        }
                     74:        return [D,M,reverse(H)];
                     75: }
                     76: extern M_NM,M_DN;
                     77:
                     78: def nf_demand(B,G,DIR,PSH)
                     79: {
                     80:        T1 = T2 = T3 = T4 = 0;
                     81:        Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
                     82:        for ( D = 0, H = []; G; ) {
                     83:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                     84:                        if ( dp_redble(G,PSH[car(L)]) > 0 ) {
                     85:                                R = bload(DIR+rtostr(car(L)));
                     86:                                H = cons(car(L),H);
                     87:                                print([length(H),dp_mag(dp_hm(G))]);
                     88:                                GCD = igcd(dp_hc(G),dp_hc(R));
                     89:                                CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
                     90:                                T0 = time()[0]; A1 = CG*G; T1 += time()[0]-T0;
                     91:                                T0 = time()[0]; A2 = dp_subd(G,R)*CR*R; T2 += time()[0]-T0;
                     92:                                T0 = time()[0]; U = A1-A2; T3 += time()[0]-T0;
                     93:                                if ( !U )
                     94:                                        return [D,reverse(H),T1,T2,T3,T4];
                     95:                                D *= CG;
                     96:                                break;
                     97:                        }
                     98:                }
                     99:                if ( U ) {
                    100:                        G = U;
                    101:                        if ( D ) {
                    102:                                if ( dp_mag(dp_hm(D)) > Mag ) {
                    103:                                        T0 = time()[0];
                    104:                                        print("ptozp2");
                    105:                                        T = dp_ptozp2(D,G); D = T[0]; G = T[1];
                    106:                                        T4 += time()[0]-T0;
                    107:                                        Mag = idiv(dp_mag(dp_hm(D))*M_NM,M_DN);
                    108:                                }
                    109:                        } else {
                    110:                                if ( dp_mag(dp_hm(G)) > Mag ) {
                    111:                                        T0 = time()[0];
                    112:                                        print("ptozp");
                    113:                                        G = dp_ptozp(G);
                    114:                                        T4 += time()[0]-T0;
                    115:                                        Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
                    116:                                }
                    117:                        }
                    118:                } else {
                    119:                        D += dp_hm(G); G = dp_rest(G);
                    120:                }
                    121:        }
                    122:        return [D,reverse(H),T1,T2,T3,T4];
                    123: }
                    124:
                    125: def nf_demand_d(B,G,DIR,NDIR,PDIR,PSH,PROC)
                    126: {
                    127:        INDEX = 0;
                    128:        T1 = T2 = 0;
                    129: /*     dp_gr_flags(["Dist",PROC]); */
                    130:        if ( PROC != [] ) {
                    131:                PROC = newvect(length(PROC),PROC);
                    132:                NPROC = size(PROC)[0];
                    133:        } else
                    134:                NPROC = 0;
                    135:        Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
                    136:        Kara = dp_set_kara()*27; /* XXX */
                    137:        D = [0,0]; R = [1,G];
                    138:        print("");
                    139:        for ( H = []; R[1]; ) {
                    140:                for ( U = 0, L = B; L != []; L = cdr(L) ) {
                    141:                        if ( dp_redble(R[1],PSH[car(L)]) > 0 ) {
                    142:                                Red = bload(DIR+rtostr(car(L)));
                    143:                                H = cons(car(L),H);
                    144: /*                             D0=dp_mag(D[0]*<<0>>); D1=dp_mag(dp_hm(D[1]));
                    145:                                R0=dp_mag(R[0]*<<0>>); R1=dp_mag(dp_hm(R[1]));
                    146:                                print([car(L),length(H),[D0,D1,dp_nt(D[1])],[R0,R1,dp_nt(R[1])]],2); */
                    147:
                    148:                                GCD = igcd(dp_hc(R[1]),dp_hc(Red));
                    149:                                CR = idiv(dp_hc(Red),GCD);
                    150:                                CRed = idiv(dp_hc(R[1]),GCD);
                    151:
                    152:                                T0 = time()[3];
                    153:                                if ( (PROC != []) && (dp_mag(dp_hm(Red)) > Kara) ) {
                    154:                                        print("d",2);
                    155:                                        rpc(PROC[0],"dp_imul_index",CRed,car(L));
                    156:                                        U = CR*R[1] - dp_subd(R[1],Red)*rpcrecv(PROC[0]);
                    157:                                } else {
                    158:                                        print("l",2);
                    159:                                        U = CR*R[1] - dp_subd(R[1],Red)*CRed*Red;
                    160:                                }
                    161:                                TT = time()[3]-T0; T1 += TT; /* print(TT); */
                    162:                                if ( !U )
                    163:                                        return [D,reverse(H),T1,T2];
                    164:                                break;
                    165:                        }
                    166:                }
                    167:                if ( U ) {
                    168:                        if ( (HMAG=dp_mag(dp_hm(U))) > Mag ) {
                    169:                                T0 = time()[3];
                    170:                                if ( HMAG > Kara ) {
                    171:                                        print("D",2);
                    172:                                        T = dp_ptozp_d(U,PROC,NPROC);
                    173:                                } else {
                    174:                                        print("L",2);
                    175:                                        T = dp_ptozp(U);
                    176:                                }
                    177:                                TT = time()[3]-T0; T2 += TT; /* print(TT); */
                    178:                                Cont = idiv(dp_hc(U),dp_hc(T));
                    179:                                D0 = kmul(CR,D[0]);
                    180:                                R0 = kmul(Cont,R[0]);
                    181:                                S = igcd(D0,R0);
                    182:                                D = [idiv(D0,S),D[1]];
                    183:                                R = [idiv(R0,S),T];
                    184:                                Mag = idiv(dp_mag(dp_hm(R[1]))*M_NM,M_DN);
                    185:                        } else {
                    186:                                D = [kmul(CR,D[0]),D[1]];
                    187:                                R = [R[0],U];
                    188:                        }
                    189:                } else {
                    190:                        C = kmul(dp_hc(R[1]),R[0]);
                    191:                        T = igcd(D[0],C);
                    192:                        D = [T,idiv(D[0],T)*D[1]+idiv(C,T)*dp_ht(R[1])];
                    193:                        R = [R[0],dp_rest(R[1])];
                    194:                }
                    195:        }
                    196:        return [D[1],reverse(H),T1,T2];
                    197: }
                    198:
                    199: extern PTOZP,DPCV,NFSDIR;
                    200:
                    201: def imulv(A,B)
                    202: {
                    203:        return A*B;
                    204: }
                    205:
                    206: def dp_imul(P,N)
                    207: {
                    208:        return N*P;
                    209: }
                    210:
                    211: def mul(A,B)
                    212: {
                    213:        return A*B;
                    214: }
                    215:
                    216: def dp_imul_index(C,INDEX)
                    217: {
                    218:        T = C*bload(NFSDIR+rtostr(INDEX));
                    219:        return T;
                    220: }
                    221:
                    222: def reg_dp_dtov(P)
                    223: {
                    224:        PTOZP = P;
                    225:        DPCV = dp_dtov(PTOZP);
                    226: }
                    227:
                    228: def reg_dp_iqr(G)
                    229: {
                    230:        N = size(DPCV)[0];
                    231:        for ( R = [], I = 0; I < N; I++ ) {
                    232:                DPCV[I] = T = irem(DPCV[I],G);
                    233:                if ( T )
                    234:                        R = cons(T,R);
                    235:        }
                    236:        return R;
                    237: }
                    238:
                    239: def reg_dp_cont(P)
                    240: {
                    241:        PTOZP = P;
                    242:        C = dp_dtov(P);
                    243:        return igcd(C);
                    244: }
                    245:
                    246: def reg_dp_idiv(GCD)
                    247: {
                    248:        return dp_idiv(PTOZP,GCD);
                    249: }
                    250:
                    251: def dp_cont_d(G,PROC,NPROC)
                    252: {
                    253:        C = dp_sep(G,NPROC);
                    254:        N = size(C)[0];
                    255:        for ( I = 0; I < N; I++ )
                    256:                rpc(PROC[I],"reg_dp_cont",C[I]);
                    257:        R = map(rpcrecv,PROC);
                    258:        GCD = igcd(R);
                    259:        return GCD;
                    260: }
                    261:
                    262: /*
                    263: def iqrv(C,D)
                    264: {
                    265:        return map(iqr,C,D);
                    266: }
                    267: */
                    268:
                    269: #if 1
                    270: def dp_ptozp_d(G,PROC,NPROC)
                    271: {
                    272:        C = dp_dtov(G); L = size(C)[0];
                    273:        T0 = time()[3];
                    274:        D0 = igcd_estimate(C);
                    275:        TT = time()[3]-T0; TG1 += TT;
                    276:        CS = sepvect(C,NPROC+1);
                    277:        N = size(CS)[0]; N1 = N-1;
                    278:        QR = newvect(N); Q = newvect(L); R = newvect(L);
                    279:        T0 = time()[3];
                    280:        for ( I = 0; I < N1; I++ )
                    281:                rpc0(PROC[I],"iqrv",CS[I],D0);
                    282:        QR[N1] = iqrv(CS[N1],D0);
                    283:        for ( I = 0; I < N1; I++ )
                    284:                QR[I] = rpcrecv(PROC[I]);
                    285:        TT = time()[3]-T0; TD += TT;
                    286:        for ( I = J = 0; I < N; I++ ) {
                    287:                T = QR[I]; M = size(T)[0];
                    288:                for ( K = 0; K < M; K++, J++ ) {
                    289:                        Q[J] = T[K][0]; R[J] = T[K][1];
                    290:                }
                    291:        }
                    292:        T0 = time()[3];
                    293:        D1 = igcd(R);
                    294:        GCD = igcd(D0,D1);
                    295:        A = idiv(D0,GCD);
                    296:        TT = time()[3]-T0; TG2 += TT;
                    297:        T0 = time()[3];
                    298:        for ( I = 0; I < L; I++ )
                    299:                Q[I] = kmul(A,Q[I])+idiv(R[I],GCD);
                    300:        TT = time()[3]-T0; TM += TT;
                    301:        print([TG1,TD,TG2,TM,dp_mag(D0*<<0>>),dp_mag(D1*<<0>>),dp_mag(GCD*<<0>>)],2);
                    302:        return dp_vtod(Q,G);
                    303: }
                    304: #endif
                    305:
                    306: #if 0
                    307: def dp_ptozp_d(G,PROC,NPROC)
                    308: {
                    309:        C = dp_sep(G,NPROC+1);
                    310:        N = size(C)[0]; N1 = N-1;
                    311:        R = newvect(N);
                    312:        for ( I = 0; I < N1; I++ )
                    313:                rpc(PROC[I],"reg_igcdv",dp_dtov(C[I]));
                    314:        R[N1] = reg_igcdv(dp_dtov(C[N1]));
                    315:        for ( I = 0; I < N1; I++ )
                    316:                R[I] = rpcrecv(PROC[I]);
                    317:        GCD = igcd(R);
                    318:        for ( I = 0; I < N1; I++ )
                    319:                rpc(PROC[I],"idivv_final",GCD);
                    320:        S = dp_vtod(idivv_final(GCD),C[N1]);
                    321:        for ( I = 0; I < N1; I++ )
                    322:                S += dp_vtod(rpcrecv(PROC[I]),C[I]);
                    323:        return S;
                    324: }
                    325: #endif
                    326:
                    327: #if 0
                    328: def dp_ptozp_d(G,PROC,NPROC)
                    329: {
                    330:        C = dp_sep(G,NPROC+1);
                    331:        N = size(C)[0]; N1 = N-1;
                    332:        R = newvect(N);
                    333:        for ( I = 0; I < N1; I++ )
                    334:                rpc(PROC[I],"reg_dp_cont",C[I]);
                    335:        R[N1] = igcd(dp_dtov(C[N1]));
                    336:        for ( I = 0; I < N1; I++ )
                    337:                R[I] = rpcrecv(PROC[I]);
                    338:        GCD = igcd(R);
                    339:        for ( I = 0; I < N1; I++ )
                    340:                rpc(PROC[I],"reg_dp_idiv",GCD);
                    341:        S = dp_idiv(C[N1],GCD);
                    342:        for ( I = 0; I < N1; I++ )
                    343:                S += rpcrecv(PROC[I]);
                    344:        return S;
                    345: }
                    346: #endif
                    347:
                    348: #if 0
                    349: def dp_ptozp_d(G,PROC,NPROC)
                    350: {
                    351:        C = dp_sep(G,NPROC);
                    352:        N = size(C)[0]; T = newvect(N);
                    353:        for ( I = 0; I < N; I++ )
                    354:                T[I] = PROC[I];
                    355:        PROC = T;
                    356:        for ( I = 0; I < N; I++ )
                    357:                rpc(PROC[I],"reg_dp_dtov",C[I]);
                    358:        A = dp_dtov(G); A = isort(A); L = size(A)[0];
                    359:        map(rpcrecv,PROC);
                    360:        while ( 1 ) {
                    361:                GCD = igcd(A[0],A[1]);
                    362:                for ( I = 2; I < L; I++ ) {
                    363:                        GT = igcd(GCD,A[I]);
                    364:                        if ( GCD == GT )
                    365:                                break;
                    366:                        else
                    367:                                GCD = GT;
                    368:                }
                    369:                if ( I == L )
                    370:                        break;
                    371:                else {
                    372:                        map(rpc,PROC,"reg_dp_iqr",GCD);
                    373:                        R = map(rpcrecv,PROC);
                    374:                        for ( J = 0, U = []; J < N; J++ )
                    375:                                U = append(R[J],U);
                    376:                        L = length(U);
                    377:                        if ( L == 0 )
                    378:                                break;
                    379:                        else
                    380:                                U = cons(GCD,U);
                    381:                        A = isort(newvect(L+1,U));
                    382:                        print(["length=",L,"minmag=",dp_mag(A[0]*<<0>>)]);
                    383:                }
                    384:        }
                    385:        for ( I = 0; I < N; I++ )
                    386:                rpc(PROC[I],"reg_dp_idiv",GCD);
                    387:        R = map(rpcrecv,PROC);
                    388:        for ( I = 0, S = 0; I < N; I++ )
                    389:                S += R[I];
                    390:        return S;
                    391: }
                    392: #endif
                    393:
                    394: def dp_ptozp2_d(D,G,PROC,NPROC)
                    395: {
                    396:        T = dp_ptozp_d(D+G,PROC,NPROC);
                    397:        for ( S = 0; D; D = dp_rest(D), T = dp_rest(T) )
                    398:                S += dp_hm(T);
                    399:        return [S,T];
                    400: }
                    401:
                    402: def genindex(N)
                    403: {
                    404:        R = [];
                    405:        for ( I = 0; I < N; I++ )
                    406:                R = cons(I,R);
                    407:        return reverse(R);
                    408: }
                    409:
                    410: def nftest_ext_mod(N1,N2,N,MOD,DIR,PSH)
                    411: {
                    412:        S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
                    413:        R = nf_mod_ext(genindex(N-1),S,MOD,DIR,PSH);
                    414:        return R;
                    415: }
                    416:
                    417: def nftest_mod(N1,N2,N,MOD,DIR,PSH)
                    418: {
                    419:        S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
                    420:        R = nf_mod(genindex(N-1),S,MOD,DIR,PSH);
                    421:        return dp_mdtod(R);
                    422: }
                    423:
                    424: def nftest(N1,N2,N,DIR,PSH)
                    425: {
                    426:        S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
                    427:        R = nf_demand(genindex(N-1),S,DIR,PSH);
                    428:        return R;
                    429: }
                    430:
                    431: def nftest_d(N1,N2,N,DIR,NDIR,PDIR,PSH,PROC)
                    432: {
                    433:        S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
                    434:        R = nf_demand_d(genindex(N-1),S,DIR,NDIR,PDIR,PSH,PROC);
                    435:        return R;
                    436: }
                    437:
                    438: def abs(A)
                    439: {
                    440:        return A >= 0 ? A : -A;
                    441: }
                    442:
                    443: def sort(A)
                    444: {
                    445:        N = size(A)[0];
                    446:        for ( I = 0; I < N; I++ ) {
                    447:                for ( M = I, J = I + 1; J < N; J++ )
                    448:                        if ( A[J] < A[M] )
                    449:                                M = J;
                    450:                if ( I != M ) {
                    451:                        T = A[I]; A[I] = A[M]; A[M] = T;
                    452:                }
                    453:        }
                    454:        return A;
                    455: }
                    456:
                    457: def igcdv(C)
                    458: {
                    459:        C = sort(C); N = size(C)[0];
                    460:        G = igcd(C[0],C[1]);
                    461:        for ( I = 2; I < N; I++ ) {
                    462:                T = time();
                    463:                G = igcd(G,C[I]);
                    464:        /*      print([I,dp_mag(G*<<0>>),time()[0]-T[0]]); */
                    465:        }
                    466:        return G;
                    467: }
                    468:
                    469: def igcdv2(C)
                    470: {
                    471:        C = sort(C); N = size(C)[0];
                    472:        G = igcd(C[0],C[1]);
                    473:        for ( I = 2; I < N; I++ ) {
                    474:                T = time();
                    475:                C[I] %= G;
                    476:                GT = igcd(G,C[I]);
                    477:                if ( G == GT ) {
                    478:                        for ( J = I+1, NZ=0; J < N; J++ ) {
                    479:                                C[J] %= G;
                    480:                                if ( C[J] )
                    481:                                        NZ = 1;
                    482:                        }
                    483:                        if ( !NZ )
                    484:                                break;
                    485:                } else
                    486:                        G = GT;
                    487:                print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
                    488:        }
                    489:        return G;
                    490: }
                    491:
                    492: def igcdv_d(C,P,NP)
                    493: {
                    494:        C = isort(C); N = size(C)[0];
                    495:        G = igcd(C[0],C[1]);
                    496:        for ( I = 2; I < N; I++ ) {
                    497:                T = time();
                    498:                C[I] %= G;
                    499:                GT = igcd(G,C[I]);
                    500:                if ( G == GT ) {
                    501:                        for ( J = I+1, NZ=0; J < N; J++ ) {
                    502:                                C[J] %= G;
                    503:                                if ( C[J] )
                    504:                                        NZ = 1;
                    505:                        }
                    506:                        if ( !NZ )
                    507:                                break;
                    508:                } else
                    509:                        G = GT;
                    510:                print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
                    511:        }
                    512:        return G;
                    513: }
                    514:
                    515: def dup(A)
                    516: {
                    517:        N = size(A)[0]; V = newvect(N);
                    518:        for ( I = 0; I < N; I++ )
                    519:                V[I] = A[I];
                    520:        return V;
                    521: }
                    522:
                    523: def idivv_init(A)
                    524: {
                    525:        IDIVV_REM = dup(A);
                    526: }
                    527:
                    528: def igcd_cofactor(D)
                    529: {
                    530:        T = map(iqr,IDIVV_REM,D);
                    531:        N = size(T)[0];
                    532:        Q = newvect(N); R = newvect(N);
                    533:        for ( I = 0; I < N; I++ ) {
                    534:                Q[I] = T[I][0]; R[I] = T[I][1];
                    535:        }
                    536:        D1 = igcdv(dup(R));
                    537:        if ( !D1 ) {
                    538:                IDIVV_HIST = [[Q,D]];
                    539:                return D;
                    540:        } else {
                    541:                Q1 = map(idiv,R,D1);
                    542:                IDIVV_HIST = [[Q,D],[Q1,D1]];
                    543:                return D1;
                    544:        }
                    545: }
                    546:
                    547: def idivv_step(D)
                    548: {
                    549:        T = map(iqr,IDIVV_REM,D);
                    550:        N = size(T)[0];
                    551:        Q = newvect(N); R = newvect(N);
                    552:        for ( I = 0, NZ = 0; I < N; I++ ) {
                    553:                Q[I] = T[I][0]; R[I] = T[I][1];
                    554:                if ( R[I] )
                    555:                        NZ = 1;
                    556:        }
                    557:        IDIVV_REM = R;
                    558:        IDIVV_HIST = cons([Q,D],IDIVV_HIST);
                    559:        return NZ;
                    560: }
                    561:
                    562: def igcdv3(C)
                    563: {
                    564:        idivv_init(C);
                    565:        C = isort(C); N = size(C)[0];
                    566:        D = igcd(C[0],C[1]);
                    567:        for ( I = 2; I < N; I++ ) {
                    568:                T = time();
                    569:                C[I] %= G;
                    570:                GT = igcd(G,C[I]);
                    571:                if ( G == GT ) {
                    572:                        NZ = idivv_step(G);
                    573:                        if ( !NZ )
                    574:                                break;
                    575:                } else
                    576:                        G = GT;
                    577:        }
                    578:        CF = idivv_final(G);
                    579:        return [G,CF];
                    580: }
                    581:
                    582: def igcdv4(C)
                    583: {
                    584:        idivv_init(C);
                    585:        C = isort(C); N = size(C)[0];
                    586:        D = igcd_estimate(C);
                    587:        G = igcd_cofactor(D);
                    588:        G = igcd(G,D);
                    589:        CF = idivv_final(G);
                    590:        return [G,CF];
                    591: }
                    592:
                    593: def igcd_estimate(C)
                    594: {
                    595:        C = isort(C); N = size(C)[0];
                    596:        M = idiv(N,2);
                    597:        for ( I = 0, S = 0; I < M; I++ )
                    598:                S += C[I]>0?C[I]:-C[I];
                    599:        for ( T = 0; I < N; I++ )
                    600:                T += C[I]>0?C[I]:-C[I];
                    601:        if ( !S && !T )
                    602:                return igcdv(C);
                    603:        else
                    604:                return igcd(S,T);
                    605: }
                    606:
                    607: def reg_igcdv(C)
                    608: {
                    609:        D = igcd_estimate(C);
                    610:        T = map(iqr,C,D); N = size(T)[0];
                    611:        Q = newvect(N); R = newvect(N);
                    612:        for ( I = 0; I < N; I++ ) {
                    613:                Q[I] = T[I][0]; R[I] = T[I][1];
                    614:        }
                    615:        D1 = igcdv(dup(R));
                    616:        if ( !D1 ) {
                    617:                IDIVV_HIST = [[Q,D]];
                    618:                return D;
                    619:        } else {
                    620:                Q1 = map(idiv,R,D1);
                    621:                IDIVV_HIST = [[Q,D],[Q1,D1]];
                    622:                return igcd(D,D1);
                    623:        }
                    624: }
                    625:
                    626: def idivv_final(D)
                    627: {
                    628:        for ( T = IDIVV_HIST, S = 0; T != []; T = cdr(T) ) {
                    629:                H = car(T); S += map(kmul,H[0],idiv(H[1],D));
                    630:        }
                    631:        return S;
                    632: }
                    633:
                    634: def dp_vtod(C,P)
                    635: {
                    636:        for ( R = 0, I = 0; P; P = dp_rest(P), I++ )
                    637:                R += C[I]*dp_ht(P);
                    638:        return R;
                    639: }
                    640:
                    641: def dp_nt(P)
                    642: {
                    643:        for ( I = 0; P; P = dp_rest(P), I++ );
                    644:        return I;
                    645: }
                    646: end$

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>