We want N divisible by only primes of this form, and squarefree, therefore we cannot have any factor twice or more.
There are only 16 primes of the form 4k + 1 < 150, which leaves us with only 2^16 = 65534 numbers to check, but they can get very larger, so we need an efficient way to check them.
I used itertools.combinations to generate all the necessary subsets, then given a subset say (5, 13) this represents the number 5*13 = 65, because these are both primes of the form 4k + 1, they have solutions 5 = 1^2 + 2^2, 13 = 2^2 + 3^2, and therefore by the above method the 2 solutions to 65 = a^2 + b^2 are 1^2 + 8^2 and 4^2 + 7^2.
For the set (5, 13, 17) I would first solve (5, 13) then with these solutions solve the next one using the solutions to 17 = 1^2 + 4^2, etc, etc
I could've built up a much faster and better solution using memoization but, it is already quite fast, ~35 seconds to get the desired answer, so I left it!
Enter a number (yourinput)
Code will output ∑S(N), for all squarefree N only divisible by primes of the form 4k+1 with 4k+1 < yourinput