// Nut.java
// Andrew Davison, ad@coe.psu.ac.th, April 2023

/* Utility methods for Number Theory operations.

   The methods are divided into the following broad categories:

     * special numbers
     * parsing
     * IO
     * series
     * gcd ops
     * primes
     * prime factors
     * continued fractions
     * vectorgram

  The functions rely upon:

     * Java Math:
          https://docs.oracle.com/en/java/javase/16/docs/api/java.base/java/lang/Math.html

     * Java BigInteger:
          https://docs.oracle.com/en/java/javase/16/docs/api/java.base/java/math/BigInteger.html

     * Java BigDecimal:
          https://docs.oracle.com/en/java/javase/16/docs/api/java.base/java/math/BigDecimal.html

     * Parts of Apache Commons Numbers:
          https://commons.apache.org/proper/commons-numbers/commons-numbers-docs/apidocs/index.html

     * java-math-library (JML): 
          https://github.com/TilmanNeumann/java-math-library

     * The BigDecimal class in Big-Math:
          https://eobermuhlner.github.io/big-math/javadoc/v2.3.2/
*/


import java.io.*;
import java.math.*;
import java.util.*;
import java.nio.file.*;
import java.util.function.*;

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

import de.tilman_neumann.jml.transcendental.Ln;
import de.tilman_neumann.jml.transcendental.Exp;
import de.tilman_neumann.jml.precision.Precision;


public class Nut
{
  public static final double EPS = 1e-5;


  // -------------- special numbers ---------------------

  public static BigDecimal PI =  // 100 dp
     new BigDecimal("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679");

  public static BigDecimal PHI =  // 100 dp
     new BigDecimal("1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374");

  public static BigDecimal E =  // 100 dp
     new BigDecimal("2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274");

  public static BigDecimal SQRT2 =  // 100 dp
     new BigDecimal("1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727");



  public static int[] toBase(long number, int base)
  // convert the number to the specified base
  {
    long basePower = base;
    int exponent = 1;
    while (number >= basePower) {
      basePower *= base;
      exponent++;
    }

    int[] rep = new int[exponent];
    for (int i = 0; i < exponent; i++) {
      basePower /= base;
      rep[i] = (int) (number / basePower);
      number %= basePower;
    }
    return rep;
  }  // end of toBase()


  public static BigInteger catalan(int n)
  // https://en.wikipedia.org/wiki/Catalan_number
  {
    BigInteger[] cats = new BigInteger[n + 1];
    cats[0] = new BigInteger("1");
    for (int i = 1; i <= n; i++)
      cats[i] = cats[i - 1].
             multiply(BigInteger.valueOf(4 * i - 2)).
             divide(BigInteger.valueOf(i + 1));
    return cats[n];
  }  // end of catalan()



  public static boolean isPerfect(long num)
  /* A perfect number is a positive integer that is equal 
    to the sum of its proper positive divisors, that is, the 
    sum of its positive divisors excluding the number itself.

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

    Equivalently, a perfect number is a number that is half the 
    sum of all of its positive divisors. 
     e.g. ( 1 + 2 + 3 + 6 )/2 = 6
  */
  {
    if (num < 0)
      return false;

    long temp = 0;
    for(long i=1; i <= num/2; i++)
      if (num%i == 0)
        temp += i;
    return (temp == num);
  }  // end of isPerfect()



  public static boolean isNarcissistic(long num)
  {  return isArmstrong(num);  }

  public static boolean isArmstrong(long num)
  /*  An Armstrong number is an n-digit number that is equal to 
      the sum of the nth powers of its digits. 
      For example, 371 is an Armstrong number, since
          3^3 + 7^3 + 1^3 = 371. 
      It is also known as a narcissistic number.
      https://en.wikipedia.org/wiki/Narcissistic_number
      https://rosettacode.org/wiki/Narcissistic_decimal_number
  */
  {
    if (num < 0)
      return false;

    long numCopy = num;
    int len = String.valueOf(num).length();
    double sum = 0;

    while (numCopy != 0) {
      int lastDigit = (int) numCopy % 10;
      sum += (long) Math.pow(lastDigit, len);
      numCopy = numCopy/10;
    }
    return (num == sum);
  }  // end of isArmstrong()


