\\ ---- PARI/GP program for A352739 -----
\\  Peter Munn, April 2022 

\\ ----------Abbreviations used----------
\\  AFG: see comment at point of use
\\  No = Number
\\  PnD = Primitive nonDeficient (see A006039 for definition)
\\  RFactor = Replacement Factor
\\  Seq = Sequence

A352739uptovalue(UpperLimit) =
{ GetPrime = 1; GetExponent = 2;  \\ Fixed values, giving names to the factor list columns.

  my (SeqIndex);
  forfactored (FactoredCandidate = 1, UpperLimit,
     my (FactorList, PrimeDivisor, Prime_sExponent, ComplementaryDivisor, maybeStronglyPnD,
         AFG, CriticalReplacementPrime, ReplacementPrime, RFactorNo, 
         Candidate_sAbundancy, ComparisonAbundancy);

     Candidate_sAbundancy = sigma(FactoredCandidate)/FactoredCandidate[1];
     if (Candidate_sAbundancy >= 2,
        maybeStronglyPnD = 1;  \\ 1 = True; the candidate is still under consideration
        FactorList = FactoredCandidate[2];

        forstep (FactorNo=omega(FactoredCandidate), 1, -1,  \\ Run through all the factor list
           PrimeDivisor = FactorList[FactorNo,GetPrime];
           ComplementaryDivisor = FactoredCandidate[1]/PrimeDivisor;

           if (sigma(ComplementaryDivisor) >= 2*ComplementaryDivisor,

              \\ This proper divisor is nondeficient, so the candidate is not primitive
              maybeStronglyPnD = 0; break

           ,  \\ otherwise: apply the replacement prime part of the stronger test

              Prime_sExponent = FactorList[FactorNo,GetExponent];
              AFG = PrimeDivisor * (PrimeDivisor^Prime_sExponent - 1) / (PrimeDivisor - 1);
              \\ AFG stands for AbundancyFactorGenerator. PrimeDivisor contributes (AFG+1)/AFG to the
              \\ abundancy of the candidate; the rest comes from the abundancy of ComplementaryDivisor.

              CriticalReplacementPrime = nextprime(AFG+1);
              while (FactoredCandidate[1] % CriticalReplacementPrime == 0,
                 CriticalReplacementPrime = nextprime(CriticalReplacementPrime+1)
              );
              \\ We need only test using the CriticalReplacementPrime and the candidate's prime divisors.

              ReplacementPrime = CriticalReplacementPrime; RFactorNo = omega(FactoredCandidate)+1;
              while (1,
                 ComparisonAbundancy
                  = sigma(ComplementaryDivisor*ReplacementPrime) / (ComplementaryDivisor*ReplacementPrime);
                 if (Candidate_sAbundancy > ComparisonAbundancy && ComparisonAbundancy >= 2,
                    maybeStronglyPnD=0; break(2)  \\ candidate fails the replacement prime test
                 );
                 if (RFactorNo == 1, break) ;
                 ReplacementPrime = FactorList[RFactorNo--,GetPrime];
              ) \\ end of while, running through replacement primes

           ) \\ endif

        ); \\ endfor, running through the candidate's factor list (of primes to be replaced)

        \\ if the candidate remains under consideration, it has passed all the tests
        if (maybeStronglyPnD, print(SeqIndex++, " ", FactoredCandidate[1]))
     )
  ) \\ end of forfactored, running through the candidates
}