/* A7896.c - Implementation of A001055 and A007896
 20240320 copied from the original in math3000.cxx
 */
#include <stdio.h>        // defines size_t
#define MAXFACTORS 27
long f1055[MAXFACTORS]; /* factors for A001055, denoms for A007896 */
long n7896[MAXFACTORS]; /* numerators for A007896 */
long c1055, c7896; /* counts, i.e. terms of A001055 and A007896 */

/* Find the GCD (greatest common divisor) of a and b by the Euler method */
long euler_gcd(long a, long b)
{
  if (a*b == 0) { return 0; }
  if (a < 0) { a = -a; }
  if (b < 0) { b = -b; }
  while (a != b) {
    if (a > b) {
      a = a-b;
    } else {
      b = b-a;
    }
  }
  return a;
}

/* Given a sequence representing a factorisation, generate the combinations
   of simple fractions whose denominators are the terms in the sequence.
   When there are two or more of the same denominator, the subsequence of
   numerators must be non-decreasing. */
void gen_7896(long seqlen, long offset, long num_min)
{
  long numer, denom, i;
  denom = f1055[offset];
  for (numer = num_min; numer < denom; numer++) {
    if (euler_gcd(numer, denom) == 1) {
      n7896[offset] = numer;
      /* recurse or print */
      if (seqlen > offset+1) {
        long next_mn = (f1055[offset+1] == denom) ? numer : 1;
        gen_7896(seqlen, offset+1, next_mn);
      } else {
        c7896++;
      }
    }
  }
} /* End of gen4.7896 */

/* Compute ("generate") A001055 by counting factorisations. The iteration and
   recursion forces enumeration of only the lexicographically first of each
   permutation for a given set of factors. */
void gen2_1055(long n, long minf, long offset)
{
  long f, r, i;
  if (n == 1) {
    /* We have a solution, and offset is the number of factors */
    /* Generate solutions for A007896 */
    gen_7896(offset, 0, 1);
    c1055++;
    return;
  }
  for (f=minf; f<=n; f++) {
    if ((n % f) == 0) {
      f1055[offset] = f;
      r = n/f;
      gen2_1055(r, f, offset+1);
    }
  }
} /* End of gen2.7896 */

void gen1_1055(long n)
{
  c1055 = c7896 = 0; /* Global counts used in recursive functions */
  if (n < 1) { return; }
  if (n == 1) { c1055 = c7896 = 1; return; }
  gen2_1055(n, 2, 0);
} /* End of gen1.7896 */

/* A001055 enumerates unordered factorisations of n (A001055),
   then for each one generates A007896 via subroutines */
void A001055(void)
{
  long n;
  printf("A007896: ");
  for (n=1; n<=66; n++) {
    gen1_1055(n);
    printf("%ld, ", c7896);
  }
  printf("...\n");
} /* End of A001055() */

int main(int argc, char** argv)
{
  A001055();
  return 0;
}