// C# program for computing terms of OEIS sequence A344422
//
// Jon E. Schoenfield, 19 Jun 2021
//
// (My apologies for any things I've done inelegantly or
//  inefficiently in the program below. I've done very
//  little programming in C# at present.)
//
// The program below will compile using the Microsoft (R)
//  Visual C# Compiler version 4.8.3752.0 for C# 5,
/// which is provided as part of the Microsoft (R)
//  .NET Framework, but only supports language versions
//  up to C# 5, which is no longer the latest version.
//
// The program compiles and runs successfully for me on a
//  PC running the Windows 10 operating system if I open a
//  Command Prompt window, modify the %path% variable by
//  running a batch file containing only the line
//
// path=%path%C:\Windows\Microsoft.NET\Framework\v4.0.30319;
//
//  and navigate to the folder where this source code file
//  is stored and run the compiler by entering the command
//
// csc /reference:system.numerics.dll A344422.cs
//
//  (where "A344422.cs" is, of course, the name of this
//  source code file).
//
// The "/reference" switch is required because the program
//  declares some variables of type BigInteger so that it
//  can handle terms larger than 2^64 (the first such term
//  is A344422(54) = 25728416344361482752).
//
// The program's output to the screen is in OEIS b-file
//  format, so the first few lines are
//
// 1 1
// 2 2
// 3 4
// 4 6
// 5 44
//
// (At certain intervals, to provide information on the
//  progress of the search, the palindrome currently being
//  factored is temporarily shown below the last line of
//  "permanent" output and is overwritten by the next such
//  update or by the next term found).
//  


using System;
using System.Numerics;

