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

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

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/lib/fff,v 1.1.1.1 1999/11/10 08:12:31 noro Exp $ */
                      2: /*
                      3:        fff : Univariate factorizer over a finite field.
                      4:
                      5:     Revision History:
                      6:
                      7:     99/05/18    noro    the first official version
                      8:     99/06/11    noro
                      9:     99/07/27    noro
                     10: */
                     11:
                     12: #include "defs.h"
                     13:
                     14: extern TPMOD,TQMOD$
                     15:
                     16: /*
                     17:   Input : a univariate polynomial F
                     18:   Output: a list [[F1,M1],[F2,M2],...], where
                     19:           Fi is a monic irreducible factor, Mi is its multiplicity.
                     20:           The leading coefficient of F is abondoned.
                     21: */
                     22:
                     23: def fctr_ff(F)
                     24: {
                     25:        F = simp_ff(F);
                     26:        F = F/LCOEF(F);
                     27:        L = sqfr_ff(F);
                     28:        for ( R = [], T = L; T != []; T = cdr(T) ) {
                     29:                S = car(T); A = S[0]; E = S[1];
                     30:                B = ddd_ff(A);
                     31:                R = append(append_mult_ff(B,E),R);
                     32:        }
                     33:        return R;
                     34: }
                     35:
                     36: /*
                     37:   Input : a list of polynomial L; an integer E
                     38:   Output: a list s.t. [[L0,E],[L1,E],...]
                     39:           where Li = L[i]/leading coef of L[i]
                     40: */
                     41:
                     42: def append_mult_ff(L,E)
                     43: {
                     44:        for ( T = L, R = []; T != []; T = cdr(T) )
                     45:                R = cons([car(T)/LCOEF(car(T)),E],R);
                     46:        return R;
                     47: }
                     48:
                     49: /*
                     50:        Input : a polynomial F
                     51:        Output: a list [[F1,M1],[F2,M2],...]
                     52:                where Fi is a square free factor,
                     53:                Mi is its multiplicity.
                     54: */
                     55:
                     56: def sqfr_ff(F)
                     57: {
                     58:        V = var(F);
                     59:        F1 = diff(F,V);
                     60:        L = [];
                     61:        /* F=H*Fq^p => F'=H'*Fq^p => gcd(F,F')=gcd(H,H')*Fq^p */
                     62:        if ( F1 != 0 ) {
                     63:                F1 = F1/LCOEF(F1);
                     64:                F2 = ugcd(F,F1);
                     65:                /* FLAT = H/gcd(H,H') : square free part of H */
                     66:                FLAT = sdiv(F,F2);
                     67:                I = 0;
                     68:                /* square free factorization of H */
                     69:                while ( deg(FLAT,V) ) {
                     70:                        while ( 1 ) {
                     71:                                QR = sqr(F,FLAT);
                     72:                                if ( !QR[1] ) {
                     73:                                        F = QR[0]; I++;
                     74:                                } else
                     75:                                        break;
                     76:                        }
                     77:                        if ( !deg(F,V) )
                     78:                                FLAT1 = simp_ff(1);
                     79:                        else
                     80:                                FLAT1 = ugcd(F,FLAT);
                     81:                        G = sdiv(FLAT,FLAT1);
                     82:                        FLAT = FLAT1;
                     83:                        L = cons([G,I],L);
                     84:                }
                     85:        }
                     86:        /* now F = Fq^p */
                     87:        if ( deg(F,V) ) {
                     88:                Char = characteristic_ff();
                     89:                T = sqfr_ff(pthroot_p_ff(F));
                     90:                for ( R = []; T != []; T = cdr(T) ) {
                     91:                        H = car(T); R = cons([H[0],Char*H[1]],R);
                     92:                }
                     93:        } else
                     94:                R = [];
                     95:        return append(L,R);
                     96: }
                     97:
                     98: /*
                     99:        Input : a polynomial F
                    100:        Output: F^(1/char)
                    101: */
                    102:
                    103: def pthroot_p_ff(F)
                    104: {
                    105:        V = var(F);
                    106:        DVR = characteristic_ff();
                    107:        PWR = DVR^(extdeg_ff()-1);
                    108:        for ( T = F, R = 0; T; ) {
                    109:                D1 = deg(T,V); C = coef(T,D1,V); T -= C*V^D1;
                    110:                R += C^PWR*V^idiv(D1,DVR);
                    111:        }
                    112:        return R;
                    113: }
                    114:
                    115: /*
                    116:        Input : a polynomial F of degree N
                    117:        Output: a list [V^Ord mod F,Tab]
                    118:                where V = var(F), Ord = field order
                    119:                Tab[i] = V^(i*Ord) mod F (i=0,...,N-1)
                    120: */
                    121:
                    122: def tab_ff(F)
                    123: {
                    124:        V = var(F);
                    125:        N = deg(F,V);
                    126:        F = F/LCOEF(F);
                    127:        XP = pwrmod_ff(F);
                    128:        R = pwrtab_ff(F,XP);
                    129:        return R;
                    130: }
                    131:
                    132: /*
                    133:        Input : a square free polynomial F
                    134:        Output: a list [F1,F2,...]
                    135:                where Fi is an irreducible factor of F.
                    136: */
                    137:
                    138: def ddd_ff(F)
                    139: {
                    140:        V = var(F);
                    141:        if ( deg(F,V) == 1 )
                    142:                return [F];
                    143:        TAB = tab_ff(F);
                    144:        for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
                    145:                lazy_lm(1);
                    146:                for ( T = 0, K = 0; K <= deg(W,V); K++ )
                    147:                        if ( C = coef(W,K,V) )
                    148:                                T += TAB[K]*C;
                    149:                lazy_lm(0);
                    150:                W = simp_ff(T);
                    151:                GCD = ugcd(F,W-V);
                    152:                if ( deg(GCD,V) ) {
                    153:                        L = append(berlekamp_ff(GCD,I,TAB),L);
                    154:                        F = sdiv(F,GCD);
                    155:                        W = urem(W,F);
                    156:                }
                    157:        }
                    158:        if ( deg(F,V) )
                    159:                return cons(F,L);
                    160:        else
                    161:                return L;
                    162: }
                    163:
                    164: /*
                    165:        Input : a polynomial
                    166:        Output: 1 if F is irreducible
                    167:                        0 otherwise
                    168: */
                    169:
                    170: def irredcheck_ff(F)
                    171: {
                    172:        V = var(F);
                    173:        if ( deg(F,V) <= 1 )
                    174:                return 1;
                    175:        F1 = diff(F,V);
                    176:        if ( !F1 )
                    177:                return 0;
                    178:        F1 = F1/LCOEF(F1);
                    179:        if ( deg(ugcd(F,F1),V) > 0 )
                    180:                return 0;
                    181:        TAB = tab_ff(F);
                    182:        for ( I = 1, W = V, L = []; 2*I <= deg(F,V); I++ ) {
                    183:                for ( T = 0, K = 0; K <= deg(W,V); K++ )
                    184:                        if ( C = coef(W,K,V) )
                    185:                                T += TAB[K]*C;
                    186:                W = T;
                    187:                GCD = ugcd(F,W-V);
                    188:                if ( deg(GCD,V) )
                    189:                        return 0;
                    190:        }
                    191:        return 1;
                    192: }
                    193:
                    194: /*
                    195:        Input : a square free (canonical) modular polynomial F
                    196:        Output: a list of polynomials [LF,CF,XP] where
                    197:                LF=the product of all the linear factors
                    198:                CF=F/LF
                    199:                XP=x^field_order mod CF
                    200: */
                    201:
                    202: def meq_linear_part_ff(F)
                    203: {
                    204:        F = simp_ff(F);
                    205:        F = F/LCOEF(F);
                    206:        V = var(F);
                    207:        if ( deg(F,V) == 1 )
                    208:                return [F,1,0];
                    209: T0 = time()[0];
                    210:        XP = pwrmod_ff(F);
                    211:        GCD = ugcd(F,XP-V);
                    212:        if ( deg(GCD,V) ) {
                    213:                GCD = GCD/LCOEF(GCD);
                    214:                F = sdiv(F,GCD);
                    215:                XP = srem(XP,F);
                    216:                R = GCD;
                    217:        } else
                    218:                R = 1;
                    219: TPMOD += time()[0]-T0;
                    220:        return [R,F,XP];
                    221: }
                    222:
                    223: /*
                    224:        Input : a square free polynomial F s.t.
                    225:                all the irreducible factors of F
                    226:                has the same degree D.
                    227:        Output: D
                    228: */
                    229:
                    230: def meq_ed_ff(F,XP)
                    231: {
                    232: T0 = time()[0];
                    233:        F = simp_ff(F);
                    234:        F = F/LCOEF(F);
                    235:        V = var(F);
                    236:
                    237:        TAB = pwrtab_ff(F,XP);
                    238:
                    239:        D = deg(F,V);
                    240:        for ( I = 1, W = V, L = []; 2*I <= D; I++ ) {
                    241:                lazy_lm(1);
                    242:                for ( T = 0, K = 0; K <= deg(W,V); K++ )
                    243:                        if ( C = coef(W,K,V) )
                    244:                                T += TAB[K]*C;
                    245:                lazy_lm(0);
                    246:                W = simp_ff(T);
                    247:                if ( W == V ) {
                    248:                        D = I; break;
                    249:                }
                    250:        }
                    251: TQMOD += time()[0]-T0;
                    252:        return D;
                    253: }
                    254:
                    255: /*
                    256:        Input : a square free polynomial F
                    257:                an integer E
                    258:             an array TAB
                    259:             where all the irreducible factors of F has degree E
                    260:             and TAB[i] = V^(i*Ord) mod F. (V=var(F), Ord=field order)
                    261:     Output: a list containing all the irreducible factors of F
                    262: */
                    263:
                    264: def berlekamp_ff(F,E,TAB)
                    265: {
                    266:        V = var(F);
                    267:        N = deg(F,V);
                    268:        Q = newmat(N,N);
                    269:        for ( J = 0; J < N; J++ ) {
                    270:                T = urem(TAB[J],F);
                    271:                for ( I = 0; I < N; I++ ) {
                    272:                        Q[I][J] = coef(T,I);
                    273:                }
                    274:        }
                    275:        for ( I = 0; I < N; I++ )
                    276:                Q[I][I] -= simp_ff(1);
                    277:        L = nullspace_ff(Q); MT = L[0]; IND = L[1];
                    278:        NF0 = N/E;
                    279:        PS = null_to_poly_ff(MT,IND,V);
                    280:        R = newvect(NF0); R[0] = F/LCOEF(F);
                    281:        for ( I = 1, NF = 1; NF < NF0 && I < NF0; I++ ) {
                    282:                PSI = PS[I];
                    283:                MP = minipoly_ff(PSI,F);
                    284:                ROOT = find_root_ff(MP); NR = length(ROOT);
                    285:                for ( J = 0; J < NF; J++ ) {
                    286:                        if ( deg(R[J],V) == E )
                    287:                                continue;
                    288:                        for ( K = 0; K < NR; K++ ) {
                    289:                                GCD = ugcd(R[J],PSI-ROOT[K]);
                    290:                                if ( deg(GCD,V) > 0 && deg(GCD,V) < deg(R[J],V) ) {
                    291:                                        Q = sdiv(R[J],GCD);
                    292:                                        R[J] = Q; R[NF++] = GCD;
                    293:                                }
                    294:                        }
                    295:                }
                    296:        }
                    297:        return vtol(R);
                    298: }
                    299:
                    300: /*
                    301:        Input : a matrix MT
                    302:                an array IND
                    303:                an indeterminate V
                    304:             MT is a matrix after Gaussian elimination
                    305:             IND[I] = 0 means that i-th column of MT represents a basis
                    306:             element of the null space.
                    307:        Output: an array R which contains all the basis element of
                    308:                 the null space of MT. Here, a basis element E is represented
                    309:                 as a polynomial P of V s.t. coef(P,i) = E[i].
                    310: */
                    311:
                    312: def null_to_poly_ff(MT,IND,V)
                    313: {
                    314:        N = size(MT)[0];
                    315:        for ( I = 0, J = 0; I < N; I++ )
                    316:                if ( IND[I] )
                    317:                        J++;
                    318:        R = newvect(J);
                    319:        for ( I = 0, L = 0; I < N; I++ ) {
                    320:                if ( !IND[I] )
                    321:                        continue;
                    322:                for ( J = K = 0, T = 0; J < N; J++ )
                    323:                        if ( !IND[J] )
                    324:                                T += MT[K++][I]*V^J;
                    325:                        else if ( J == I )
                    326:                                T -= V^I;
                    327:                R[L++] = simp_ff(T);
                    328:        }
                    329:        return R;
                    330: }
                    331:
                    332: /*
                    333:        Input : a polynomial P, a polynomial F
                    334:        Output: a minimal polynomial MP(V) of P mod F.
                    335: */
                    336:
                    337: def minipoly_ff(P,F)
                    338: {
                    339:        V = var(P);
                    340:        P0 = P1 = simp_ff(1);
                    341:        L = [[P0,P0]];
                    342:        while ( 1 ) {
                    343:                /* P0 = V^K, P1 = P^K mod F */
                    344:                P0 *= V;
                    345:                P1 = urem(P*P1,F);
                    346:                /*
                    347:                NP0 = P0-c1L1_0-c2L2_0-...,
                    348:             NP1 is a normal form w.r.t. [L1_1,L2_1,...]
                    349:                    NP1 = P1-c1L1_1-c2L2_1-...,
                    350:             NP0 represents the normal form computation.
                    351:             */
                    352:                L1 = lnf_ff(P0,P1,L,V); NP0 = L1[0]; NP1 = L1[1];
                    353:                if ( !NP1 )
                    354:                        return NP0;
                    355:                else
                    356:                        L = lnf_insert([NP0,NP1],L,V);
                    357:        }
                    358: }
                    359:
                    360: /*
                    361:        Input ; a list of polynomials [P0,P1] = [V^K,P^K mod F]
                    362:                a sorted list L=[[L1_0,L1_1],[L2_0,L2_1],...]
                    363:                of previously computed pairs of polynomials
                    364:                an indeterminate V
                    365:        Output: a list of polynomials [NP0,NP1]
                    366:                where NP1 = P1-c1L1_1-c2L2_1-...,
                    367:                      NP0 = P0-c1L1_0-c2L2_0-...,
                    368:             NP1 is a normal form w.r.t. [L1_1,L2_1,...]
                    369:             NP0 represents the normal form computation.
                    370:                [L1_1,L_2_1,...] is sorted so that it is a triangular
                    371:                linear basis s.t. deg(Li_1) > deg(Lj_1) for i<j.
                    372: */
                    373:
                    374: def lnf_ff(P0,P1,L,V)
                    375: {
                    376:        NP0 = P0; NP1 = P1;
                    377:        for ( T = L; T != []; T = cdr(T) ) {
                    378:                Q = car(T);
                    379:                D1 = deg(NP1,V);
                    380:                if ( D1 == deg(Q[1],V) ) {
                    381:                        H = -coef(NP1,D1,V)/coef(Q[1],D1,V);
                    382:                        NP0 += Q[0]*H;
                    383:                        NP1 += Q[1]*H;
                    384:                }
                    385:        }
                    386:        return [NP0,NP1];
                    387: }
                    388:
                    389: /*
                    390:        Input : a pair of polynomial P=[P0,P1],
                    391:                a list L,
                    392:                an indeterminate V
                    393:        Output: a list L1 s.t. L1 contains P and L
                    394:                and L1 is sorted in the decreasing order
                    395:                w.r.t. the degree of the second element
                    396:                of elements in L1.
                    397: */
                    398:
                    399: def lnf_insert(P,L,V)
                    400: {
                    401:        if ( L == [] )
                    402:                return [P];
                    403:        else {
                    404:                P0 = car(L);
                    405:                if ( deg(P0[1],V) > deg(P[1],V) )
                    406:                        return cons(P0,lnf_insert(P,cdr(L),V));
                    407:                else
                    408:                        return cons(P,L);
                    409:        }
                    410: }
                    411:
                    412: /*
                    413:        Input : a square free polynomial F s.t.
                    414:                all the irreducible factors of F
                    415:                has the degree E.
                    416:        Output: a list containing all the irreducible factors of F
                    417: */
                    418:
                    419: def c_z_ff(F,E)
                    420: {
                    421:        Type = field_type_ff();
                    422:        if ( Type == 1 || Type == 3 )
                    423:                return c_z_lm(F,E);
                    424:        else
                    425:                return c_z_gf2n(F,E);
                    426: }
                    427:
                    428: /*
                    429:        Input : a square free polynomial P s.t. P splits into linear factors
                    430:        Output: a list containing all the root of P
                    431: */
                    432:
                    433: def find_root_ff(P)
                    434: {
                    435:        V = var(P);
                    436:        L = c_z_ff(P,1);
                    437:        for ( T = L, U = []; T != []; T = cdr(T) ) {
                    438:                S = car(T)/LCOEF(car(T)); U = cons(-coef(S,0,V),U);
                    439:        }
                    440:        return U;
                    441: }
                    442:
                    443: /*
                    444:        Input : a square free polynomial F s.t.
                    445:                all the irreducible factors of F
                    446:                has the degree E.
                    447:        Output: an irreducible factor of F
                    448: */
                    449:
                    450: def c_z_one_ff(F,E)
                    451: {
                    452:        Type = field_type_ff();
                    453:        if ( Type == 1 || Type == 3 )
                    454:                return c_z_one_lm(F,E);
                    455:        else
                    456:                return c_z_one_gf2n(F,E);
                    457: }
                    458:
                    459: /*
                    460:        Input : a square free polynomial P s.t. P splits into linear factors
                    461:        Output: a list containing a root of P
                    462: */
                    463:
                    464: def find_one_root_ff(P)
                    465: {
                    466:        V = var(P);
                    467:        LF = c_z_one_ff(P,1);
                    468:        U = -coef(LF/LCOEF(LF),0,V);
                    469:        return [U];
                    470: }
                    471:
                    472: /*
                    473:        Input : an integer N; an indeterminate V
                    474:        Output: a polynomial F s.t. var(F) = V, deg(F) < N
                    475:                and its coefs are random numbers in
                    476:                the ground field.
                    477: */
                    478:
                    479: def randpoly_ff(N,V)
                    480: {
                    481:        for ( I = 0, S = 0; I < N; I++ )
                    482:                S += random_ff()*V^I;
                    483:        return S;
                    484: }
                    485:
                    486: /*
                    487:        Input : an integer N; an indeterminate V
                    488:        Output: a monic polynomial F s.t. var(F) = V, deg(F) = N-1
                    489:                and its coefs are random numbers in
                    490:                the ground field except for the leading term.
                    491: */
                    492:
                    493: def monic_randpoly_ff(N,V)
                    494: {
                    495:        for ( I = 0, N1 = N-1, S = 0; I < N1; I++ )
                    496:                S += random_ff()*V^I;
                    497:        return V^N1+S;
                    498: }
                    499:
                    500: /* GF(p) specific functions */
                    501:
                    502: /*
                    503:        Input : a square free polynomial F s.t.
                    504:                all the irreducible factors of F
                    505:                has the degree E.
                    506:        Output: a list containing all the irreducible factors of F
                    507: */
                    508:
                    509: def c_z_lm(F,E)
                    510: {
                    511:        V = var(F);
                    512:        N = deg(F,V);
                    513:        if ( N == E )
                    514:                return [F];
                    515:        M = field_order_ff();
                    516:        K = idiv(N,E);
                    517:        L = [F];
                    518:        while ( 1 ) {
                    519:                W = monic_randpoly_ff(2*E,V);
                    520:                T = generic_pwrmod_ff(W,F,idiv(M^E-1,2));
                    521:                W = T-1;
                    522:                if ( !W )
                    523:                        continue;
                    524:                G = ugcd(F,W);
                    525:                if ( deg(G,V) && deg(G,V) < N ) {
                    526:                        L1 = c_z_lm(G,E);
                    527:                        L2 = c_z_lm(sdiv(F,G),E);
                    528:                        return append(L1,L2);
                    529:                }
                    530:        }
                    531: }
                    532:
                    533: /*
                    534:        Input : a square free polynomial F s.t.
                    535:                all the irreducible factors of F
                    536:                has the degree E.
                    537:        Output: an irreducible factor of F
                    538: */
                    539:
                    540: def c_z_one_lm(F,E)
                    541: {
                    542:        V = var(F);
                    543:        N = deg(F,V);
                    544:        if ( N == E )
                    545:                return F;
                    546:        M = field_order_ff();
                    547:        K = idiv(N,E);
                    548:        while ( 1 ) {
                    549:                W = monic_randpoly_ff(2*E,V);
                    550:                T = generic_pwrmod_ff(W,F,idiv(M^E-1,2));
                    551:                W = T-1;
                    552:                if ( W ) {
                    553:                        G = ugcd(F,W);
                    554:                        D = deg(G,V);
                    555:                        if ( D && D < N ) {
                    556:                                if ( 2*D <= N ) {
                    557:                                        F1 = G; F2 = sdiv(F,G);
                    558:                                } else {
                    559:                                        F2 = G; F1 = sdiv(F,G);
                    560:                                }
                    561:                                return c_z_one_lm(F1,E);
                    562:                        }
                    563:                }
                    564:        }
                    565: }
                    566:
                    567: /* GF(2^n) specific functions */
                    568:
                    569: /*
                    570:        Input : a square free polynomial F s.t.
                    571:                all the irreducible factors of F
                    572:                has the degree E.
                    573:        Output: a list containing all the irreducible factors of F
                    574: */
                    575:
                    576: def c_z_gf2n(F,E)
                    577: {
                    578:        V = var(F);
                    579:        N = deg(F,V);
                    580:        if ( N == E )
                    581:                return [F];
                    582:        K = idiv(N,E);
                    583:        L = [F];
                    584:        while ( 1 ) {
                    585:                W = randpoly_ff(2*E,V);
                    586:                T = tracemod_gf2n(W,F,E);
                    587:                W = T-1;
                    588:                if ( !W )
                    589:                        continue;
                    590:                G = ugcd(F,W);
                    591:                if ( deg(G,V) && deg(G,V) < N ) {
                    592:                        L1 = c_z_gf2n(G,E);
                    593:                        L2 = c_z_gf2n(sdiv(F,G),E);
                    594:                        return append(L1,L2);
                    595:                }
                    596:        }
                    597: }
                    598:
                    599: /*
                    600:        Input : a square free polynomial F s.t.
                    601:                all the irreducible factors of F
                    602:                has the degree E.
                    603:        Output: an irreducible factor of F
                    604: */
                    605:
                    606: def c_z_one_gf2n(F,E)
                    607: {
                    608:        V = var(F);
                    609:        N = deg(F,V);
                    610:        if ( N == E )
                    611:                return F;
                    612:        K = idiv(N,E);
                    613:        while ( 1 ) {
                    614:                W = randpoly_ff(2*E,V);
                    615:                T = tracemod_gf2n(W,F,E);
                    616:                W = T-1;
                    617:                if ( W ) {
                    618:                        G = ugcd(F,W);
                    619:                        D = deg(G,V);
                    620:                        if ( D && D < N ) {
                    621:                                if ( 2*D <= N ) {
                    622:                                        F1 = G; F2 = sdiv(F,G);
                    623:                                } else {
                    624:                                        F2 = G; F1 = sdiv(F,G);
                    625:                                }
                    626:                                return c_z_one_gf2n(F1,E);
                    627:                        }
                    628:                }
                    629:        }
                    630: }
                    631:
                    632: /*
                    633:        Input : an integer D
                    634:        Output: an irreducible polynomial F over GF(2)
                    635:                of degree D.
                    636: */
                    637:
                    638: def defpoly_mod2(D)
                    639: {
                    640:        return gf2ntop(irredpoly_up2(D,0));
                    641: }
                    642:
                    643: def dummy_time() {
                    644:        return [0,0,0,0];
                    645: }
                    646: end$

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