  public static boolean isTriangular(int num)
  // https://en.wikipedia.org/wiki/Triangular_number
  {
    if (num < 0)
      return false;
    
    int sum = 0;
    for (int i = 1; sum <= num; i++) {
      sum = sum + i;
      if (sum == num)
        return true;
    }
    return false;
  }  // end of isTriangular()


  public static boolean isPentagonal(int num)
  // https://en.wikipedia.org/wiki/Pentagonal_number
  {
    if (num < 0)
      return false;

    int i = 1;
    int pent;
    do {
      pent = (3*i*i - i)/2;
      i += 1;
    }
    while (pent < num);
    return (pent == num);
  }  // end of isPentagonal()


  public static boolean areAmicable(int n1, int n2)
  // https://en.wikipedia.org/wiki/Amicable_numbers
  {
    if ((n1 < 0) || (n2 < 0))
      return false;

    return ((n2 == sumDivisors(n1)) &&
            (n1 == sumDivisors(n2)) );
  }


  private static int sumDivisors(int num)
  // sum all the divisors of num beside itself
  {
    int div = num-1;
    int sum = 0;
    for (int i=div; i > 0; i--)
      if (num % i == 0)
        sum += i;
    return sum;
  } /* sumDivisors */


  public static boolean isPerfectSquare(int n)
  {
    int sqrt = (int) Math.sqrt(n);
    return (sqrt * sqrt == n);
  }


  public static double log(double base, double x) 
  {   return Math.log(x)/Math.log(base); }


  public static BigDecimal log(double base, BigDecimal x, int numDigits) 
  {
    Precision pre = Precision.valueOf(numDigits);
    return Ln.ln(x, pre).divide(
            Ln.ln(BigDecimal.valueOf(base), pre),
                  RoundingMode.HALF_UP);
  }


  public static long pow(long x, long y)
  {
    long res = 1; 
    while (y > 0) {
      // If y is odd, multiply x with the rest
      if ((y & 1) != 0)
        res = res * x;
      // n must be even now
      y = y >> 1;   // y = y/2
      x = x * x;
    }
    return res;
  }

  /* Note: for more versions of pow(), see Java's Math,
     BigInteger, and BigDecimal, and Apache Commons 
     Numbers' ArithmeticUtils class.
  */

  public static BigDecimal pow(BigDecimal a, BigDecimal b)
  {  return pow(a, b, 20);  }


  public static BigDecimal pow(BigDecimal a, BigDecimal b, int numDigits)
  // use identity a^b == e^(b * log a)
  {
    Precision pre = Precision.valueOf(numDigits);
    return Exp.exp(b.multiply( Ln.ln(a, pre)), pre);
  }  // pow()



  public static BigInteger stirling(int n)
  // factorial using stirling approximation
  // https://en.wikipedia.org/wiki/Stirling%27s_approximation
  // n! ~ sqrt(2*pi*n) * pow((n/e), n) 
  {
    if (n == 1)
      return BigInteger.ONE;
     
    MathContext mc = new MathContext(100);
    BigDecimal nBig = BigDecimal.valueOf(n);

    BigDecimal lhs = (BigDecimal.valueOf(2).multiply(Nut.PI.multiply(nBig))).sqrt(mc);
    BigDecimal rhs = (nBig.divide(Nut.E, mc)).pow(n, mc);

    return lhs.multiply(rhs).toBigInteger();   // will discard any fractional part
  } // end of stirling()



  public static double sigmoid(double x)
  // https://en.wikipedia.org/wiki/Sigmoid_function
  {  return 1/(1 + Math.exp(-x));  }

  public static double invSigmoid(double sig)
  { return Math.log( sig/(1 - sig) ); }


  // ---------------- parsing -------------------

  public static int parseInt(String s)
  {
    if (s == null)
      return 0;
    try {
      return Integer.parseInt(s);
    }
    catch (NumberFormatException ex) { 
      System.out.println(s + " could not be parsed as an int; using 0");
      return 0;
    }
  }  // end of parseInt()


