stergios/temp/calc_i.c

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// briskei ola ta i

typedef __int128 int128_t;

int128_t power(int128_t b, int128_t ex)
{
  int128_t res = (int128_t)1;
  while (ex > 0) {
    // If the exponent is odd, multiply the result by the current base value
    if (ex % 2 == 1) { // Same as (exponent & 1)
      res = (int128_t)res * b;
    }

    // Square the base for the next iteration
    b = (int128_t)b * b;
    // Halve the exponent (integer division)
    ex /= (int128_t)2; // Same as (exponent >>= 1)
  }
  return res;
}

void print_int128(__int128 n)
{
  if (n == 0) {
    putchar('0');
    return;
  }
  if (n < 0) {
    putchar('-');
    n = -n;
  }

  // A buffer to hold the digits (max digits ~40)
  char buf[40];
  int i = 0;

  // Extract digits in reverse order
  while (n > 0) {
    buf[i++] = (char)((n % 10) + '0');
    n /= 10;
  }

  // Print the digits in the correct order
  while (i > 0) {
    putchar(buf[--i]);
  }
}

const int primes[] = {2, 3, 5};
int exponents[] = {0, 0, 0};
const int NUM_PRIMES = 3;

// Global variable to store the number N (calculated in main)
int128_t N_value = 1;
int128_t a = 1;
int128_t b = 1;
// --- Function to Recursively Generate Factor Pairs ---

/**
 * @brief Recursively generates the first factor (n1) and prints the pair (n1, n2)
 * such that n1 * n2 = N. The generation stops when n1 exceeds sqrt(N).
 *
 * @param prime_index The index of the current prime factor being considered.
 * @param n1 The factor built so far.
 */
void generate_factor_pairs(int prime_index, int128_t n1, int f)
{
  // Base Case: If we have considered all unique prime factors,
  // n1 is a complete factor of N.
  if (prime_index == NUM_PRIMES) {
    // Optimization: Stop if n1 exceeds the square root of N.
    // This prevents printing the pair (n2, n1) after (n1, n2) has been printed.
    // Note: For large N, we should compare n1 * n1 > N_value to avoid
    // issues with long long or the sqrt() function if it's less precise.
    int128_t n2 = (int128_t)N_value / (int128_t)n1;
    if ((int128_t)n1 > (int128_t)n2) {
      return;
    }
    if (((int128_t)n1 + (int128_t)n2) % 2 == 0) {
      int128_t i = (int128_t)a + (int128_t)b + (((int128_t)n1 + (int128_t)n2) >> 1);
      int128_t k = (((int128_t)n2 - (int128_t)n1) >> 1);
      if (f) {
        printf("i: ");
        print_int128(i);
        printf(" k: ");
        print_int128(k);
        printf("   (");
        print_int128(n1);
        printf(" , ");
        print_int128(n2);
        printf(")");
        printf("\n");
      }
      else {
        print_int128(i);
        printf("\n");
      }

      if ((int128_t)a + (int128_t)b - (((int128_t)n1 + (int128_t)n2) >> 1) >= 0) {
        int128_t i = (int128_t)a + (int128_t)b - (((int128_t)n1 + (int128_t)n2) >> 1);
        if (f) {
          printf("i: ");
          print_int128(i);
          printf(" k: ");
          print_int128(k);
          printf("   (");
          print_int128(n1);
          printf(" , ");
          print_int128(n2);
          printf(")");
          printf("\n");
        }
        else {
          print_int128(i);
          printf("\n");
        }
      }
    }
    return;
  }

  // Get the current prime and its exponent
  int prime = primes[prime_index];
  int exponent = exponents[prime_index];
  int128_t power_of_prime = 1;

  // Recursive Step: Iterate through all possible powers of the current prime
  for (int i = 0; i <= exponent; i++) {
    // The new factor n1 is the old n1 multiplied by (prime^i)
    int128_t next_n1 = (int128_t)n1 * (int128_t)power_of_prime;

    // Recursive call for the next prime factor
    generate_factor_pairs(prime_index + 1, next_n1, f);

    // Calculate the next power of the current prime: prime^(i+1)
    if (i < exponent) {
      power_of_prime *= (int128_t)prime;
    }
  }
}

int main(int argc, char **argv)
{
  // long long triple[] = {}
  int p;
  sscanf(argv[1], "%d", &p);
  int f;
  sscanf(argv[2], "%d", &f);
  a = (int128_t)power(10, p);
  b = (int128_t)power(6, p);
  N_value = (int128_t)4 * a * b;
  exponents[0] = 2 * p + 2;
  exponents[1] = p;
  exponents[2] = p;
  generate_factor_pairs(0, 1, 0);
}