I have been exploring this particular prime 13238717 which is sum of two squares and has Pythagorean triple.
I found this interesting property and so far I haven’t found any texts about what I wrote below.

This is just my conjecture, I have no formal proof and I have only tested few small primes.
I haven’t found any counter examples yet, but I have checked only few dozen primes and couple composites by hand.

I modified ChatGPT script which lists numbers that have both forms P^2=a^2+b^2 and P=c^2+d^2 and it appears to generate the exact same sequence as: A004431
5 10 13 17 20 25 26 29 34 37 40 41 45 50...

Numbers P seem to always have both following forms Hypotenuse numbers (Pythagorean triples) A009003: P^2=a^2+b^2
Numbers that are the sum of 2 squares A001481: P=c^2+d^2

Wolfram notes, “-- one side of every Pythagorean triple is divisible by 3, another by 4, and another by 5.”
I noticed if P is prime and have both forms: one of the Pythagorean sides (a or b) whichever is divisible by 4 has the exact factors that construct both of the square sum components ‘c’ and ‘d’, with the exception of extra factor ‘2’.

Here’s the conjecture put out more formally based solely on my observations:
Pythagorean sides: a,b and square sum componentsc,d are natural numbers and n#, m# are prime factors.

a^2 + b^2 = P^2 (pyth. triple)  
c^2 + d^2 = P (sum of two squares) 

(a or b) mod 4 = 0
(a or b) factors are = 2* (n1*n2*n3...) * (m1*m2*m3...)
c = n1*n2*n3...
d = m1*m2*m3...
(a or b) = 2*(c*d)

Here’s couple examples:
primes 2 and 5 are trivial exceptions as 1 isn’t a prime factor.

1^2+2^2=5  
1^2+1^2=2

Prime: 13 (first non-trivial prime case)

5^2 + 12^2 = 13^2 (pyth. triple) 
2^2 + 3^2 = 13 (sum of two squares)  

12 factors are 2 2 3
12 mod 4 = 0
c=2
d=3
12 = 2* (2*3)

Prime: 821

429^2 + 700^2 = 821^2 (pyth. triple)  
14^2 + 25^2 = 821 (sum of two squares)  
  
700 mod 4 = 0
700 factors are 2 (5 5) (2 7)
c = 5*5 = 25
d = 2*7 = 14
700 = 2* (25*14)

prime: 13238717

1315508^2 + 13173195^2 = 13238717^2 (pyth. triple)
181^2 + 3634^2 = 13238717 (sum of two squares)  
  
1315508 mod 4 = 0
1315508 factors are 2 (2 23 79) (181)
c=181
d=2*23*79=3634
1315508 = 2* (181*3634)

Some composites have multiple ways to write sum of two squares, which each have different (a or b) counterpart, but not necessarily divisible by 4. composite: 260

(four valid pythagorean side pairs)  
132^2+224^2 = 64^2+252^2 = 100^2+240^2 = 156^2+208^2 = 260^2 
(two valid square sums)
8^2+14^2 = 2^2+16^2 = 260

8^2+14^2:
224 mod 4 = 0
224 factors 2 (2 2 2) (2 7)
c= 2*2*2 = 8
d= 2*7 = 14
224/(8*14) = 2
  
2^2+16^2:
64 mod 4 = 0
64 factors 2 (2) (2 2 2 2)
c = 2
d = 2^4 = 16
64 = 2* (2*16)

composite: 58

40^2+42^2=58^2 (pyth. triple)  
3^2+7^2=58  (sum of two squares) 

42 mod 4 = 2 
Not 0 mod 4 congruent, unlike primes. 
Might be result of both c, d being odd.
42 factors 2 3 7
c=3
d=7
58 = 2* (3*7)

Generally it seems that there’s always at least one Pythagorean component where (a or b) = c*d*2, but I haven’t quite figured why this is the case.

I reckon it has something to do with the fact that mod 4 congruence of 4k+1 doesn’t change when you square it: (4k+1)^2 = 8k*(2k+1)+1
Additionally the fact that when sum of two squares is prime or odd, exactly one of the components is always odd, which may explain why (a or b) isn’t always divisible by 4 with composites.

ChatGPT wrote a counterexample finder I checked up to 100 000.

To my knowledge there isn’t straight up equation that would spit out a Pythagorean triple or sum of square solutions for any integer.
There might be some other way to prove or disprove that a or b = c*d*2 when P is prime, but so far I couldn’t think of any.
It might be something obvious that I am missing, or it’s simply all about congruence rules of additions and multiplications.

Again, this is just what I’ve found from my few tests and I don’t have any formal proof. I couldn’t find any papers or posts specifically about this. This is nothing too important, but I found it interesting enough to share.
I just write bad python code out of interest in number theory without an university degree.

---

Thanks for the solution goes to: @0v0
Here’s summary how I understood this:
Brahmagupta–Fibonacci identity

P=c²+d²
Squaring 'P' results in:
P²=(c²-d²)² + (2cd)²

These 'P²' sum components are equivalent to the Pythagorean legs 'a' and 'b':
a=c²-d²
b=2*c*d

Hence, '(a or b)' always contains the factors '2' and all factors of 'c' and 'd'. 

Additionally, the 2*c*d divisibility by 4 is result of P being odd.
2*c*d being divisible by 4 is true for any odd number as one of the sum of square components ‘c’ or ‘d’ must be even.

c=2*k
b = 2*c*d = 4*k*d