  public static long parseLong(String s)
  {
    if (s == null)
      return 0;
    try {
      return Long.parseLong(s);
    }
    catch (NumberFormatException ex) { 
      System.out.println(s + " could not be parsed as a long; using 0");
      return 0;
    }
  }  // end of parseLong()



  public static double parseDouble(String s)
  { 
    if (s == null)
      return 0;
    try {
      return Double.parseDouble(s);
    }
    catch (NumberFormatException ex){ 
      System.out.println(s + " could not be parsed as a double; using 0");
      return 0;
    }
  }  // end of parseDouble()


  // ------------------ IO ------------------------


  public static void pause(int ms)
  {
    try {
      Thread.sleep(ms);
    }
    catch(InterruptedException ex) { }
  }  // end of pause()


  public static void waitEnter() 
  {
    Console c = System.console();
    if (c != null) {
      c.format("Press <enter> to continue...\n");
      c.readLine();
    }
  }


  public static String loadNumberString(String fnm)
  /* The format of the input lines is a single real 
     decimal number separated into digit blocks with spaces
        and newlines
     and blank lines and comment lines.

     The number is returned as a single string, that could be 
     passed as an argument to BigDecimal for instance.
  */
  { 
    File file = new File(fnm);
    if (!file.isFile())
      return null;

    StringBuilder sb = new StringBuilder();
    try {
      BufferedReader br = new BufferedReader( new FileReader(fnm));
      String line;
      boolean isFirst = true;
      while((line = br.readLine()) != null) {
        line = line.trim();
        if (line.length() == 0)  // blank line
          continue;
        if (line.startsWith("//"))   // comment
          continue;

        sb.append( checkLine(line, isFirst));
        isFirst = false;
      }
      br.close();
    } 
    catch (IOException e) 
    {  System.out.println("  " + e);  }

    return sb.toString();
  }  // end of loadNumberString()



  private static String checkLine(String line, boolean isFirst)
  {
    line = line.replaceAll("\\s+","");
    for(int i=0; i < line.length(); i++) {
      char ch = line.charAt(i);
      if (!Character.isDigit(ch)) {
        if (!isFirst)
          System.out.println("Illegal character detected: " + ch);
        else if (isFirst) {
          if (ch != '.')
            System.out.println("Illegal character detected: " + ch);
        }
      }
    }
    return line;
  }  // end of checkLine()



  // ------------------ series ------------------

  public static long fib(int n)
  // https://en.wikipedia.org/wiki/Fibonacci_sequence
  {
    if (n == 0)
      return 0;
    long prev = 0, res = 1, next;
    for (int i = 2; i <= n; i++) {
      next = Math.addExact(prev, res);
      prev = res;
      res = next;
    }
    return res;
  } /* fib */


  public static BigInteger fibBig(int n)
  {
    if (n == 0)
      return BigInteger.ZERO;

    BigInteger prev = BigInteger.valueOf(0);
    BigInteger res = BigInteger.valueOf(1);
    BigInteger next;
    for (int i = 2; i <= n ; i++) {
      next =  prev.add(res);
      prev = res;
      res = next;
    }
    return res;
  }  // end of fibBig()



  public static long lucas(int n)
  /* Calculate nth number of lucas series
     (2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, ....)
     https://en.wikipedia.org/wiki/Lucas_number
  */
  {
    long prev = 2;
    long curr = 1;
    long next;
    for (int i = 1; i <= n; i++) {
      next = Math.addExact(prev, curr);
      prev = curr;
      curr = next;
    }
    return prev;
  } /* lucas */


  public static BigInteger lucasBig(int n)
  {
    BigInteger prev = BigInteger.valueOf(2);
    BigInteger curr = BigInteger.valueOf(1);
    BigInteger next;
    for (int i = 1; i <= n ; i++) {
      next = prev.add(curr);
      prev = curr;
      curr = next;
    }
    return prev;
  }  // end of lucasBig()


  // ------------------- gcd ops -------------------


  public static int gcd(int a, int b)
  // https://en.wikipedia.org/wiki/Greatest_common_divisor
  {
    if ((a < 0) || (b < 0))
      throw new IllegalArgumentException("Negative number");

    return (b == 0) ? a : gcd(b, a % b);
  }


