Partiamo con il dire che è un problema molto complesso quello che cerchi di risolvere, per nulla semplice. Per controllare che \(N\) sia scrivibile come somma di due quadrati (attenzione non primi necessariamente) basta controllare che nella decomposizione in fattori primi di \(N\) non ci siano fattori \(p^k \) dove \(p\) è un numero primo della forma \(4n+3\) e \(k\) è dispari
https://en.wikipedia.org/wiki/Sum_of_two_squares_theorem. Questo è il Teorema della somma dei due quadrati che generalizza il Teorema di Fermat sulla somma dei due quadrati
https://en.wikipedia.org/wiki/Fermat%27s_theorem_on_sums_of_two_squares. Già le tecniche per dimostrare questo fatto utilizzano metodi di matematica relativamente avanzata. Un metodo sicuramente funzionante ma inefficiente per trovare una decomposizione del numero \(N\) come somma di due quadrati, nota che richiede tempi esponenziali, è controllare se \( \sqrt{N - n^2} \) è un intero con \( n \leq \sqrt{N}\). Se invece richiedi che entrambi i numeri siano pure primi beh allora un metodo ancora più inefficiente è crivellare tutte le soluzioni trovate con i primi. So che esistono algoritmi più efficienti con tempistiche polinomiali (per la semplice somma di due quadrati) ma temo di non conoscerli sinceramente.
Se supponiamo che \(N\) sia un primo della forma \(4M+1\) allora per trovare un coppia di interi \( a,b\) non necessariamente primi, tale per cui \(N^2=a^2+b^2\) (e lo so che questo non è il problema che hai chiesto, ma nota che questo è un problema molto più semplice di quello che hai chiesto te), ci sono dei modi che richiedono già la conoscenza di concetti avanzati in Teoria algebrica dei numeri come ad esempio la Jacobsthal sums, per ulteriori dettagli puoi vedere questo link
https://math.stackexchange.com/questions/435139/given-the-norm-of-a-gaussian-integer-how-to-find-the-original-gaussian-integer in cui c'è un link al Paper originale in tedesco. Ma nota che è matematica molto complessa (che esce sicuramente da questa sezione del forum) e richiede molta conoscenza per essere in grado di capire cosa c'è scritto (anche sapendo il tedesco) e in generale per cercare di risolvere questo problema.
Per quanto riguarda la tua domanda originale, ovvero dato \(N \) un intero positivo scrivibile come somma di quadrati, trovare due primi (se esistono) \(p,q\) tale per cui \(N = p^2 + q^2 \), non sono a conoscenza di nessun metodo purtroppo che risolve questo problema, ma in ogni caso sono molto sicuro che se dovesse esistere non puoi dedurlo semplicemente dalla decomposizione in fattori primi di \(N\), ma anzi chiederebbe dei metodi avanzati di Teoria Algebrica dei numeri. Un idea che non ho provato ad esplorare sarebbe quella di partire dalla soluzione di Jacobsthal e potenziarla in altri modi. Ad esempio riconducendoti in qualche modo al caso in cui \(N\) è un primo della forma \(4M+1\). Ma non credo comunque che la Jacobsthal sums dia necessariamente due interi primi se esistono, per cui dovresti ancora farci un lavoro sopra e non mi è molto chiaro come.
Sono stato molto vago perché non sapendo il tuo livello di matematica non sapevo bene come prendere questa domanda, anche perché non sono in grado di darti molti più dettagli a riguardo.