public class A344422
{
   public static void Main()
   {
      int nDgtsMax = 20; // find all terms < 10^nDgtsMax

      int n = 0;

      BigInteger[] Pwr10 = new BigInteger[nDgtsMax+1];
      Pwr10[0] = 1;

      for (int j = 1; j <= nDgtsMax; j++)
      {
         Pwr10[j] = 10*Pwr10[j-1];
      }

      int numPrimesGenerate = 1000; // this will be much more than enough

      int[] prime = new int[numPrimesGenerate+1];
      prime[1] = 2; // (for primes, I prefer to start at index 1)

      int numPrimesStored = 1;

      for (int pTest = 3; numPrimesStored < numPrimesGenerate; pTest += 2)
      {
         for (int iChk = 1; iChk <= numPrimesStored; iChk++)
         {
            int pDivTest = prime[iChk];

            if (pTest % pDivTest == 0)
            {
               break;
            }

            if (pDivTest * pDivTest > pTest)
            {
               numPrimesStored++;
               prime[numPrimesStored] = pTest;
               break;
            }
         }
      }

      int tauMaxFound = 0;


      // Find all 1-digit terms

      for (BigInteger k = 1; k < 10; k++)
      {
         BigInteger quot = k;
         int tau = 1;

         for (int primeIndex = 1; primeIndex < numPrimesStored+1; primeIndex++)
         {
            if (quot % prime[primeIndex] == 0)
            {
               int multiplicity = 0;

               while (quot % prime[primeIndex] == 0)
               {
                  quot = quot / prime[primeIndex];
                  multiplicity++;
               }

               tau *= multiplicity + 1;

               if (quot == 1)
               {
                  break;
               }
            }
         }

         if (tau > tauMaxFound)
         {
            tauMaxFound = tau;
            n++;
            Console.WriteLine(n + " " + k);
         }
      }


      int minNumPrimeFactorsReqd = 0;
      int temp = tauMaxFound;

      while (temp > 0)
      {
         temp = temp/2;
         minNumPrimeFactorsReqd++;
      }

      // Every number having more than tauMaxFound divisors
      //  will have at least minNumPrimeFactorsReqd prime
      //  factors (not necessarily distinct)


      for (int nDgts = 2; nDgts <= nDgtsMax; nDgts++)
      {
         int nDgts1stHalf = (nDgts + 1)/2; //  (if nDgts is odd, the middle digit is considered
                                           //   the 1's digit of the "first half")

         BigInteger[] M = new BigInteger[nDgts1stHalf];

         // M[j] = multiplier (for use in iterating through nDgts-digit palindromes)
         //        for (10^j)'s digit of "first half" of palindrome; e.g.,
         //
         // at nDgts = 5:
         //
         //    nDgts1stHalf = 3;
         //          M[0] =   100,
         //          M[1] =  1010,
         //          M[2] = 10001;
         //
         // at nDgts = 6:
         //
         //    numDgts1stHalf = 3;
         //          M[0] =   1100,
         //          M[1] =  10010,
         //          M[2] = 100001.

         if (nDgts % 2 == 0)
         {
            M[0] = 11 * Pwr10[nDgts1stHalf-1];
         }
         else
         {
            M[0] = Pwr10[nDgts1stHalf-1];
         }

         for (int m = 1; m < nDgts1stHalf; m++)
         {
            if (nDgts % 2 == 0)
            {
               M[m] = Pwr10[nDgts1stHalf+m] + Pwr10[nDgts1stHalf-1-m];
            }
            else
            {
               M[m] = Pwr10[nDgts1stHalf+m-1] + Pwr10[nDgts1stHalf-1-m];
            }
         }

         BigInteger k1stHalfMin = Pwr10[nDgts1stHalf-1];
         BigInteger k1stHalfMax = Pwr10[nDgts1stHalf] - 1;

         BigInteger k = Pwr10[nDgts-1] - 10; // initializing k to this value
                                             //  will make 10^(nDgts-1) + 1
                                             //  the first palindrome k tested
                                             //  in the loop below

         for (BigInteger k1stHalf = k1stHalfMin; k1stHalf <= k1stHalfMax; k1stHalf++)
         {
            // update the value of k so that its "1st half" is k1stHalf
            //  and k remains palindromic

            BigInteger t = k1stHalf;
            int m = 0;

            while (t % 10 == 0)
            {
               k -= 9 * M[m];
               t = t / 10;
               m++;
            }

            k += M[m];

            if (m >= 5)
            {
               Console.Write("\r");
               Console.Write("   " + k); // show progress
            }

            BigInteger quot = k; // initialize the quotient remaining to be factored
            int tau = 1; // initialize partial product; when full, will be number of divisors of k
            int minNumPrimeFactorsReqdLeft = minNumPrimeFactorsReqd;

            int maxDivisorToCheck = 1 + (int) Math.Pow((double) quot, 1.0/minNumPrimeFactorsReqdLeft);
            // "maxDivisorToCheck" is quot^(1/minNumPrimeFactorsReqdLeft),
            //  truncated, plus 1 (added out of an abundance of caution, probably unnecessary)

            // The loop below uses trial division by each prime until it either
            //  finds a prime that divides quot or has tested all primes
            //  less than or equal to maxDivisorToCheck:

            for (int primeIndex = 1; primeIndex <= numPrimesStored; primeIndex++)
            {
               if (prime[primeIndex] > maxDivisorToCheck)
               {
                  break; // quot cannot have enough remaining divisors to make k a term of A344422
               }
               if (quot % prime[primeIndex] == 0) // primeIndex-th prime divides quot
               {
                  int multiplicity = 0;

                  while (quot % prime[primeIndex] == 0)
                  {
                     quot = quot / prime[primeIndex];
                     multiplicity++;
                  }

                  tau *= multiplicity + 1;

                  if (quot == 1) // no more prime factors in quot
                  {
                     break;
                  }

                  int pNext = prime[primeIndex+1];
                  if (quot < pNext * pNext) // quot is a prime
                  {
                     tau *= 2;
                     break;
                  }

                  int tauReqdLeft = tauMaxFound / tau;
                  minNumPrimeFactorsReqdLeft = 0;

                  while (tauReqdLeft > 0)
                  {
                     tauReqdLeft = tauReqdLeft / 2;
                     minNumPrimeFactorsReqdLeft++;
                  }
                  maxDivisorToCheck = 1 + (int) Math.Pow((double) quot, 1.0/minNumPrimeFactorsReqdLeft);
               }
            }

            if (tau > tauMaxFound) // new term found
            {
               tauMaxFound = tau;
               n++;
               Console.Write("\r");
               Console.Write("                                        ");
               Console.Write("\r");
               Console.WriteLine(n + " " + k);

               minNumPrimeFactorsReqd = 0;
               temp = tauMaxFound;

               while (temp > 0)
               {
                  temp=temp/2;
                  minNumPrimeFactorsReqd++;
               }
            }
         }
      }

      Console.Write("\r");
      Console.Write("                                        ");
      Console.Write("\r");
      Console.WriteLine("(no more terms through 10^" + nDgtsMax + ")");
   }
}