// FareySeq.java
// Fanchen Bao, 12/13/2018
// https://raw.githubusercontent.com/FanchenBao/MATH_Farey_sequence/master/FareySeq.java
// Slightly modified by Andrew Davison, ad@coe.psu.ac.th, April 2023


/* Produce a Farey sequence with a upper denominator maxDem 
   and in a specified range (the default is [0, 1]). 

   The algorithm is based on three principles:
   1. Two reduced fractions 0 <= a/b < c/d <= 1 are neighbors in a
      Farey sequence iff b*c - a*d == 1.

   2. If b*c - a*d == 1, then the unique fraction between a/b and c/d with
      the smallest denominator is (a+b)/(c+d).

   3. If a/b < c/d < e/f are three consecutive Farey sequence fractions,
      then c/d == (a+e)/(b+f)

   https://en.wikipedia.org/wiki/Farey_sequence

   Usage:
     > compile FareySeq.java
     > run FareySeq

   See the main() test rig for how to use this class.
   This class is also used in FareyGrid.java and FordCircles.java.
*/


import java.util.*;

import org.apache.commons.numbers.fraction.*;
import org.apache.commons.numbers.core.*;


public class FareySeq
{
  private ArrayList<Fraction> fareys;  


  public FareySeq(int maxDem)
  /* Construct a Farey sequence with maxDem being the max 
     denominator in the range [0, 1]
  */
  { fareys = new ArrayList<Fraction>();
    makeSeq(maxDem, Fraction.of(0, 1), Fraction.of(1, 1));
  }


  public FareySeq(int maxDem, int beginN, int beginD, int endN, int endD)
  /* Construct a Farey sequence with maxDem being the max denominator 
     in the range [beginN/beginD, endN/endD]. 
     The boundaries must satisfy the conditions: 
     0 <= beginN/beginD <= endN/endD <= 1; 
     Also, after putting both boundaries into proper reduced fraction 
     form, then max(beginD, endD) <= maxDem
  */
  { fareys = new ArrayList<Fraction>();

    // turn begin and end into proper fraction
    int beginGCD = Nut.gcd(beginN, beginD);
    int endGCD = Nut.gcd(endN, endD);
    Fraction lower = Fraction.of(beginN/beginGCD, beginD/beginGCD);
    Fraction upper = Fraction.of(endN/endGCD, endD/endGCD);

    if (!checkBound(maxDem, lower, upper))
      System.exit(0);

    makeSeq(maxDem, lower, upper);
  }  // end of FareySeq()


  public ArrayList<Fraction> getFareySeq()
  {  return fareys;   }


  public int getSize()
  {  return fareys.size();   }


  public long[] getFrac(int pos)
  /* Return the [numerator, denominator] pair of the
     fraction in the Farey sequence at this position
  */
  {
    if ((pos >= fareys.size()) || (pos < 0)) {
      System.out.println("Error; index out of range.");
      return null;
    }

    long[] res = new long[2];
    res[0] = fareys.get(pos).getNumerator();
    res[1] = fareys.get(pos).getDenominator();
    return res;
  } /* getFrac */



  private boolean checkBound(int maxDem, Fraction lower,
                                         Fraction upper)
  /* Return true if upper and lower bound are in the
     range of [0, 1], lower bound < upper bound,
     and the denominator in the most reduced form of both
     upper and lower bound does not exceed maxDem
  */
  {
    if (lower.compareTo(upper) > 0) {
      System.out.println("Lower bound must be <= upper bound");
      return false;
    }

    if (lower.compareTo(Fraction.of(0, 1)) <= 0) {
      System.out.println("Lower bound must be > 0");
      return false;
    }

    if (upper.compareTo(Fraction.of(1, 1)) > 0) {
      System.out.println("Upper bound must be <= 1");
      return false;
    }

    if ((lower.getDenominator() > maxDem) ||
        (upper.getDenominator() > maxDem) ) {
      System.out.println("Denominators of lower or upper cannot exceed " + maxDem);
      return false;
    }

    return true;
  } /* checkBound */



  private void makeSeq(int maxDem, Fraction lower, Fraction upper)
  /* Produce a Farey sequence.
     Its fractions will range from lower to upper,
     with a maximum denominator == maxDem
  */
  {
    fareys.add(lower);

    int a, b;   // left-most fraction in the sequence
    if (lower.getNumerator() == 0) {    
        // when lower is 0, we must use the next
        // fraction as the starting point
      fareys.add(Fraction.of(1, maxDem));
      a = 1;
      b = maxDem;
    }
    else {
      a = lower.getNumerator();
      b = lower.getDenominator();
    }

    int c0 = a;
    int d0 = b - 1/a;   
    /* Use any fraction that satisfies b*c0 - a*d0 == 1 and 
       c0 <= a, d0 <= b. We choose c0 = a, d0 == b - 1/a for
       convenience. The purpose of c0/d0 is to
       calculate the first c/d next to a/b. Once we have 
       the first adjacent pair of a/b and c/d, then all 
       the remaining fractions can be generated
    */

    /* Generate the first c/d, the next fraction to 
       the right of a/b. The equation is acquired from 
       b*c0 - a*d0 == 1 == b*c - a*d, given gcd(a, b) == 1 
       and d <= maxDem < b+d
    */
    int k = (maxDem - d0) / b;     // integer division
    int c = k * a + c0;
    int d = k * b + d0;

    /* Generate the rest of the sequence with the given maxDem and
       each next fraction is next to the previous one
       until c = 1, d = 2 (end condition). 
    */

    int e, f;     // the next fraction to c/d

    /* Since a/b < c/d < e/f, and all three are adjacent in the 
       current sequence, then we have
       (a + e)/(b + f) == c/d. Since gcd(c,d) == 1, then
       we have a+e == k*c and b+f == k*d.
       Since c/d and e/f are adjacent, then the mediant between 
       them must have a denominator (d + f) > maxDem.
	  
       Thus, (b + f) == k*d <= b + maxDem < (b + d + f) == 
             d + k*d = (k + 1)*d
       k*d <= b + maxDem < (k + 1)*d
       k <= (b + maxDem) / d < k + 1
       So k == (b + maxDem) / d with integer division
    */

    while (c != upper.getNumerator() || 
           d != upper.getDenominator()) {
      fareys.add(Fraction.of(c, d));
      k = (b + maxDem) / d; // integer division
      e = k * c - a;       // next c
      f = k * d - b;       // next d
      a = c;    // update a/b and c/d, and move them to the next pos
      b = d;
      c = e;
      d = f;
    }

    fareys.add(upper);     // last fraction
  } /* makeSeq */


  // ------------------- test rig -------------

  public static void main(String[] args)
  {
    // Create a sequence with max denominator 8 
    // in the default range [0, 1]
    System.out.println("Farey Sequence; den=8; [0 - 1]:");
    FareySeq fs = new FareySeq(8); 
    System.out.println(fs.getFareySeq() + "\n");

    // Create a sequence with max denominator 50 
    // in the range [1/5, 1/2]
    System.out.println("Farey Sequence; den=50; [1/5 - 1/2]:");
    fs = new FareySeq(50, 1, 5, 1, 2);
    // fs.printSeq();
    
    // more compact print-out
    for (int i = 0; i < fs.getSize(); i++) { 
      long[] f = fs.getFrac(i);
      System.out.print(f[0] + "/" + f[1] + " ");
    }
    System.out.println("\n");
  }  // end of main()


}  // end of FareySeq class