  public static long gcd(long a, long b)
  {
    if ((a < 0) || (b < 0))
      throw new IllegalArgumentException("Negative number");

    while (b != 0) {
      long temp = a % b;
      a = b;
      b = temp;
    }
    return a;
  }


  public static BigInteger gcd(BigInteger a, BigInteger b)
  {  return a.gcd(b);  }


  public static long lcm(long a, long b)
  // https://en.wikipedia.org/wiki/Least_common_multiple
  {  return a * (b / gcd(a, b));  }


  public static BigInteger lcm(BigInteger a, BigInteger b)
  {  return a.multiply(b.divide(a.gcd(b)));  }




  public static long[] egcd(long a, long b)
  /* This function performs the Extended Euclidean algorithm
     to find the GCD of a and b.  This function
     returns three values in an array: {gcd, x, y} such that
             gcd = x*a + y*b
  */
  {
    if ((a < 0) || (b < 0))
      throw new IllegalArgumentException("Negative number");

    long[] vals = new long[3];
    if (b == 0)  { 
      vals[0] = a;
      vals[1] = 1;
      vals[2] = 0;
    }
    else {     //  Otherwise, recurse
      long q = a/b;
      vals = egcd(b, a % b);
      long temp = vals[1] - vals[2]*q;
      vals[1] = vals[2];
      vals[2] = temp;
    }
    return vals;
  }  // end of egcd()
 

  public static long modInv(long a, long n) 
  // Compute b such that ab mod n = 1.
  // Returns n if gcd(a,n) != 1, since b exists
  {
    long[] res = egcd(a, n);   // {gcd, x, y}
    long x = res[1];
    long gcd = res[0];
    if (gcd != 1) 
      return -1;
    x = x % n;
    if (x < 0) 
      x += n;
    return x;
  }  // end of modInv()


  public static ArrayList<Long> modSolve(long a, long b, long n) 
  /* Find all x solutions to ax mod n = b and return them in a list.
     Solutions exist iff gcd(a,n) divides b, and there will
     be gcd(a,n) of them.
  */
  { long gcd = gcd(a,n);
    if( (b % gcd) == 0 ) {
      ArrayList<Long> solns = new ArrayList<>();
        long x = modInv(a/gcd, n/gcd);
        x = (x * b/gcd) % (n/gcd);
        for(long k = 0; k < gcd; k++)
            solns.add((x + k*(n/gcd)) % n);
        return solns;
    }
    return null;
  }  // end of modSolve()



  public static long crt3(long m, long a, long n, long b, 
                          long p, long c) 
  /* Chinese remainder theorem for three congruences.
     Given m, a, n, b, p, c return z which satisfies both
         z mod m = a  
         z mod n = b   and
         z mod p = c

     https://en.wikipedia.org/wiki/Chinese_remainder_theorem
  */
  {
     long z1 = crt(m, a, n, b);
     if (z1 == -1)
       return -1;
     else {
       long z1Mod = lcm(n, m);
       return crt(z1Mod, z1, p, c);
     }
  }  // end of crt3()




  public static long crt(long m, long a, long n, long b) 
  /* Chinese remainder theorem.
     Given m, a, n, b, return z which satisfies
          z mod m = a
     and  z mod n = b.

    z exists iff a == b (mod gcd(m,n)).
    It is unique mod lcm(n,m).
    If z does not exist, then return -1.
  */
  {
    long[] res = egcd(m, n);      // {gcd, s, t}
    long gcd = res[0];
    long s = res[1];
    long t = res[2];

    long l = m/gcd*n;
    long r = a % gcd;

    if ((b % gcd) != r) 
      return -1;

    if (gcd == 1) {
      s = s % l; 
      if (s < 0) 
        s += l;

      t = t % l; 
      if (t < 0) 
        t += l;

      long r1 = (s * b) % l;
      r1 = (r1 * m) % l;

      long r2 = (t * a) % l;
      r2 = (r2 * n) % l;

      return (r1 + r2) % l;
    }
    else
      return gcd*crt(m/gcd, a/gcd, n/gcd, b/gcd) + r;
  }  // end of crt()




