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