![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to const 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 const CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/lib/const,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
2: def cat(D) {
3: tstart;
4: for ( S = T = P = idiv(10^D,2), I = 1, J = 3; T; I++, J += 2 ) {
5: P = idiv(P*I,J); T = idiv(T*I+P,J); S += T;
6: }
7: tstop;
8: return S;
9: }
10:
11: def e(D,N)
12: {
13: for ( F = 1, S = 1, I = 1; I <= N; I++ ) {
14: S = S*I + 1;
15: F *= I;
16: }
17: T = red(S/F);
18: return idiv(nm(T)*10^D,dn(T));
19: }
20:
21: def at0(X,D)
22: {
23: for ( S = T = idiv(D,X), I = 1, Y = X^2, Sgn = -1;
24: T;
25: I += 2, Sgn *= -1 ) {
26: T = idiv(T*I,Y*(I+2)); S += (Sgn*T);
27: }
28: return S;
29: }
30:
31: def pi(D)
32: {
33: tstart; Y = 10^D; X = 16*at0(5,Y)-4*at0(239,Y); tstop;
34: return X;
35: }
36:
37: def at1(M,D) {
38: for (N = 1, SGN = 1, MM = M*M, A = 0, XN = idiv(D,M);
39: XN;
40: N += 2, XN = idiv(XN,MM), SGN *= -1)
41: A += (SGN*idiv(XN,N));
42: return A;
43: }
44:
45: def pi1(D) {
46: tstart; Y = 10^D; X = 16*at1(5,Y)-4*at1(239,Y); tstop;
47: return X;
48: }
49:
50: def pi2(D) {
51: tstart; Y = 10^D;
52: X = 48*at1(49,Y)+128*at1(57,Y)-20*at1(239,Y)+48*at1(110443,Y);
53: tstop;
54: return X;
55: }
56:
57:
58: def bn(N)
59: {
60: B = newvect(N+1); C = c2(N+1);
61: for ( I = 1, B[0] = 1; I <= N; I++ ) {
62: for ( D = C[I+1], J = 0, S = 0; J < I; J++ )
63: S += D[J]*B[J];
64: B[I] = red(-S/(I+1));
65: }
66: return [B,C];
67: }
68:
69: def bp(N,B,C,V)
70: {
71: for ( I = 0, S = 0; I <= N; I++ )
72: S += C[I]*B[N-I]*V^I;
73: return S;
74: }
75:
76: /*
77: * sum(N) = 1^N+2^N+...+n^N
78: */
79:
80: def sum(N)
81: {
82: L = bn(N+1);
83: R = car(L); C = car(cdr(L));
84: S = bp(N+1,R,C[N+1],n);
85: return red((subst(S,n,n+1)-subst(S,n,1))/(N+1));
86: }
87:
88: def c(N,I)
89: {
90: for ( M = 1, J = 0; J < I; J++ )
91: M *= N-J;
92: return red(M/f(I));
93: }
94:
95: def c1(N)
96: {
97: A = newvect(N+1); B = newvect(N+1); A[0] = 1;
98: for ( K = 1; K <= N; K++ ) {
99: B[0] = B[K] = 1;
100: for ( J = 1; J < K; J++ ) B[J] = A[J-1]+A[J];
101: T = A; A = B; B = T;
102: }
103: return A;
104: }
105:
106: def c2(N)
107: {
108: A = newvect(N+1); A[0] = B = newvect(1); B[0] = 1;
109: for ( K = 1; K <= N; K++ ) {
110: A[K] = B = newvect(K+1); B[0] = B[K] = 1;
111: for ( P = A[K-1], J = 1; J < K; J++ )
112: B[J] = P[J-1]+P[J];
113: }
114: return A;
115: }
116:
117: def f(N)
118: {
119: for ( I = 1, M = 1; I <= N; I++ )
120: M *= I;
121: return M;
122: }
123:
124: def sumd(N,M)
125: {
126: for ( I = 1, S = 0; I <= M; I++ )
127: S += I^N;
128: return S;
129: }
130:
131: #if 0
132: def sqrt(A,N) {
133: for ( I = 0, X = 1, B = A; I < N; I++, B *= 100, X *= 10 ) {
134: while ( 1 ) {
135: T = idiv(idiv(B,X) + X,2);
136: /*
137: if ((Y = T - X)== 0)
138: if ( B == X^2) return (X/(10^I));
139: else break;
140: else if ( (Y == 1) || (Y == -1) ) break;
141: */
142: if ( ( (Y = T - X) == 0 ) || (Y == 1) || (Y == -1) ) break;
143: X = T;
144: }
145: }
146: return (X/(10^I));
147: }
148: #endif
149:
150: def sqrt(A) {
151: for ( J = 0, T = A; T >= 2^27; J++ ) {
152: T = idiv(T,2^27)+1;
153: }
154: for ( I = 0; T >= 2; I++ ) {
155: S = idiv(T,2);
156: if ( T = S+S )
157: T = S;
158: else
159: T = S+1;
160: }
161: X = (2^27)^idiv(J,2)*2^idiv(I,2);
162: while ( 1 ) {
163: if ( (Y=X^2) < A )
164: X += X;
165: else if ( Y == A )
166: return X;
167: else
168: break;
169: }
170: while ( 1 )
171: if ( (Y = X^2) <= A )
172: return X;
173: else
174: X = idiv(A + Y,2*X);
175: }
176: end$