/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C GENERATE FAREY SERIES C 06/10/07 (Darrell Cox) C C This C program computes the first half of a Farey series of order n (the C fractions in a Farey series are complementary around the fraction 1/2, C that is, if 1/2 is the i th fraction in the Farey series, then F(i+1) = C 1 - F(i-1), F(i+2) = 1 - F(i-2), etc.). Temporary results are pushed on C the stack ("STACK"). The size of this array must be at least n+1. The C numerators and denominators of the fractions in the first half of the C Farey series are output in the array "R". The size of this array must C be as large as the number of fractions in half a Farey series (about C .1520*n*n). C C This algorithm for generating a Farey series uses Haros' theorem C that if h/k and h'/k' are successive fractions in a Farey series, then C kh' - hk' = 1. The numerator and denominator of a current fraction C ("NUM" and "DEN") and the numerators and denominators of fractions on a C stack are maintained with the numerator and denominator of the fraction C on the top of the stack being denoted by "STACK(SP,1)" and "STACK(SP,2)" C respectively (where "SP" denotes the stack pointer). The initial values C of "NUM" and "DEN" are 1 and N respectively where N is the order of the C Farey series. The initial values of the numerators and denominators of C the fractions on the stack are 1 for the numerators and N-1, N-2, N-3,..., C 2 for the denominators (with N-1 being on the top of the stack). As the C numerators and denominators of the successive fractions in the Farey series C are computed, they are stored into the array "R(M,1)" and "R(M,2)" where C "M" is the current index into the array. (Actually, as implemented here, C the numerators and denominators of the fractions 1/N, 1/(N-1), 1/(N-2),..., C 1/J where J=N-[N/2] are first stored in the array "R", the index "M" is set C to [N/2]+1, "NUM" and "DEN" are set to 1 and "J" respectively, and the C numerators and denominators of the rest of the fractions with a numerator C of 1 are pushed on the stack. The only reason for doing this is that the C processing for the first few fractions can be simplified and thus a little C execution time saved.) Since the numerators of the fractions initially on C the stack are 1 and the denominators are numerically consecutive and C decreasing, an adjacent pair of fractions on the stack, say n/d and n'/d', C satisfies d*n' - n*d' = 1. Also, initially DEN*STACK(SP,1)-NUM*STACK(SP,2) C = 1. If also DEN + STACK(SP,2) > N (using the theorem that the sum of the C denominators of two successive fractions in a Farey series is greater than C the order of the Farey series), then by Haros' theorem (that two successive C fractions h/k and h'/k' in a Farey series satisfy k*h' - h*k' = 1), C "STACK(SP,1)" and "STACK(SP,2)" are the numerator and denominator of the C next fraction in the Farey series. (There are infinitely many solutions C (x,y) of k*x - h*y = 1, (h,k)=1, but only one solution where x <= N, C y <= N, and y + k > N.) The modulus function is used to find a linear C combination of "DEN" and STACK(SP,2)", denoted by "NEWDEN", such that C N - NEWDEN < DEN. If "NEWDEN" is set to "STACK(SP,2) + [(N - STACK(SP,2))/ C DEN]*DEN", then N - NEWDEN < DEN (since "N - NEWDEN" is "(N - STACK(SP,2)) C modulus DEN)". Denote "[(N - STACK(SP,2))/DEN]" (where the brackets denote C the greatest integer function) by "J" and denote "STACK(SP,1) + NUM*J" C by "NEWNUM". Adding "DEN*NUM*J - NUM*DEN*J" (= 0) to both sides of DEN* C STACK(SP,1) - NUM*STACK(SP,2) = 1 gives DEN*NEWNUM - NUM*NEWDEN = 1 and C hence "NEWNUM/NEWDEN" is the next fraction in the Farey series. (If J=0, C "NEWNUM" and "NEWDEN" are "STACK(SP,1)" and "STACK(SP,2)". These values C are popped off the stack and stored in the "R" array. "NUM" and "DEN" C are set to these values. It is not possible that the stack [even after it C has been popped] be empty when "J" equals 0 [this will be elaborated later], C so checking the stack pointer is not necessary.) If J > 1, "NEWNUM" and C "NEWDEN" are stored in the "R" array. If J > 1, [STACK(SP,2) + DEN*I]* C [STACK(SP,1) + NUM*(I-1)] - [STACK(SP,1) + NUM*I]*[STACK(SP,2) + DEN* C (I-1)] = 1, (I=J,J-1,J-2,...,2), (since DEN*STACK(SP,1) - NUM*STACK(SP,2) C = 1). Let "OLDSP" equal "SP". If J > 1, the numerators and denominators C "STACK(OLDSP,1) + NUM*I" and "STACK(OLDSP,2) + DEN*I", (I=1,2,3,...,(J-1)), C are pushed on the stack (so that the fraction "[STACK(OLDSP,1) + NUM*(J-1)]/ C [STACK(OLDSP,2) + DEN*(J-1)]" is on the top of the stack). [STACK(OLDSP,2) C + DEN]*STACK(OLDSP,1) - [STACK(OLDSP,1) + NUM]*STACK(OLDSP,2) = 1, so any C adjacent pair of fractions on the stack, say n/d and n'/d', still satisfies C d*n' - n*d' = 1. If J > 1, DEN <= DEN*(J-1) and hence STACK(OLDSP,2) + C DEN*J + DEN < 2*STACK(OLDSP,2) + DEN*J + DEN <= 2*STACK(OLDSP,2) + DEN*J + C DEN*(J-1). Since N - [STACK(OLDSP,2) + DEN*J] < DEN, N < STACK(OLDSP,2) + C DEN*J + DEN, and hence N < 2*STACK(OLDSP,2) + DEN*J + DEN*(J-1), that is, C N - [STACK(OLDSP,2) + DEN*(J-1)] < STACK(OLDSP,2) + DEN*J. Therefore if J > C 1, NEWDEN*STACK(SP,1) - NEWNUM*STACK(SP,2) = 1 and N - STACK(SP,2) < NEWDEN, C and hence the fraction on the top of the stack is the next fraction in the C Farey series. (In practice, it is not necessary to push this fraction on C the stack in the first place since it is known that it will be the next C fraction in the series. If J=2, it is not necessary to push any fractions C on the stack.) If J=1, NEWNUM = STACK(SP,1) + NUM and NEWDEN = STACK(SP,2) C + DEN. "NEWNUM" and "NEWDEN" are then the numerator and denominator of the C next fraction in the Farey series and are stored in the "R" array. Denote C the increase in the denominator by "DELDEN". If J=1, "DELDEN" equals C "STACK(SP,2)" (since NEWDEN - DEN = STACK(SP,2)). Let K = [(N - NEWDEN)/ C DELDEN] ("K" is the maximum number of times "NEWDEN" can be increased by C "DELDEN" and still remain less than or equal to "N"). If J=1 and K > 0, C [STACK(SP,2) + DEN*I]*[STACK(SP,1) + NUM*(I+1)] - [STACK(SP,1) + NUM*I]* C [STACK(SP,2) + DEN*(I+1)] = 1, (I=1,2,3,...,K), (since DEN*STACK(SP,1) - C NUM*STACK(SP,2) = 1). Since STACK(SP,2) < DEN and N - [STACK(SP,2) + DEN] < C DEN, N - [STACK(SP,2) + DEN*(I+1)] < [STACK(SP,2) + DEN*I], (I=1,2,3,...,K), C and hence (if K>0) "[STACK(SP,1) + NUM*(I+1)]" and "[STACK(SP,2) + DEN* C (I+1)]", (I=1,2,3,...,K) are the numerators and denominators of the next "K" C successive fractions in the Farey series. Now denote "DEN + STACK(SP,2)* C (K+1)" by "NEWDEN" and "NUM + STACK(SP,1)*(K+1)" by "NEWNUM". (Note that if C K=0, the definition of "NEWDEN" and "NEWNUM" is the same as before.) Adding C "STACK(SP,2)*(K+1)*STACK(SP,1) - STACK(SP,1)*(K+1)*STACK(SP,2)" (=0) to both C sides of DEN*STACK(SP,1) - NUM*STACK(SP,2) = 1 gives NEWDEN*STACK(SP,1) - C NEWNUM*STACK(SP,2) = 1. N - STACK(SP,2) = DEN + [N - (STACK(SP,2) + DEN)] < C DEN + STACK(SP,2)*{[(N - (STACK(SP,2) + DEN))/STACK(SP,2)] + 1} = DEN + C STACK(SP,2)*(K+1) and hence N - STACK(SP,2) < NEWDEN. Therefore if J=1, the C numerator and denominator of the next fraction in the Farey series are C "STACK(SP,1)" and "STACK(SP,2)". "NUM" is set to "STACK(SP,1)" and "DEN" C is set to "STACK(SP,2)" and these values are popped off the stack and stored C in the "R" array. The stack pointer is then checked to see if the fraction C 1/2 has been reached (the fraction 1/2 can only be reached when J=1, so C checking the stack pointer in other instances is unnecessary). If not, C another iteration of the algorithm is performed. The algorithm can be C iterated since the stack is always maintained so that two adjacent fractions C on the stack, say n/d and n'/d', satisfy d*n' - n*d' = 1. The execution C time of the algorithm is proportional to the number of fractions in the C Farey series since all the operations result in a successive fraction being C generated (and the number of operations required to generate a fraction is C less than some constant) or a fraction being pushed on the stack (not all C the fractions are pushed on the stack and if they are, only once). Other C than a relatively small array of [(n - 1)/2]*2 elements (where the brackets C denote the greatest integer function), the algorithm requires no more C memory than that required to hold the fractions. This technique can also C be used to compute the fractions in a Farey series between the fractions C 1/a and 1/b where a>b. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/ unsigned int farey(unsigned int N, unsigned int *R, unsigned int *STACK) { unsigned int M,NUM,DEN,DELNUM,DELDEN,I,J,K,SP; // // STORE THE FIRST FRACTION // R[0]=0; R[1]=1; // // STORE THE FRACTIONS BEFORE THE FIRST INCREASE IN DENOMINATOR // K=(N/2)+2; J=N; for (I=1; I<K; I++) { R[2*I]=1; R[2*I+1]=J; J=J-1; } M=K-1; NUM=R[2*M]; DEN=R[2*M+1]; // // PUSH THE REMAINING FRACTIONS WITH A NUMERATOR OF 1 ON THE STACK // K=N-K+1; if(K<2) goto L600; SP=0; for (I=2; I<=K; I++) { SP=SP+1; STACK[2*SP]=1; STACK[2*SP+1]=I; } /***************/ /* BEGIN LOOP */ /***************/ // // CHECK FOR AN INCREASE IN DENOMINATOR // L100: J=(N-STACK[2*SP+1])/DEN; if(J!=0) goto L200; // // DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK // NUM=STACK[2*SP]; DEN=STACK[2*SP+1]; SP=SP-1; M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; goto L100; // // INCREASE IN DENOMINATOR // L200: if(J==1) goto L400; DELNUM=NUM; DELDEN=DEN; NUM=STACK[2*SP]+NUM*J; DEN=STACK[2*SP+1]+DEN*J; M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; NUM=STACK[2*SP]+DELNUM; DEN=STACK[2*SP+1]+DELDEN; if(J==2) goto L300; // // PUSH FRACTIONS ON THE STACK // K=J-2; for (I=1; I<=K; I++) { SP=SP+1; STACK[2*SP]=NUM; STACK[2*SP+1]=DEN; NUM=NUM+DELNUM; DEN=DEN+DELDEN; } // // DECREASE IN DENOMINATOR // L300: M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; goto L100; // // INCREASE IN DENOMINATOR // L400: DELNUM=STACK[2*SP]; DELDEN=STACK[2*SP+1]; NUM=DELNUM+NUM; DEN=DELDEN+DEN; K=(N-DEN)/DELDEN; if(K==0) goto L500; for (I=1; I<=K; I++) { M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; NUM=NUM+DELNUM; DEN=DEN+DELDEN; } L500: M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; // // DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK // NUM=STACK[2*SP]; DEN=STACK[2*SP+1]; SP=SP-1; M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; if(SP>0) goto L100; /***************/ /* END LOOP */ /***************/ L600: return(M+1); }