/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C GENERATE FAREY SERIES C 06/10/07 (DKC) C C This subroutine generates the fractions in a Farey series segment after C a given initial fraction. The numerator and denominator of a current C fraction ("NUM" and "DEN") and the numerators and denominators of fractions C on a stack are maintained with the numerator and denominator of the C fraction on the top of the stack being denoted by "STACK(SP,1)" and C "STACK(SP,2)" respectively (where "SP" denotes the stack pointer). 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. An adjacent pair of fractions C initially on the stack, say n/d and n'/d', must satisfy d*n' - n*d' = 1. C Also, initially DEN*STACK(SP,1) - NUM*STACK(SP,2) must equal 1. If C 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 then set to these values.) If J > 1, NEWNUM and NEWDEN are C 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 last C fraction has been reached. If not, another iteration of the algorithm is C performed. The algorithm can be iterated since the stack is always C maintained so that two adjacent fractions on the stack, say n/d and n'/d', C satisfy d*n' - n*d' = 1. The execution time of the algorithm is C proportional to the number of fractions in the Farey series since all the C operations result in a successive fraction being generated (and the number C of operations required to generate a fraction is less than some constant) C or a fraction being pushed on the stack (not all the fractions are pushed C on the stack and if they are, only once). Other than a relatively small C array of (n - 1)*2 elements, the algorithm requires no more memory than C that required to hold the fractions. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/ unsigned int frank(unsigned int N, unsigned int M, unsigned int *R, unsigned int NUM, unsigned int DEN, unsigned int SP, unsigned int *STACK) { unsigned int I,J,K,DELNUM,DELDEN; // // CHECK FOR AN INCREASE IN DENOMINATOR // L100: J=(N-STACK[2*SP+1])/DEN; if(J==0) goto L500; // // INCREASE IN DENOMINATOR // if(J==1) goto L300; 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 L200; // // PUSH FRACTIONS ON THE STACK // K=J-2; for (I=0; I<K; I++) { SP=SP+1; STACK[2*SP]=NUM; STACK[2*SP+1]=DEN; NUM=NUM+DELNUM; DEN=DEN+DELDEN; } // // DECREASE IN DENOMINATOR // L200: M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; goto L100; // // INCREASE IN DENOMINATOR // L300: DELNUM=STACK[2*SP]; DELDEN=STACK[2*SP+1]; NUM=DELNUM+NUM; DEN=DELDEN+DEN; K=(N-DEN)/DELDEN; if(K==0) goto L400; for (I=0; I<K; I++) { M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; NUM=NUM+DELNUM; DEN=DEN+DELDEN; } L400: M=M+1; R[2*M]=NUM; R[2*M+1]=DEN; // // DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK // L500: 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; // STACK[0]=NUM; STACK[1]=DEN; return(M); }