Project Euler 233 on Hackerrank

Question

I've been struggling with this for over a month now. The coding challenge on HR is similar to the official Project Euler 233, but instead of finding $f(N)=420$, your code may need to solve for $f(N)=12$. Also, instead of the result being the sum of the $n$'s, it is just the count of the $n$'s.

Following is the official challenge:

Let $f(n)$ be the number of points with integer coordinates that are on a circle passing through $(0,0)$, $(0,n)$, $(n,0)$ and $(n,n)$.

It can be shown that $f(10000)=36$.

Given two integers $N$ and $m$, what is the number of all positive integers $n \le N$ such that $f(n)=4m$?

So, I must be going about this the wrong way, and that is basically what the HR community told me, but I really do not see any other way to generate the correct solution.

To use one of the smallest examples I can think of, if $N=13$... then we know that:

$f(5) = 12$, $f(10) = 12$, and $f(13) = 12$, so the answer would be $3$.

I can see that $f(n) = 12$ when the number has one and only one prime of the form $4k+1$.

So my approach would generally involve producing all primes $4k+1 \le N$ in order to determine the answer. But the problem with that approach is that when $f(n)=12$, $N$ can be as large as $5 \times 10^{10}$. That is far too many primes to generate for the code to run in under a couple seconds.

I can only imagine that there must be some way to come up with the answer $3$ if you know that $N=13$ without generating all primes $\le 13$, but I am really clueless.

I guess my actual question would be that: "Is there a way to find the solution without knowing all primes $\le N$ or would you, in fact, need all primes $\le N$ to solve for $f(n)=12$?"

If the answer to that question is "Yes, there is a way to find the solution without knowing all primes $\le N$", then my next question would be "can you offer some clues and/or point me to online material I can research that could help me recognize this solution?"

Thanks in advance!

Answer 1

I found a PDF, "COUNTING PRIMES IN RESIDUE CLASSES" , which shows what appears to be a fast way to count the number of primes congruent to $l$ modulo $k$ up to $x$.

So I guess I'm saying this is kind of an answer (but I still have questions)

The equations are quite complex and I have not yet translated them into layman programmer terms in order to produce a C (my language of choice) algorithm.

I think what I really need is a person who speaks both "math" and "code" to translate this math for me, and then I can implement the algorithm and use it to solve the overall challenge.

I think I might understand $P_2(x,y,l)$ but I am unsure about a few elements in the following sections:

From - 2. Proof of Theorem 1

We now explain the method we used to compute $π(x, k, l)$ for large values of $x$. From now on, we assume that $k$ is fixed and we write $π(x, l)$ instead of $π(x, k, l)$. Let $y$ be a real positive number and let $T(x, y, l)$ be the set of positive integers $n$ such that \begin{cases} n ≤ x,\\ n ≡ l \text{(mod k)},\\ p | n ⇒ p > y \end{cases} Assume that $y$ is such that $x^{1/3} ≤ y ≤ x^{1/2}$. Then each element $n$ of $T(x, y, l)$ has at most two (not necessarily distinct) prime factors. Thus we can split this set into three disjoint subsets $T_0(x, y, l)$, $T_1(x, y, l)$, and $T_2(x, y, l)$, according to the number of (not necessarily distinct) prime factors.

Then the function $P_2(x,y,l)$ is later reduced to:

From - 2.1 Computation of $P_2(x,y,l)$

$$P_2(x,y,l)=\sum_{y<p\le x^{1/2}}\pi(x/p,lp^{-1})-\sum_{y<p\le x^{1/2}}\pi(p-1,lp^{-1})$$

Following are some of what I'm confused about:

I don't understand the significance of $y$.

Could the summation iteration not just say $x^{1/3} < p ≤ x^{1/2}$?

I am not sure what the term $lp^{-1}$ does exactly.

I think that any number $n^{-1}$ is the same as $\frac1n$ which would make it $\frac{l}{p}$ but that is a fraction and $pi(x,l)$ and $P_2(x,l)$ expect $l$ to be a whole number (I assume), so I get pretty lost there.

My best code translation of these summation formulas (which is very incomplete atm) is:

/****

I assume k = 4.

I'm not defining pi(x,l) yet, as it involves more terms than just P2(x,l).

pi(x,l) calls P2(x,l) and P2(x,l) calls pi(x,l) so I declare them both
and then I attempt to define P2(x,l) with my best interpretation but I get stuck.

****/

typedef uint64_t u64;

// The maximum number I'll be working with, 10^11, 100 billion.
#define MAX_N 100000000000

u64 *primes // all primes up to sqrt(MAX_N)
  , *p_last // the last prime up to sqrt(MAX_N)
  , *p_end // pointer that is 1 greater than p_last, so p_curr should always be < p_end
  ;

void init(); // this initializes primes up to sqrt(MAX_N) via generic sieve.

u64 P2(u64 x,u64 l); // declaration needed by pi(x,l)

// place-holder for other terms needed by pi(x,l)

u64 pi(u64 x, u64 l); // declaration needed by P2(x,l)

u64 P2(u64 x, u64 l) { // definition
  u64 sum = 0
    , cbrtx = cbrt(x)
    , sqrtx = sqrt(x)
    , *p_curr = primes
  ;
  while(*p_curr <= cbrtx) ++p_curr; // get to the first prime > cbrtx

  /*
    Note: I am not sure what lp^-1 is in this case.
          I am also not positive that you are supposed to start at cbrtx but it says
          to assume y >= to cbrtx and it says that p > y, so I think it is right.
  */
  while(*p_curr <= sqrtx) sum += pi(x/*p_curr, ???) - pi(*p_curr - 1, ???);

  return sum;

}

/*
  And then pi(x,l) would call one of 2 other functions, pi_e or pi_b:
    if(x > *p_last) {
      return pi_e(x,l), use the equations (of course I have to correctly interpret
        P2 and all the other terms in the equations in order to do this)
    } else {
      pi_b(x,l), use brute iteration through primes already sieved
    }
*/

int main() {
  init();
  // should print the number of primes congruent to 1 mod 4 <= 10^11
  printf("%ld\n", pi(10^11,1)); 
}