[BACK]Return to ntl.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / ox_ntl

Annotation of OpenXM/src/ox_ntl/ntl.rr, Revision 1.2

1.2     ! iwane       1: /* $OpenXM: OpenXM/src/ox_ntl/ntl.rr,v 1.1 2003/11/03 03:11:21 iwane Exp $ */
1.1       iwane       2:
                      3: module ntl;
                      4: localf factor$
1.2     ! iwane       5: localf lll$
1.1       iwane       6: localf ex_data$
                      7: localf ex_data_tmp$
1.2     ! iwane       8: localf mat2list$
        !             9: localf list2mat$
1.1       iwane      10:
                     11:
                     12: /* static variables */
                     13:
                     14: /* extern variables */
                     15:
                     16: #if 1
                     17: localf test$
                     18:
                     19: /*&usage begin: ntl.test(PID, POLY)
                     20:  compare on ox_NTL and on asir.
                     21:
                     22: example:
                     23: [1028] F=ntl.ex_data(4);
                     24: x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225
                     25: [1029] F = F * subst(F, x, x + 1)$
                     26: [1030] ntl.factor(PID, F);
                     27: [[1,1],[x^16+16*x^15-16*x^14-1344*x^13-4080*x^12+32576*x^11+157376*x^10-255232*x^9-2062624*x^8-249088*x^7+10702080*x^6+9126912*x^5-18643712*x^4-24167424*x^3+2712576*x^2+10653696*x+2324736,1],[x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225,1]]
                     28: [1031] ntl.test(PID, F);
                     29: [CPU,0.121539,0.001354,GC,0.0222,0]
                     30:
                     31:
                     32: end: */
                     33: def test(PID, F)
                     34: {
                     35:        T0 = time();
                     36:        fctr(F);
                     37:        T1 = time();
                     38:        ntl.factor(PID, F);
                     39:        T2 = time();
                     40:
                     41:        return (["CPU", T1[0]-T0[0],T2[0]-T1[0],"GC",T1[1]-T0[1],T2[1]-T1[1]]);
                     42: }
                     43:
                     44: #endif
                     45:
                     46: /*&usage begin: ntl.factor(PID, POLY)
                     47:
                     48:   Factorize polynomial {POLY} over the rationals.
                     49:
                     50: {RETURN}
                     51:        list
                     52: {POLY}
                     53:        univariate polynomial with rational coefficients
                     54:
                     55: description:
                     56:   Factorizes polynomial {POLY} over the rationals.
                     57:   The result is represented by a list, whose elements are a pair represented as
                     58:
                     59: [[num,1],[factor1,multiplicity1],[factor2,multiplicity2],...].
                     60:
                     61: Products of all factor^multiplicity and num is equal to {POLY}.
                     62:
                     63: The number {num} is determined so that ({POLY}/{num}) is an integral polynomial
                     64: and its content (GCD of all coefficients) is 1.
                     65:
                     66: author: H.Iwane <iwane@math.sci.kobe-u.ac.jp>
                     67:
                     68: example:
                     69: [1282] F=(x^5-1)*(x^3+1)^2;
                     70: x^11+2*x^8-x^6+x^5-2*x^3-1
                     71: [1283] ntl.factor(PID, F);
                     72: [[1,1],[x^4+x^3+x^2+x+1,1],[x-1,1],[x+1,2],[x^2-x+1,2]]
                     73: [1284] fctr(F);
                     74: [[1,1],[x-1,1],[x^4+x^3+x^2+x+1,1],[x+1,2],[x^2-x+1,2]]
                     75:
                     76: ref: fctr
                     77: end:
                     78: */
                     79: def factor(PID, POLY)
                     80: {
                     81:        local F, C, LIST, STR, RET, RET_NTL, VAR, I;
                     82:        local TYPE;
                     83:
                     84:        /* 入力チェック */
                     85:        TYPE = type(POLY);
                     86:        if (TYPE == 0 || TYPE == 1) {
                     87:                return ([[POLY,1]]);
                     88:        } else if (TYPE != 2) {
                     89:                error("ntl.factor: invalid argument");
                     90:        }
                     91:
                     92:
                     93:        LIST = vars(POLY);
                     94:        if (length(LIST) >= 2) { /* 一変数多項式のみ */
                     95:                error("ntl.factor: invalid argument");
                     96:        }
                     97:
                     98:
                     99:        /* NTL で 有理係数多項式は不可 */
                    100:        F = ptozp(POLY);
                    101:
                    102:        C = sdiv(POLY, F);
                    103:
                    104:        ox_cmo_rpc(PID, "fctr", F);
                    105:
                    106:        RET_NTL = ox_pop_cmo(PID);
                    107:
                    108:        /* ERROR Check */
                    109:        if (type(RET_NTL) != 4 || length(RET_NTL) < 2) {
1.2     ! iwane     110:                return (RET_NTL);
1.1       iwane     111:        }
                    112:
                    113:        RET = cons([RET_NTL[0] * C, 1], RET_NTL[1]);
                    114:
                    115:        return (RET);
                    116: }
                    117:
                    118: /*&usage begin: ex_data_tmp(F, N)
                    119:   Generate sample irreducible polynomial
                    120:
                    121: example:
                    122: [1032] F=t^3-p;
                    123: -p+t^3
                    124: [1033] ntl.ex_data_tmp(F,3);
                    125: x^27-90*x^24+1089*x^21-62130*x^18+105507*x^15-16537410*x^12-30081453*x^9-1886601330*x^6+73062900*x^3-6859000
                    126: [1034] ntl.ex_data_tmp(F,4);
                    127: -x^81+459*x^78-76896*x^75+7538094*x^72-347721147*x^69+3240161703*x^66+1032617170332*x^63-37499673798288*x^60+784360767442050*x^57+150576308695750650*x^54+771023617441694964*x^51+67248913649472410184*x^48+13913995714637027898294*x^45+270221527865987051874714*x^42-8828542741395296724347412*x^39-154971101776040822743879716*x^36+4343529580943017469231383983*x^33+4648027555241630173815780123*x^30-1072436585643253024332438894564*x^27+16237394255218510503554781142602*x^24-134104542851048701593527668875195*x^21+727430949790032393675174790142991*x^18-2727255031466780569416130788693624*x^15+7102683996190585423335589883738868*x^12-12524463688445776069953452180105904*x^9+14119870779369458313271232460210576*x^6-8591747198480108022372636571451136*x^3+2346749360904699138972190279765184
                    128:
                    129: ref: ex_data
                    130:
                    131: end: */
                    132: def ex_data_tmp(F, N)
                    133: {
                    134:        FP = subst(F, p, prime(0));
                    135:        GP = subst(FP, t, x);
                    136:
                    137:        for (I = 1; I < N; I++) {
                    138:                FP = subst(F, p, prime(I));
                    139:                GP = res(t, subst(GP, x, x-t), FP);
                    140:        }
                    141:
                    142:        return (GP);
                    143: }
                    144:
                    145: /*&usage begin: ntl.ex_data(N)
                    146:   Generate sample irreducible polynomial
                    147:
                    148: example:
                    149: [1028] ntl.ex_data(3);
                    150: x^8-40*x^6+352*x^4-960*x^2+576
                    151: [1029] ntl.ex_data(4);
                    152: x^16-136*x^14+6476*x^12-141912*x^10+1513334*x^8-7453176*x^6+13950764*x^4-5596840*x^2+46225
                    153:
                    154: ref:
                    155: ex_data_tmp
                    156:
                    157: end: */
                    158: def ex_data(N)
                    159: {
                    160:        if (type(N) != 1 && type(N) != 10) {
                    161:                print("invalid argument");
                    162:                return (0);
                    163:        }
                    164:
                    165:        return (ex_data_tmp(t^2-p, N));
1.2     ! iwane     166: }
        !           167:
        !           168:
        !           169: def mat2list(M)
        !           170: {
        !           171:        A = size(M);
        !           172:
        !           173:        ROW=A[0];
        !           174:        COL=A[1];
        !           175:
        !           176:        for (I = 0; I < ROW; I++) {
        !           177:                for (J = 0; J < COL; J++) {
        !           178:                        A = append(A, [M[I][J]]);
        !           179:                }
        !           180:        }
        !           181:
        !           182:        return (A);
        !           183: }
        !           184:
        !           185: def list2mat(L)
        !           186: {
        !           187:        if (type(L) != 4) {
        !           188:                return ("Invalid Argument");
        !           189:        }
        !           190:
        !           191:        ROW = L[0];
        !           192:        if (type(ROW) == 10)
        !           193:                ROW = int32ton(ROW);
        !           194:        COL = L[1];
        !           195:        if (type(COL) == 10)
        !           196:                COL = int32ton(COL);
        !           197:
        !           198:        A = newmat(2, 2); /*, [[1, 0],[0, 1]]); /* COL, COL); */
        !           199:
        !           200:        C = 2;
        !           201:        for (I = 0; I < ROW; I++) {
        !           202:                for (J = 0; J < COL; J++) {
        !           203:                        A[I][J] = L[C];
        !           204:                        C++;
        !           205:                }
        !           206:        }
        !           207:
        !           208:
        !           209:        return (A);
        !           210: }
        !           211:
        !           212:
        !           213:
        !           214: /*&usage begin: ntl.lll(PID, MAT)
        !           215:   the basics of LLL reducation.
        !           216:
        !           217: {M}
        !           218:        Matrix which element is Integer.
        !           219:
        !           220: example:
        !           221: [1081] M=newmat(2,2,[[10,0],[-7,3]]);
        !           222: [ 10 0 ]
        !           223: [ -7 3 ]
        !           224: [1082] ntl.lll(PID, M);
        !           225: [ 3 3 ]
        !           226: [ 4 -6 ]
        !           227: [1083] pari(lll, M);
        !           228: [ 0 1 ]
        !           229: [ 1 2 ]                 <== why ?
        !           230:
        !           231:
        !           232: ref:
        !           233: pari(lll)
        !           234:
        !           235: end: */
        !           236: def lll(PID, M)
        !           237: {
        !           238:        /* parameter check */
        !           239:        TYPE = type(M);
        !           240:        if (TYPE != 6) {        /* matrix */
        !           241:                error("ntl.lll: invalid argument");
        !           242:        }
        !           243:
        !           244:        A = mat2list(M);
        !           245:
        !           246:        ox_cmo_rpc(PID, "lll", A);
        !           247:
        !           248:        RET_NTL = ox_pop_cmo(PID);
        !           249:
        !           250:        /* return value check */
        !           251:        if (type(RET_NTL) != 4) { /* list */
        !           252:                error("ntl.lll: error");
        !           253:        }
        !           254:
        !           255:        R = list2mat(RET_NTL);
        !           256:
        !           257:        return (R);
1.1       iwane     258: }
                    259:
                    260:
                    261: endmodule;
                    262:
                    263:
                    264: end$

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