  public static long[] linearDiophantine(long a, long b, long c) 
  /* Compute x and y such that that ax + by = c.
     The solution, when it exists, is unique up to the transformation
         x -> x + kb/g
         y -> y - ka/g
     for integers k, where g = gcd(a,b).
     The solution exists iff gcd(a,b) divides c.
     The return value is {x, y} iff the solution exists.
  */
  { long[] res = egcd(a,b);   // {gcd, s, t}
    long gcd = res[0];
    long s = res[1];
    long t = res[2];
    if((c % gcd) != 0) {
      System.out.println("No solution for linear diophantine");
        return null;
    }
    long x = c/gcd * s; 
    long y = c/gcd * t;
    return new long[] {x, y};
  }  // end of linearDiophantine()



  public static int modPow(int a, int i, int m)
  // calculates a^i (mod m)
  // or could use de.tilman_neumann.jml.modular.ModularPower
  {
    BigInteger ai = BigInteger.valueOf(a);
    BigInteger mod = BigInteger.valueOf(m);
    BigInteger mp = ai.modPow(BigInteger.valueOf(i), mod);
    return mp.intValue();
  }


  // ----------------- primes --------------------
  // https://en.wikipedia.org/wiki/Prime_number


  public static boolean isPrime(long n)
  // uses Math.sqrt()
  {
    if ((n <= 0) || (n == 1))
      return false;

    if (n == 2)
      return true;

    if (n % 2 == 0)
      return false;

    int max = (int) Math.sqrt(n) + 1;
    for (int i = 3; i < max; i += 2)
      if (n % i == 0)
        return false;

    return true;
  } /* isPrime */


  public static boolean isPrime(BigInteger n)
  // isPrime() without using sqrt() to compensate for
  // using BigInteger
  {
    if (n.compareTo(BigInteger.ZERO) < 0)
      return false;

    if (n.equals(BigInteger.ZERO) || 
        n.equals(BigInteger.ONE))
      return false;

    BigInteger two = new BigInteger("2");
    BigInteger three = new BigInteger("3");
    if (n.equals(two) || n.equals(three))
      return true;

    if (n.mod(two).equals(BigInteger.ZERO) ||
        n.mod(three).equals(BigInteger.ZERO))
      // n%2 == 0 or n%3 == 0
      return false;

    BigInteger i = new BigInteger("5");
    while(i.multiply(i).compareTo(n) <= 0) {
      if (n.mod(i).equals(BigInteger.ZERO) ||
          n.mod(i.add(two)).equals(BigInteger.ZERO))
        return false;
      i = i.add(BigInteger.valueOf(6));  // i=i+6
    }
    return true;
  } /* isPrime */


  public static boolean isPrime(int num, int[] primes)
  {
    if ((num <= 0) || (num == 1))
      return false;

    for (int i=0; i < primes.length; i++) {
      if (primes[i] == num)
        return true;
      if (primes[i] > num)
        return false;
    }
    System.out.println(num + " out of range (2-" +
              primes[primes.length-1] + ")");
    return false;
  } // end of isPrime()



  public static boolean[] isPrimes(int n)
  /* Returns a boolean array isPrime[i] which indicates
     whether i is prime, for 0 <= i <= n.
     For a large batch of queries, this is faster 
     than calling isPrime() for each integer.
  */
  { if (n < 0)
      throw new IllegalArgumentException("Negative array size");

    boolean[] isPrime = new boolean[n+1];
    if (n >= 2)
      isPrime[2] = true;

    for (int i = 3; i <= n; i += 2)
      isPrime[i] = true;

    // Sieve of Eratosthenes
    int max = (int) Math.sqrt(n) + 1;
    for (int i = 3; i < max; i += 2) {
      if (isPrime[i]) {
        int incr = i*2;
        for (int j = i*i; j <= n; j += incr)
          isPrime[j] = false;
      }
    }
    return isPrime;
  } /* isPrimes */


  public static int[] getPrimes(int n)
  /* Returns all the primes <= n, in ascending
     order. For example: getPrimes(97) = 
     {2, 3, 5, 7, 11, 13, 17, 19, 23, 29,..., 83, 89, 97}.
  */
  {
    boolean[] isPrimes = isPrimes(n);

    ArrayList<Integer> primesList = new ArrayList<>();
    for (int i = 0; i < isPrimes.length; i++)
      if (isPrimes[i])
        primesList.add(i);

    int[] primesArr = new int[primesList.size()];
    for(int i=0; i < primesArr.length; i++)
      primesArr[i] = primesList.get(i);

    return primesArr;
  } /* getPrimes */


