![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to nf CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to nf CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.2 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/nf,v 1.1.1.1 1999/12/03 07:39:11 noro Exp $
! 49: */
1.1 noro 50: extern IDIVV_REM,IDIVV_HIST;
51:
52: def nf_mod(B,G,MOD,DIR,PS)
53: {
54: setmod(MOD);
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: R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
59: CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
60: U = G-CR*R;
61: print(".",2);
62: if ( !U )
63: return D;
64: break;
65: }
66: }
67: if ( U )
68: G = U;
69: else {
70: D += dp_hm(G); G = dp_rest(G); print(!D,2);
71: }
72: }
73: return D;
74: }
75:
76: def nf_mod_ext(B,G,MOD,DIR,PS)
77: {
78: setmod(MOD);
79: for ( D = 0, H = []; G; ) {
80: for ( U = 0, L = B; L != []; L = cdr(L) ) {
81: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
82: R = dp_mod(bload(DIR+rtostr(car(L))),MOD,[]);
83: CR = inv(dp_hc(R),MOD)*dp_hc(G)*dp_mod(dp_subd(G,R),MOD,[]);
84: U = G-CR*R;
85: H = cons([CR,strtov("t"+rtostr(car(L)))],H);
86: print([length(H)]);
87: if ( !U )
88: return [D,reverse(H)];
89: break;
90: }
91: }
92: if ( U )
93: G = U;
94: else {
95: D += dp_hm(G); G = dp_rest(G);
96: }
97: }
98: return [D,reverse(H)];
99: }
100:
101: def nf(B,G,M,PS)
102: {
103: for ( D = 0, H = []; G; ) {
104: for ( U = 0, L = B; L != []; L = cdr(L) ) {
105: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
106: GCD = igcd(dp_hc(G),dp_hc(R));
107: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
108: U = CG*G-dp_subd(G,R)*CR*R;
109: H = cons(car(L),H);
110: if ( !U )
111: return [D,M,reverse(H)];
112: D *= CG; M *= CG;
113: break;
114: }
115: }
116: if ( U )
117: G = U;
118: else {
119: D += dp_hm(G); G = dp_rest(G);
120: }
121: }
122: return [D,M,reverse(H)];
123: }
124: extern M_NM,M_DN;
125:
126: def nf_demand(B,G,DIR,PSH)
127: {
128: T1 = T2 = T3 = T4 = 0;
129: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
130: for ( D = 0, H = []; G; ) {
131: for ( U = 0, L = B; L != []; L = cdr(L) ) {
132: if ( dp_redble(G,PSH[car(L)]) > 0 ) {
133: R = bload(DIR+rtostr(car(L)));
134: H = cons(car(L),H);
135: print([length(H),dp_mag(dp_hm(G))]);
136: GCD = igcd(dp_hc(G),dp_hc(R));
137: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
138: T0 = time()[0]; A1 = CG*G; T1 += time()[0]-T0;
139: T0 = time()[0]; A2 = dp_subd(G,R)*CR*R; T2 += time()[0]-T0;
140: T0 = time()[0]; U = A1-A2; T3 += time()[0]-T0;
141: if ( !U )
142: return [D,reverse(H),T1,T2,T3,T4];
143: D *= CG;
144: break;
145: }
146: }
147: if ( U ) {
148: G = U;
149: if ( D ) {
150: if ( dp_mag(dp_hm(D)) > Mag ) {
151: T0 = time()[0];
152: print("ptozp2");
153: T = dp_ptozp2(D,G); D = T[0]; G = T[1];
154: T4 += time()[0]-T0;
155: Mag = idiv(dp_mag(dp_hm(D))*M_NM,M_DN);
156: }
157: } else {
158: if ( dp_mag(dp_hm(G)) > Mag ) {
159: T0 = time()[0];
160: print("ptozp");
161: G = dp_ptozp(G);
162: T4 += time()[0]-T0;
163: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
164: }
165: }
166: } else {
167: D += dp_hm(G); G = dp_rest(G);
168: }
169: }
170: return [D,reverse(H),T1,T2,T3,T4];
171: }
172:
173: def nf_demand_d(B,G,DIR,NDIR,PDIR,PSH,PROC)
174: {
175: INDEX = 0;
176: T1 = T2 = 0;
177: /* dp_gr_flags(["Dist",PROC]); */
178: if ( PROC != [] ) {
179: PROC = newvect(length(PROC),PROC);
180: NPROC = size(PROC)[0];
181: } else
182: NPROC = 0;
183: Mag = idiv(dp_mag(dp_hm(G))*M_NM,M_DN);
184: Kara = dp_set_kara()*27; /* XXX */
185: D = [0,0]; R = [1,G];
186: print("");
187: for ( H = []; R[1]; ) {
188: for ( U = 0, L = B; L != []; L = cdr(L) ) {
189: if ( dp_redble(R[1],PSH[car(L)]) > 0 ) {
190: Red = bload(DIR+rtostr(car(L)));
191: H = cons(car(L),H);
192: /* D0=dp_mag(D[0]*<<0>>); D1=dp_mag(dp_hm(D[1]));
193: R0=dp_mag(R[0]*<<0>>); R1=dp_mag(dp_hm(R[1]));
194: print([car(L),length(H),[D0,D1,dp_nt(D[1])],[R0,R1,dp_nt(R[1])]],2); */
195:
196: GCD = igcd(dp_hc(R[1]),dp_hc(Red));
197: CR = idiv(dp_hc(Red),GCD);
198: CRed = idiv(dp_hc(R[1]),GCD);
199:
200: T0 = time()[3];
201: if ( (PROC != []) && (dp_mag(dp_hm(Red)) > Kara) ) {
202: print("d",2);
203: rpc(PROC[0],"dp_imul_index",CRed,car(L));
204: U = CR*R[1] - dp_subd(R[1],Red)*rpcrecv(PROC[0]);
205: } else {
206: print("l",2);
207: U = CR*R[1] - dp_subd(R[1],Red)*CRed*Red;
208: }
209: TT = time()[3]-T0; T1 += TT; /* print(TT); */
210: if ( !U )
211: return [D,reverse(H),T1,T2];
212: break;
213: }
214: }
215: if ( U ) {
216: if ( (HMAG=dp_mag(dp_hm(U))) > Mag ) {
217: T0 = time()[3];
218: if ( HMAG > Kara ) {
219: print("D",2);
220: T = dp_ptozp_d(U,PROC,NPROC);
221: } else {
222: print("L",2);
223: T = dp_ptozp(U);
224: }
225: TT = time()[3]-T0; T2 += TT; /* print(TT); */
226: Cont = idiv(dp_hc(U),dp_hc(T));
227: D0 = kmul(CR,D[0]);
228: R0 = kmul(Cont,R[0]);
229: S = igcd(D0,R0);
230: D = [idiv(D0,S),D[1]];
231: R = [idiv(R0,S),T];
232: Mag = idiv(dp_mag(dp_hm(R[1]))*M_NM,M_DN);
233: } else {
234: D = [kmul(CR,D[0]),D[1]];
235: R = [R[0],U];
236: }
237: } else {
238: C = kmul(dp_hc(R[1]),R[0]);
239: T = igcd(D[0],C);
240: D = [T,idiv(D[0],T)*D[1]+idiv(C,T)*dp_ht(R[1])];
241: R = [R[0],dp_rest(R[1])];
242: }
243: }
244: return [D[1],reverse(H),T1,T2];
245: }
246:
247: extern PTOZP,DPCV,NFSDIR;
248:
249: def imulv(A,B)
250: {
251: return A*B;
252: }
253:
254: def dp_imul(P,N)
255: {
256: return N*P;
257: }
258:
259: def mul(A,B)
260: {
261: return A*B;
262: }
263:
264: def dp_imul_index(C,INDEX)
265: {
266: T = C*bload(NFSDIR+rtostr(INDEX));
267: return T;
268: }
269:
270: def reg_dp_dtov(P)
271: {
272: PTOZP = P;
273: DPCV = dp_dtov(PTOZP);
274: }
275:
276: def reg_dp_iqr(G)
277: {
278: N = size(DPCV)[0];
279: for ( R = [], I = 0; I < N; I++ ) {
280: DPCV[I] = T = irem(DPCV[I],G);
281: if ( T )
282: R = cons(T,R);
283: }
284: return R;
285: }
286:
287: def reg_dp_cont(P)
288: {
289: PTOZP = P;
290: C = dp_dtov(P);
291: return igcd(C);
292: }
293:
294: def reg_dp_idiv(GCD)
295: {
296: return dp_idiv(PTOZP,GCD);
297: }
298:
299: def dp_cont_d(G,PROC,NPROC)
300: {
301: C = dp_sep(G,NPROC);
302: N = size(C)[0];
303: for ( I = 0; I < N; I++ )
304: rpc(PROC[I],"reg_dp_cont",C[I]);
305: R = map(rpcrecv,PROC);
306: GCD = igcd(R);
307: return GCD;
308: }
309:
310: /*
311: def iqrv(C,D)
312: {
313: return map(iqr,C,D);
314: }
315: */
316:
317: #if 1
318: def dp_ptozp_d(G,PROC,NPROC)
319: {
320: C = dp_dtov(G); L = size(C)[0];
321: T0 = time()[3];
322: D0 = igcd_estimate(C);
323: TT = time()[3]-T0; TG1 += TT;
324: CS = sepvect(C,NPROC+1);
325: N = size(CS)[0]; N1 = N-1;
326: QR = newvect(N); Q = newvect(L); R = newvect(L);
327: T0 = time()[3];
328: for ( I = 0; I < N1; I++ )
329: rpc0(PROC[I],"iqrv",CS[I],D0);
330: QR[N1] = iqrv(CS[N1],D0);
331: for ( I = 0; I < N1; I++ )
332: QR[I] = rpcrecv(PROC[I]);
333: TT = time()[3]-T0; TD += TT;
334: for ( I = J = 0; I < N; I++ ) {
335: T = QR[I]; M = size(T)[0];
336: for ( K = 0; K < M; K++, J++ ) {
337: Q[J] = T[K][0]; R[J] = T[K][1];
338: }
339: }
340: T0 = time()[3];
341: D1 = igcd(R);
342: GCD = igcd(D0,D1);
343: A = idiv(D0,GCD);
344: TT = time()[3]-T0; TG2 += TT;
345: T0 = time()[3];
346: for ( I = 0; I < L; I++ )
347: Q[I] = kmul(A,Q[I])+idiv(R[I],GCD);
348: TT = time()[3]-T0; TM += TT;
349: print([TG1,TD,TG2,TM,dp_mag(D0*<<0>>),dp_mag(D1*<<0>>),dp_mag(GCD*<<0>>)],2);
350: return dp_vtod(Q,G);
351: }
352: #endif
353:
354: #if 0
355: def dp_ptozp_d(G,PROC,NPROC)
356: {
357: C = dp_sep(G,NPROC+1);
358: N = size(C)[0]; N1 = N-1;
359: R = newvect(N);
360: for ( I = 0; I < N1; I++ )
361: rpc(PROC[I],"reg_igcdv",dp_dtov(C[I]));
362: R[N1] = reg_igcdv(dp_dtov(C[N1]));
363: for ( I = 0; I < N1; I++ )
364: R[I] = rpcrecv(PROC[I]);
365: GCD = igcd(R);
366: for ( I = 0; I < N1; I++ )
367: rpc(PROC[I],"idivv_final",GCD);
368: S = dp_vtod(idivv_final(GCD),C[N1]);
369: for ( I = 0; I < N1; I++ )
370: S += dp_vtod(rpcrecv(PROC[I]),C[I]);
371: return S;
372: }
373: #endif
374:
375: #if 0
376: def dp_ptozp_d(G,PROC,NPROC)
377: {
378: C = dp_sep(G,NPROC+1);
379: N = size(C)[0]; N1 = N-1;
380: R = newvect(N);
381: for ( I = 0; I < N1; I++ )
382: rpc(PROC[I],"reg_dp_cont",C[I]);
383: R[N1] = igcd(dp_dtov(C[N1]));
384: for ( I = 0; I < N1; I++ )
385: R[I] = rpcrecv(PROC[I]);
386: GCD = igcd(R);
387: for ( I = 0; I < N1; I++ )
388: rpc(PROC[I],"reg_dp_idiv",GCD);
389: S = dp_idiv(C[N1],GCD);
390: for ( I = 0; I < N1; I++ )
391: S += rpcrecv(PROC[I]);
392: return S;
393: }
394: #endif
395:
396: #if 0
397: def dp_ptozp_d(G,PROC,NPROC)
398: {
399: C = dp_sep(G,NPROC);
400: N = size(C)[0]; T = newvect(N);
401: for ( I = 0; I < N; I++ )
402: T[I] = PROC[I];
403: PROC = T;
404: for ( I = 0; I < N; I++ )
405: rpc(PROC[I],"reg_dp_dtov",C[I]);
406: A = dp_dtov(G); A = isort(A); L = size(A)[0];
407: map(rpcrecv,PROC);
408: while ( 1 ) {
409: GCD = igcd(A[0],A[1]);
410: for ( I = 2; I < L; I++ ) {
411: GT = igcd(GCD,A[I]);
412: if ( GCD == GT )
413: break;
414: else
415: GCD = GT;
416: }
417: if ( I == L )
418: break;
419: else {
420: map(rpc,PROC,"reg_dp_iqr",GCD);
421: R = map(rpcrecv,PROC);
422: for ( J = 0, U = []; J < N; J++ )
423: U = append(R[J],U);
424: L = length(U);
425: if ( L == 0 )
426: break;
427: else
428: U = cons(GCD,U);
429: A = isort(newvect(L+1,U));
430: print(["length=",L,"minmag=",dp_mag(A[0]*<<0>>)]);
431: }
432: }
433: for ( I = 0; I < N; I++ )
434: rpc(PROC[I],"reg_dp_idiv",GCD);
435: R = map(rpcrecv,PROC);
436: for ( I = 0, S = 0; I < N; I++ )
437: S += R[I];
438: return S;
439: }
440: #endif
441:
442: def dp_ptozp2_d(D,G,PROC,NPROC)
443: {
444: T = dp_ptozp_d(D+G,PROC,NPROC);
445: for ( S = 0; D; D = dp_rest(D), T = dp_rest(T) )
446: S += dp_hm(T);
447: return [S,T];
448: }
449:
450: def genindex(N)
451: {
452: R = [];
453: for ( I = 0; I < N; I++ )
454: R = cons(I,R);
455: return reverse(R);
456: }
457:
458: def nftest_ext_mod(N1,N2,N,MOD,DIR,PSH)
459: {
460: S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
461: R = nf_mod_ext(genindex(N-1),S,MOD,DIR,PSH);
462: return R;
463: }
464:
465: def nftest_mod(N1,N2,N,MOD,DIR,PSH)
466: {
467: S = dp_mod(dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2)))),MOD,[]);
468: R = nf_mod(genindex(N-1),S,MOD,DIR,PSH);
469: return dp_mdtod(R);
470: }
471:
472: def nftest(N1,N2,N,DIR,PSH)
473: {
474: S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
475: R = nf_demand(genindex(N-1),S,DIR,PSH);
476: return R;
477: }
478:
479: def nftest_d(N1,N2,N,DIR,NDIR,PDIR,PSH,PROC)
480: {
481: S = dp_ptozp(dp_sp(bload(DIR+rtostr(N1)),bload(DIR+rtostr(N2))));
482: R = nf_demand_d(genindex(N-1),S,DIR,NDIR,PDIR,PSH,PROC);
483: return R;
484: }
485:
486: def abs(A)
487: {
488: return A >= 0 ? A : -A;
489: }
490:
491: def sort(A)
492: {
493: N = size(A)[0];
494: for ( I = 0; I < N; I++ ) {
495: for ( M = I, J = I + 1; J < N; J++ )
496: if ( A[J] < A[M] )
497: M = J;
498: if ( I != M ) {
499: T = A[I]; A[I] = A[M]; A[M] = T;
500: }
501: }
502: return A;
503: }
504:
505: def igcdv(C)
506: {
507: C = sort(C); N = size(C)[0];
508: G = igcd(C[0],C[1]);
509: for ( I = 2; I < N; I++ ) {
510: T = time();
511: G = igcd(G,C[I]);
512: /* print([I,dp_mag(G*<<0>>),time()[0]-T[0]]); */
513: }
514: return G;
515: }
516:
517: def igcdv2(C)
518: {
519: C = sort(C); N = size(C)[0];
520: G = igcd(C[0],C[1]);
521: for ( I = 2; I < N; I++ ) {
522: T = time();
523: C[I] %= G;
524: GT = igcd(G,C[I]);
525: if ( G == GT ) {
526: for ( J = I+1, NZ=0; J < N; J++ ) {
527: C[J] %= G;
528: if ( C[J] )
529: NZ = 1;
530: }
531: if ( !NZ )
532: break;
533: } else
534: G = GT;
535: print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
536: }
537: return G;
538: }
539:
540: def igcdv_d(C,P,NP)
541: {
542: C = isort(C); N = size(C)[0];
543: G = igcd(C[0],C[1]);
544: for ( I = 2; I < N; I++ ) {
545: T = time();
546: C[I] %= G;
547: GT = igcd(G,C[I]);
548: if ( G == GT ) {
549: for ( J = I+1, NZ=0; J < N; J++ ) {
550: C[J] %= G;
551: if ( C[J] )
552: NZ = 1;
553: }
554: if ( !NZ )
555: break;
556: } else
557: G = GT;
558: print([I,dp_mag(G*<<0>>),time()[0]-T[0]]);
559: }
560: return G;
561: }
562:
563: def dup(A)
564: {
565: N = size(A)[0]; V = newvect(N);
566: for ( I = 0; I < N; I++ )
567: V[I] = A[I];
568: return V;
569: }
570:
571: def idivv_init(A)
572: {
573: IDIVV_REM = dup(A);
574: }
575:
576: def igcd_cofactor(D)
577: {
578: T = map(iqr,IDIVV_REM,D);
579: N = size(T)[0];
580: Q = newvect(N); R = newvect(N);
581: for ( I = 0; I < N; I++ ) {
582: Q[I] = T[I][0]; R[I] = T[I][1];
583: }
584: D1 = igcdv(dup(R));
585: if ( !D1 ) {
586: IDIVV_HIST = [[Q,D]];
587: return D;
588: } else {
589: Q1 = map(idiv,R,D1);
590: IDIVV_HIST = [[Q,D],[Q1,D1]];
591: return D1;
592: }
593: }
594:
595: def idivv_step(D)
596: {
597: T = map(iqr,IDIVV_REM,D);
598: N = size(T)[0];
599: Q = newvect(N); R = newvect(N);
600: for ( I = 0, NZ = 0; I < N; I++ ) {
601: Q[I] = T[I][0]; R[I] = T[I][1];
602: if ( R[I] )
603: NZ = 1;
604: }
605: IDIVV_REM = R;
606: IDIVV_HIST = cons([Q,D],IDIVV_HIST);
607: return NZ;
608: }
609:
610: def igcdv3(C)
611: {
612: idivv_init(C);
613: C = isort(C); N = size(C)[0];
614: D = igcd(C[0],C[1]);
615: for ( I = 2; I < N; I++ ) {
616: T = time();
617: C[I] %= G;
618: GT = igcd(G,C[I]);
619: if ( G == GT ) {
620: NZ = idivv_step(G);
621: if ( !NZ )
622: break;
623: } else
624: G = GT;
625: }
626: CF = idivv_final(G);
627: return [G,CF];
628: }
629:
630: def igcdv4(C)
631: {
632: idivv_init(C);
633: C = isort(C); N = size(C)[0];
634: D = igcd_estimate(C);
635: G = igcd_cofactor(D);
636: G = igcd(G,D);
637: CF = idivv_final(G);
638: return [G,CF];
639: }
640:
641: def igcd_estimate(C)
642: {
643: C = isort(C); N = size(C)[0];
644: M = idiv(N,2);
645: for ( I = 0, S = 0; I < M; I++ )
646: S += C[I]>0?C[I]:-C[I];
647: for ( T = 0; I < N; I++ )
648: T += C[I]>0?C[I]:-C[I];
649: if ( !S && !T )
650: return igcdv(C);
651: else
652: return igcd(S,T);
653: }
654:
655: def reg_igcdv(C)
656: {
657: D = igcd_estimate(C);
658: T = map(iqr,C,D); N = size(T)[0];
659: Q = newvect(N); R = newvect(N);
660: for ( I = 0; I < N; I++ ) {
661: Q[I] = T[I][0]; R[I] = T[I][1];
662: }
663: D1 = igcdv(dup(R));
664: if ( !D1 ) {
665: IDIVV_HIST = [[Q,D]];
666: return D;
667: } else {
668: Q1 = map(idiv,R,D1);
669: IDIVV_HIST = [[Q,D],[Q1,D1]];
670: return igcd(D,D1);
671: }
672: }
673:
674: def idivv_final(D)
675: {
676: for ( T = IDIVV_HIST, S = 0; T != []; T = cdr(T) ) {
677: H = car(T); S += map(kmul,H[0],idiv(H[1],D));
678: }
679: return S;
680: }
681:
682: def dp_vtod(C,P)
683: {
684: for ( R = 0, I = 0; P; P = dp_rest(P), I++ )
685: R += C[I]*dp_ht(P);
686: return R;
687: }
688:
689: def dp_nt(P)
690: {
691: for ( I = 0; P; P = dp_rest(P), I++ );
692: return I;
693: }
694: end$