Kaelygon

I really hoped Matt Parker would have shown how to construct these divisibility rules, so I came up with my own method. Find prime P and natural number N such that P*N = a *10^b + c
The smaller a and c are, the better. b determines where you split the number.

Example: P=313 N=16 a=5 b=3 c=8

313*16=5008  
313*16=5 *10^3 + 8  

Now we can test for 223795=313*13*11*5
Split 223795 after the b:th number, 3rd in this case.
223795/10^3=223.795 decimal point separates the components A, B

Multiply the A=223 by c and subtract from the rest B=795 multiplied by a

B*a-A*c=  
795*5-223*8=2191  

Repeat if needed till you get to small enough number.

2191/10^3=2.191  
191*5-2*8=939  

which is easy to see that's 3*313
Some bad combinations don't reduce the starting number but they are at least always divisible by P. Those cases could be called Parker divisibility rules.

You can also see that when c is negative, A*c is added rather than subtracted, which explains why some method add or subtract like Vsauce vs James.
This is just a funny trick to simulate division and modulus by 10^b to get smaller number, while preserving the congruence.

It's possible I made mistake somewhere, but was able to get correct answers with other few examples.

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