  public static int[] getPrimes(int n, int[] primes)
  /* Returns all the primes <= n, in ascending
     order by employing a supplied primes[] array. 
     For example: getPrimes(97, ...) = 
     {2, 3, 5, 7, 11, 13, 17, 19, 23, 29,..., 83, 89, 97}.
  */
  {
    if (n > primes[primes.length-1]) {
      System.out.println("Supplied primes array may be too small");
      return primes;
    }

    ArrayList<Integer> primesList = new ArrayList<>();
    int j = 0;
    while (n >= primes[j]) {
      primesList.add(primes[j]);
      j++;
    }

    int[] primesArr = new int[primesList.size()];
    for(int i=0; i < primesArr.length; i++)
      primesArr[i] = primesList.get(i);

    return primesArr;
  } /* getPrimes */


  /* primes1M.txt contains the first 1 million primes,
     and is loaded by loadPrimes().

     Text files of prime numbers can be obtained from:
       https://www.utm.edu/staff/caldwell/primes/millions/
       https://t5k.org/lists/
       http://www.prime-numbers.org/
  */

  public static int[] loadPrimes(String fnm)
  {
    System.out.println("Loading primes from " + fnm);
    ArrayList<Integer> primes = new ArrayList<>();
    try {
      Scanner scanner = new Scanner(Paths.get(fnm));
      while (scanner.hasNext()) {
        if (scanner.hasNextInt()) {
          primes.add(scanner.nextInt());
        } 
        else
          scanner.next();
      }
      scanner.close();
    } 
    catch (Exception e) 
    {  System.out.println("Unable to load " + fnm);
       System.exit(1);
    } 

    int[] primesArr = new int[primes.size()];
    for(int i = 0; i < primes.size(); i++) 
      primesArr[i] = primes.get(i);

    return primesArr;
  }  // end of loadPrimes()


  public static int[] countsOfPrimes(int n, int[] primes)
  /* Return an array of counts where counts[i] is the
     number of primes <= i, up to a value for counts[n]
  */
  {
    if (n > primes.length-1)  // limit size of n
      n = primes.length-1;

    int[] counts = new int[n+1];  // so can include n

    int primePos = 0;
    int total = 0;
    for(int i = 2; i < n+1; i++) {
      if ((i < primes[primes.length-1]) && (i == primes[primePos])) {
        total++;
        primePos++;
      }
      counts[i] = total;
    }

    return counts;
  }  // end of countPrimes()


  public static BigDecimal sqrt(BigDecimal v)
  {
    return v.sqrt(new MathContext(40));   // precision
  }



  // ------------- prime factors --------------------
  // Use JML or Apache Common Numbers; see Primefactors.java


  public static List<Integer> primeFactors(int num,
                                  int[] primes)
  {
    List<Integer> factors = new ArrayList<Integer>();
    int i = 0;
    int p = primes[i];
               // using p = 2, 3, 5, ...
    while ((num != 1) && (p*p <= num)) {
      while (num % p == 0) {
        num /= p;
        factors.add(p);
      } 
      p = primes[++i]; // only consider primes!
    }
    if (num != 1)
      factors.add(num); 
            // special case if num is actually a prime
    return factors;
  }  // end of primeFactors()

  
  public static int totient(int num, int[] primes)
  /* Euler's totient function
     https://en.wikipedia.org/wiki/Euler%27s_totient_function
     Returns the number of positive integers <= num that are 
     relatively prime to num;
     See BigInteger and gcd() solution in Totient.java.
     This version is using Euler's product formula.

     For example, totient(12) = 4 because these integers 
     are coprime with 12: 1, 5, 7, 11.
  */
  {
    if (num <= 0)
      throw new IllegalArgumentException("Totient of non-positive integer");

    int i = 0;
    int p = primes[i];
    int ans = num;
    while ((num != 1) && (p*p <= num)) {
      if (num % p == 0)
        ans -= ans / p; // only count unique factor
      while (num % p == 0)
        num /= p;
      p = primes[++i];
    }
    if (num != 1)
      ans -= ans / num; // last factor
    return ans;
  }  // end of totient()


