1 #include "Global.h"
2 #include "Rat.h"
3
rI(Long a)4 Rat rI(Long a) /* a -> a/1 */
5 { Rat c; c.N=a; c.D=1; return c;
6 }
rR(Long a,Long b)7 Rat rR(Long a, Long b) /* (a,b) -> a/b */
8 { Rat c; register Long g=Fgcd(a,b);
9 if( (c.D=(b/g)) < 0 ) {c.D=-c.D; c.N=-a/g;} else c.N=a/g;
10 #ifdef TEST
11 if(c.D<=0) {Rpr(c); puts(" *** c.D<=0 in rR! ***"); exit(0);}
12 #endif
13 return c;
14 }
irP(Long i,Rat c)15 Rat irP(Long i, Rat c) /* (int) * (Rat) */
16 { register Long g=Fgcd(i,c.D); if(g < 0) g=-g;
17 c.N *= (i/g); c.D/=g;
18 #ifdef TEST
19 if(c.D<=0) {Rpr(c); puts(" *** c.D<=0 in irP! ***"); exit(0);}
20 #endif
21 return c;
22 }
rS(Rat a,Rat b)23 Rat rS(Rat a, Rat b) /* a + b */
24 { Rat c; register Long g=Fgcd(a.D,b.D);
25 g = Fgcd(c.N=a.N*(b.D/g)+b.N*(a.D/g), c.D=a.D*(b.D/g));
26 if(g < 0) g=-g;
27 c.N/=g; c.D/=g; /* LONG in line above ... */
28 #ifdef TEST
29 if(c.D<=0) {Rpr(c); puts(" *** c.D<=0 in rS! ***"); exit(0);}
30 #endif
31 return c;
32 }
rD(Rat a,Rat b)33 Rat rD(Rat a, Rat b) /* a - b */
34 { Rat c; register Long g=Fgcd(a.D,b.D);
35 g = Fgcd(c.N=a.N*(b.D/g)-b.N*(a.D/g), c.D=a.D*(b.D/g));
36 if(g < 0) g=-g;
37 c.N/=g; c.D/=g; /* LONG in line above ... */
38 #ifdef TEST
39 if(c.D<=0) {Rpr(c); puts(" *** c.D<=0 in rD! ***"); exit(0);}
40 #endif
41 return c;
42 }
rP(Rat a,Rat b)43 Rat rP(Rat a, Rat b) /* a * b */
44 { register Long g=Fgcd(a.N,b.D); register Long h=Fgcd(b.N,a.D);
45 Rat c; c.N=(a.N/g)*(b.N/h);
46 if((c.D=(a.D/h)*(b.D/g)) < 0) {c.D=-c.D; c.N=-c.N;} return c;
47 }
rQ(Rat a,Rat b)48 Rat rQ(Rat a, Rat b) /* a / b */
49 { register Long g=NNgcd(a.N,b.N); register Long h=Fgcd(b.D,a.D);
50 Rat c; c.N=(a.N/g)*(b.D/h);
51 if((c.D=(a.D/h)*(b.N/g)) < 0) {c.N=-c.N; c.D=-c.D;} return c;
52 }
rC(Rat a,Rat b)53 int rC(Rat a, Rat b) /* Compare = [1 / 0 / -1] iff a [gt/eq/lt] b */
54 { register Long g=Fgcd(a.D,b.D);
55 if((g=a.N*(b.D/g)-b.N*(a.D/g))) {if(g>0) return 1; else return -1;}
56 else return 0;
57 }
Rpr(Rat c)58 void Rpr(Rat c) /* write "c.N/c.D" -> outFN */
59 { fprintf(outFILE,"%d/%d",(int) c.N, (int) c.D);
60 }
61 #ifdef TEST
TEST_Rat()62 void TEST_Rat()
63 { Rat a, b; int i, j, k, l;
64 puts(" type 4 numbers: "); scanf("%d%d%d%d",&i,&j,&k,&l);
65 Rpr(rS(rR(i,j),rR(k,l))); puts(" i/j + k/l");
66 Rpr(rD(rR(i,j),rR(k,l))); puts(" i/j - k/l");
67 }
68 #endif
69 /* ========== END of Rational Operations ========== */
70
71
72
73 /* ========== (Extended) Greatest Common Divisor ========= */
74 /*
75 * Fgcd() computes the GCD for two positive integers (F -> fast)
76 * NNgcd() computes the non-negative GCD for two arbitrary integers
77 *
78 * Egcd() computes the EGCD for two integers
79 * REgcd() recursively computes the EGCD for a number of integers
80 */
Fgcd(register Long a,register Long b)81 Long Fgcd(register Long a, register Long b) /* Fast greatest common div */
82 {
83 while( a %= b ) if ( !(b %= a) ) return a;
84 return b;
85 }
NNgcd(register Long a,register Long b)86 Long NNgcd(register Long a, register Long b) /* NonNegative gcd handling 0 */
87 {
88 a = (a<0) ? -a:a; b = (b<0) ? -b:b; if (!b) return a;
89 while( a %= b ) if ( !(b %= a) ) return a;
90 return b;
91 }
92 /* ========== Egcd(a,b,*x,*y)=a *x+b *y; REgcd(Vin, *dim, Vout) = Vin.Vout;
93 *
94 * a_i+1=a_i-1 % a_i, a_I+1==0 => egcd(a_0,a_1)=a_I; Vout={x_I,y_I}
95 * a_i==x_ia_0+y_ia_1a_i+1 and a_i+1=a_i-1 -(a_i-1 /a_i)a_i =>
96 * x_i+1=x_i-1 -(a_i-1 /a_i)x_i with x_0=1; x_1=0; y_I=(a_I-a_0 x_I)/a_1
97 * */
Egcd(register Long A0,register Long A1,Long * Vout0,Long * Vout1)98 Long Egcd(register Long A0, register Long A1, Long *Vout0, Long *Vout1)
99 { register Long V0=A0, V1=A1, A2, X0=1, X1=0, X2=0;
100 while((A2 = A0 % A1)) { X2=X0-X1*(A0/A1); A0=A1; A1=A2; X0=X1; X1=X2; }
101 *Vout0=X1, *Vout1=(A1-(V0) * X1)/ (V1); return A1;
102 }
103 /*
104 * g = v0 x0 + ... + v(n-1) x(n-1), G= A g+ B vn =>
105 * G = v0 x0 + ... + vn xn with xn=B, xi*=A,
106 * */
REgcd(Long * Vin,int * _n,Long * Vout)107 Long REgcd(Long *Vin, int *_n, Long *Vout) /* recursive Egcd(a_1,...,a_n) */
108 { Long Ain[2], Aout[2], gcd; int N=*_n-1, i;
109 if (*_n==2) return Egcd(*Vin,(Vin[1]),Vout,&(Vout[1]));
110 *Ain=REgcd(Vin,&N,Vout); Ain[1]=Vin[N]; gcd=
111 Egcd(*Ain,(Ain[1]),Aout,&(Aout[1])); Vout[N]=Aout[1];
112 for(i=0; i<N; i++) Vout[i] *= (*Aout);
113 return gcd;
114 }
115 /* ========== END of (Extended) Greatest Common Divisor ========= */
116
117
118 /* ========== GLZ-matrix=(Egcd, B_1, B_2, ...) E.W=g, B.W=0 ========= */
119 /*
120 * assumes W[i]!=0, lines 1 thru d-1 are lower triangular and *
121 * their diagonal entries are positive if all W[i] are positive. *
122 * GLZ has to be an initialized list of pointers. returns gcd(W_i) *
123 * */
FloorQ(Long N,Long D)124 Long FloorQ(Long N,Long D) /* F <= N/D < F+1 */
125 { Long F; if(D<0) {D=-D; N=-N;} F=N/D; return (F*D>N) ? F-1 : F;
126 }
CeilQ(Long N,Long D)127 Long CeilQ(Long N,Long D) { return -FloorQ(-N,D); }
RoundQ(Long N,Long D)128 Long RoundQ(Long N,Long D)
129 { Long F; if(D<0) {D=-D; N=-N;} F=N/D; return F+(2*(N-F*D))/D;
130 }
W_to_GLZ(Long * W,int * d,Long ** GLZ)131 Long W_to_GLZ(Long *W, int *d, Long **GLZ)
132 { int i, j; Long G, *E=*GLZ, *B=GLZ[1]; for(i=0;i<*d;i++) assert(W[i]!=0);
133 for(i=1;i<*d;i++)for(j=0;j<*d;j++)GLZ[i][j]=0;
134 G=Egcd(W[0],W[1],&E[0],&E[1]); B[0]=-W[1]/G; B[1]=W[0]/G;
135 for(i=2;i<*d;i++)
136 { Long a, b, g=Egcd(G,W[i],&a,&b); B=GLZ[i];
137 B[i]= G/g; G=W[i]/g; for(j=0;j<i;j++) B[j]=-E[j]*G; /* B=Base-Line */
138 for(j=0;j<i;j++) E[j]*=a;
139 E[j]=b; /* next Egcd */
140 for(j=i-1;0<j;j--) /* I M P R O V E M E N T */
141 { int n; Long *Y=GLZ[j], rB=RoundQ(B[j],Y[j]), rE=RoundQ(E[j],Y[j]);
142 /* rB=CeilQ(B[j],Y[j]), rE=CeilQ(E[j],Y[j]); */
143 /* printf(" [%d/%d -> %d] ",(int)B[j],(int)Y[j],(int)rB);
144 printf(" [%d/%d -> %d] ",(int)E[j],(int)Y[j],(int)rE); */
145 for(n=0;n<=j;n++) { B[n] -= rB*Y[n]; E[n] -= rE*Y[n]; }
146 } G=g;
147 }
148 return G;
149 }
150 /* Map Permutations: Do "ArgFun" for all permutations pi of *d elements */
151 #ifdef __cplusplus
152 #define ARG_FUN void (*ArgFun)(int *d,int *pi,int *pinv,void *info)
153 #else
154 #define ARG_FUN void ( ArgFun() )
155 #endif
Map_Permut(int * d,int * pi,int * pinv,ARG_FUN,void * AuxPtr)156 void Map_Permut(int *d,int *pi,int *pinv,ARG_FUN,void *AuxPtr)
157 { int i, j, n_rem_perm, n_perm=1, a, b, perm_j;
158 for (i=1;i<=*d;i++) n_perm*=i;
159 for (i=0;i<n_perm;i++)
160 { b=i; n_rem_perm=n_perm; for (j=0;j<*d;j++) pinv[j]=-1;
161 for (j=*d;j>0;j--)
162 { n_rem_perm/=j; a=b/n_rem_perm; b=b%n_rem_perm; perm_j=*d;
163 while (a>=0) { perm_j--; if (pinv[perm_j]<0) a--;}
164 pi[j-1]=perm_j; pinv[perm_j]=j-1;
165 } (ArgFun)(d,pi,pinv,AuxPtr);
166 }
167 }
168 /* ========== END of GLZ-matrix=(Egcd, B_1, B_2, ...) ========= */
169
170
171
172
LrI(LLong a)173 LRat LrI(LLong a) /* a -> a/1 */
174 { LRat c; c.N=a; c.D=1; return c;
175 }
LrR(LLong a,LLong b)176 LRat LrR(LLong a, LLong b) /* (a,b) -> a/b */
177 { LRat c; register LLong g=LFgcd(a,b);
178 if( (c.D=(b/g)) < 0 ) {c.D=-c.D; c.N=-a/g;} else c.N=a/g;
179 #ifdef TEST
180 if(c.D<=0) {LRpr(c); puts(" *** c.D<=0 in rR! ***"); exit(0);}
181 #endif
182 return c;
183 }
LirP(LLong i,LRat c)184 LRat LirP(LLong i, LRat c) /* (int) * (LRat) */
185 { register LLong g=LFgcd(i,c.D); if(g < 0) g=-g;
186 c.N *= (i/g); c.D/=g;
187 #ifdef TEST
188 if(c.D<=0) {LRpr(c); puts(" *** c.D<=0 in irP! ***"); exit(0);}
189 #endif
190 return c;
191 }
LrS(LRat a,LRat b)192 LRat LrS(LRat a, LRat b) /* a + b */
193 { LRat c; register LLong g=LFgcd(a.D,b.D);
194 g = LFgcd(c.N=a.N*(b.D/g)+b.N*(a.D/g), c.D=a.D*(b.D/g));
195 if(g < 0) g=-g;
196 c.N/=g; c.D/=g; /* LLong in line above ... */
197 #ifdef TEST
198 if(c.D<=0) {LRpr(c); puts(" *** c.D<=0 in rS! ***"); exit(0);}
199 #endif
200 return c;
201 }
LrD(LRat a,LRat b)202 LRat LrD(LRat a, LRat b) /* a - b */
203 { LRat c; register LLong g=LFgcd(a.D,b.D);
204 g = LFgcd(c.N=a.N*(b.D/g)-b.N*(a.D/g), c.D=a.D*(b.D/g));
205 if(g < 0) g=-g;
206 c.N/=g; c.D/=g; /* LLong in line above ... */
207 #ifdef TEST
208 if(c.D<=0) {LRpr(c); puts(" *** c.D<=0 in rD! ***"); exit(0);}
209 #endif
210 return c;
211 }
LrP(LRat a,LRat b)212 LRat LrP(LRat a, LRat b) /* a * b */
213 { register LLong g=LFgcd(a.N,b.D); register LLong h=LFgcd(b.N,a.D);
214 LRat c; c.N=(a.N/g)*(b.N/h);
215 if((c.D=(a.D/h)*(b.D/g)) < 0) {c.D=-c.D; c.N=-c.N;} return c;
216 }
LrQ(LRat a,LRat b)217 LRat LrQ(LRat a, LRat b) /* a / b */
218 { register LLong g=LNNgcd(a.N,b.N); register LLong h=LFgcd(b.D,a.D);
219 LRat c; c.N=(a.N/g)*(b.D/h);
220 if((c.D=(a.D/h)*(b.N/g)) < 0) {c.N=-c.N; c.D=-c.D;} return c;
221 }
LrC(LRat a,LRat b)222 int LrC(LRat a, LRat b) /* Compare = [1 / 0 / -1] iff a [gt/eq/lt] b */
223 { register LLong g=LFgcd(a.D,b.D);
224 if((g=a.N*(b.D/g)-b.N*(a.D/g))) {if(g>0) return 1; else return -1;}
225 else return 0;
226 }
LRpr(LRat c)227 void LRpr(LRat c) /* write "c.N/c.D" -> outFN */
228 { fprintf(outFILE,"%lld/%lld",(LLong) c.N, (LLong) c.D);
229 }
230 #ifdef TEST
TEST_LRat()231 void TEST_LRat()
232 { LRat a, b; int i, j, k, l;
233 puts(" type 4 numbers: "); scanf("%d%d%d%d",&i,&j,&k,&l);
234 LRpr(LrS(rR(i,j),LrR(k,l))); puts(" i/j + k/l");
235 LRpr(LrD(LrR(i,j),LrR(k,l))); puts(" i/j - k/l");
236 }
237 #endif
238 /* ========== END of LRational OpeLRations ========== */
239
240
241
242 /* ========== (Extended) Greatest Common Divisor ========= */
243 /*
244 * LFgcd() computes the GCD for two positive integers (F -> fast)
245 * LNNgcd() computes the non-negative GCD for two arbitrary integers
246 *
247 * LEgcd() computes the LEgcd for two integers
248 * LREgcd() recursively computes the LEgcd for a number of integers
249 */
LFgcd(register LLong a,register LLong b)250 LLong LFgcd(register LLong a, register LLong b) /* Fast greatest common div */
251 {
252 while( a %= b ) if ( !(b %= a) ) return a;
253 return b;
254 }
LNNgcd(register LLong a,register LLong b)255 LLong LNNgcd(register LLong a, register LLong b) /* NonNegative gcd handling 0 */
256 {
257 a = (a<0) ? -a:a; b = (b<0) ? -b:b; if (!b) return a;
258 while( a %= b ) if ( !(b %= a) ) return a;
259 return b;
260 }
261 /* ========== LEgcd(a,b,*x,*y)=a *x+b *y; LREgcd(Vin, *dim, Vout) = Vin.Vout;
262 *
263 * a_i+1=a_i-1 % a_i, a_I+1==0 => LEgcd(a_0,a_1)=a_I; Vout={x_I,y_I}
264 * a_i==x_ia_0+y_ia_1a_i+1 and a_i+1=a_i-1 -(a_i-1 /a_i)a_i =>
265 * x_i+1=x_i-1 -(a_i-1 /a_i)x_i with x_0=1; x_1=0; y_I=(a_I-a_0 x_I)/a_1
266 * */
LEgcd(register LLong A0,register LLong A1,LLong * Vout0,LLong * Vout1)267 LLong LEgcd(register LLong A0, register LLong A1, LLong *Vout0, LLong *Vout1)
268 { register LLong V0=A0, V1=A1, A2, X0=1, X1=0, X2=0;
269 while((A2 = A0 % A1)) { X2=X0-X1*(A0/A1); A0=A1; A1=A2; X0=X1; X1=X2; }
270 *Vout0=X1, *Vout1=(A1-(V0) * X1)/ (V1); return A1;
271 }
272 /*
273 * g = v0 x0 + ... + v(n-1) x(n-1), G= A g+ B vn =>
274 * G = v0 x0 + ... + vn xn with xn=B, xi*=A,
275 * */
LREgcd(LLong * Vin,int * _n,LLong * Vout)276 LLong LREgcd(LLong *Vin, int *_n, LLong *Vout) /* recursive LEgcd(a_1,...,a_n) */
277 { LLong Ain[2], Aout[2], gcd; int N=*_n-1, i;
278 if (*_n==2) return LEgcd(*Vin,(Vin[1]),Vout,&(Vout[1]));
279 *Ain=LREgcd(Vin,&N,Vout); Ain[1]=Vin[N]; gcd=
280 LEgcd(*Ain,(Ain[1]),Aout,&(Aout[1])); Vout[N]=Aout[1];
281 for(i=0; i<N; i++) Vout[i] *= (*Aout);
282 return gcd;
283 }
284 /* ========== END of (Extended) Greatest Common Divisor ========= */
285