  public static BigInteger totient(BigInteger n)
  // also known as Euler's phi function
  // see primes-based soln above
  {
    BigInteger phi = new BigInteger("1");   // phi = 1
    BigInteger i = new BigInteger("2");
    while (i.compareTo(n) < 0) {   // while (i <= n-1)
      if ((i.gcd(n)).equals(BigInteger.ONE))   // if gcd(i,n) == 1
        phi = phi.add(BigInteger.ONE);         // phi++
      i = i.add(BigInteger.ONE);          // i++
    }  
    return phi;
  }  // end of totient()


  public static List<Integer> getDivisors(long num)
  // Note: for BigIntegers use JML's Divisors class
  {
    List<Integer> divisors = new ArrayList<>();
    for(int i=1; i <= num; i++)
       if (num % i == 0)
         divisors.add(i);
    return divisors;
  }  // end of getDivisors()




  // ----------- continued fractions -----------
  // https://en.wikipedia.org/wiki/Continued_fraction


  public static List<Integer> cf(Fraction frac)
  {  return cf(frac.getNumerator(), 
               frac.getDenominator(), 20);  }


  public static List<Integer> cf(Fraction frac, int maxTerms)
  {  return cf(frac.getNumerator(), 
               frac.getDenominator(), maxTerms);  }


  public static List<Integer> cf(double val)
  {  return cf(val, 1, 20);  }


  public static List<Integer> cf(double numer, 
                                 double denom)
  {  return cf(numer, denom, 20);  }


  public static List<Integer> cf(double numer, 
                       double denom, int maxTerms)
  // Return the terms of the continued fraction for 
  // the numerator and denominator as a list
  {
    ArrayList<Integer> terms = new ArrayList<>();
    int quotient = (int)Math.floor(numer/denom);  
                              // integer quotient
    double remainder = numer - quotient*denom;
    int count = 0;

    while ((remainder != 0) && (count < maxTerms)) {
      terms.add(quotient);
      count++;
      double prev_r = remainder;
      quotient = (int)Math.floor(denom/prev_r);
      remainder = denom - quotient*prev_r;
      denom = prev_r;
    }
    terms.add(quotient);
    return terms;
  }  // end of cf()


  public static List<Integer> cf(
               BigDecimal numer, 
               BigDecimal denom, int maxTerms)
  // Return the terms of the continued fraction for 
  // the numerator and denominator as a list
  {
    ArrayList<Integer> terms = new ArrayList<>();

    BigDecimal quot = numer.divideToIntegralValue(denom);
    int quotient = quot.intValue(); // integer quotient
    BigDecimal remainder = 
          numer.subtract(quot.multiply(denom));
    int count = 0;

    while ((remainder.compareTo(BigDecimal.ZERO) != 0) && 
           (count < maxTerms)) {
      terms.add(quotient);
      count++;
      BigDecimal prev_r = remainder;

      quot = denom.divideToIntegralValue(prev_r);
      quotient = quot.intValue();
      remainder = denom.subtract(quot.multiply(prev_r));
      denom = prev_r;
    }
    terms.add(quotient);
    return terms;
  }  // end of cf()


  public static Fraction cf2Frac(List<Integer> seq)
  // Convert the simple continued fraction in `seq` 
  // into a fraction, num / den
  {
    long num = 1;
    long denom = 0;
    for(int i= seq.size()-1; i >= 0; i--) {
      long oldNum = num;
      num = denom + oldNum*seq.get(i);
      denom = oldNum;
    }
    int n = Math.toIntExact(num);   // may get an Arithmetic Exception
    int d = Math.toIntExact(denom); 
    return Fraction.of(n, d);
  }  // end of cf2Frac()


  public static BigFraction cf2BigFrac(List<Integer> seq)
  // Convert the simple continued fraction in `seq` 
  // into a big fraction, num / den
  {
    BigInteger num = BigInteger.ONE;
    BigInteger denom = BigInteger.ZERO;
    for(int i= seq.size()-1; i >= 0; i--) {
      BigInteger oldNum = num;
      num = denom.add( oldNum.multiply(BigInteger.valueOf(seq.get(i))));
      denom = oldNum;
    }
    return BigFraction.of(num, denom);
  }  // end of cf2BigFrac()


  public static List<Fraction> convergents(
                  Fraction convFrac, int maxTerms)
  {
    return convergents(convFrac.getNumerator(), 
                        convFrac.getDenominator(), 
                        maxTerms);
  }


  public static List<Fraction> convergents(
            double val, int maxTerms)
  {  return convergents(val, 1, maxTerms); }


  public static List<Fraction> convergents(
            double numerCon, double denomCon, 
                                  int maxTerms)
  {
    List<Fraction> fracs = new ArrayList<>();

    Fraction prevFrac = Fraction.ZERO;
    for(int i=0; i < maxTerms; i++) {
      Fraction f = Nut.cf2Frac(
                 Nut.cf(numerCon, denomCon, i));
      if (f.equals(prevFrac))
        break;
      fracs.add(f);
      prevFrac = f;
    }
    return fracs;
  }  // end of convergents()


  public static Supplier<Coefficient> 
                        simpleCF(double... coefs) 
  {
    return new Supplier<Coefficient>() 
    {
      private int n;   // iteration count
      private double[] bs = coefs.clone();
                      // denominator terms

      public Coefficient get() 
      {
        if (n != bs.length)
          // Return next term
          return Coefficient.of(1, bs[n++]);
        return Coefficient.of(0, 1);
      }
    };
  } // end of simpleCF()



  public static String getRecurringSeq(int num, int den)
  // https://www.geeksforgeeks.org/find-recurring-sequence-fraction/
  /* Build the decimal sequence for num/den, and return the
     recuring part of the sequence as a 
     string, or "" if there is no such recurrence.
  */
  {
    String seq = "";  // the decimal sequence
    HashMap<Integer, Integer> map = new HashMap<>();
        // remainder --> where it first appears in the sequence 

    int remainder = num % den;
    while((remainder != 0)  &&  (!map.containsKey(remainder))) {
      // Keep calculating remainder until it becomes 0 or repeats
      map.put(remainder, seq.length());
      remainder = remainder * 10;
      int resPart = remainder / den;
      seq = seq + resPart;   // increase sequence
      remainder = remainder % den;
    }

    // System.out.println("seq: " + seq);
    // System.out.println(Arrays.asList(map));

    if (remainder == 0)
      return  "";
    else   // remainder has repeated
      return seq.substring( map.get(remainder) );
          /* extract the sequence from the first appearance
             of the remainder to the end (where the remainder
             repeats) */
  } /* getRecurringSeq */



  // ------------------- vectorgram ------------

  /* Inspired by the vectorgram visualization in:
     "Wonders of Numbers", Clifford A. Pickover, OUP, 2001
     Chapter 55, Audioacfive Decay, p.134-138
  */

  public static void vectorgram(int numDirs, String numStr)
  {  vectorgram(numDirs, numStr, 1);  }


  public static void vectorgram(int numDirs, String numStr,
                                      double userScale)
  {
    double angleIncr = 2*Math.PI/numDirs;
    int len = numStr.length();

    double sf = ((double)len)/(4*numDirs) * userScale;
          // a bit of a hack
    System.out.println("sf: " + sf);
    double scaleFactor = (sf > 25) ? sf : 25;   // more hackery
    System.out.println("scaleFactor: " + scaleFactor);

    StdDraw.setScale(0, scaleFactor);

    // starting position
    double x1 = scaleFactor/2;  
    double y1 = scaleFactor/2;
    double x2, y2;

    StdDraw.setPenColor(StdDraw.BLACK);
    for(char ch: numStr.toCharArray()) {
      if (Character.isDigit(ch)) {
        double angle = ((ch - '0')%numDirs)*angleIncr;
            // the modulo ensures that all the digits are
            // mapped to a direction
        x2 = x1 + Math.cos(angle);
        y2 = y1 + Math.sin(angle);
        StdDraw.line(x1, y1, x2, y2);
        x1 = x2; y1 = y2;
      }
    }
  }  // end of vectorgram()


}  // end of Nut class