I wrote a program (randorder32v4.cpp) which randomly generated a file of numbers satisfying the following conditions;

- each line of the file consisted of 32 random element values
- no two element values in a line are the same
- each of those 32 element values is squared, then the squares are summed, and the sum is a square number, otherwise the element values are rejected as a solution
- none of the randomly generated values is larger than 2505
- the smallest value is 1
- no two solution lines have all the same values
- the GCD of all values in any solution line is 1
- the total number of solution lines is 144161
- an extra 3 numbers were added at the start of each line, the order (32) and the square root of the sum of the squares, twice.

The program randorder32v4 took about an hour to run

Conditions 1/,2/, 3/, 4/, 5/, 8/ and 9/ were produced by the program.
Condition 1/ simulates a simple perfect squared square (SPSS) with 32 elements, ie order 32.
Condition 2/ simulates 'perfection' in simple 'perfect' squared squares, ie no two elements are the same.
Condition 3/ simulates a perfect squared square, dissected into 32 squares, ie the sum of the constituent squares equals the area of the larger square (no attempt was made to try and tile the larger square).
Condition 4/ simulates the actual enumerated order 32 SPSSs, in that no element in order 32 is larger than 2505.
Condition 5/ simulates squared squares in that 0 element solutions can appear (and are visible as a cross in the tiling) but these lower the order and so are not allowed.
Conditions 6/ and 8/ simulate simple perfect squared square order 32 which has exactly 144161 solutions. Condition 6/ was confirmed by using sort -u n in linux , ie duplicate lines eliminated (there were no duplicates).
Condition 7/ simulates squared squares where the GCD of elements is always 1. The file was loaded into a spreadsheet as a CSV file and the GCD function was used to test condition 7/, all lines had a GCD of 1.
Condition 9/ was done to simulate the format of Bouwkampcode (without parentheses) so my prime decomposition program would work.

My bk2prime program processes the 144161 lines and outputs them in a modified form, where the composite numbers only are replaced with 0's and the order, sum of composites and sum of primes for each solution line of elements is given at the end of each line.
The prime decomposition file is then loaded into a spreadsheet and using the COUNTIF function the number of lines with a given number of primes is arrived at.

The prime counts ranged from 0 to 15 and the total frequencies produced for each was 806, 4433, 12609, 21581, 27742, 27619, 21781, 14228, 7696, 3566, 1426, 476, 146, 44, 6, 2. (I copied and pasted the values this time to eliminate typing error).
The frequencies add to 144161.

I also run the program, with a slight change, so that I was summing the 32 values, without squaring each value first. The sum was still checked to see if it was a square. This version of the program ran much faster, only taking a few minutes.
The distribution of primes I arrived at looked much the same , ie the prime counts ranged from 0 to 16 and the total frequencies produced for each was 798, 4429, 12279, 21875, 27874, 27726, 21760, 14157, 7621, 3553, 1421, 498, 131, 30, 7, 1, 1

The distributions look quite similar to the order 32 SPSS distribution. The hypothesis that the prime elements are distributed randomly among the elements has some confirmation and has not been falsified by the 2 program runs.
Tiling squares with smaller squares imposes a topological structure on the elements, in practise this would eliminate most, if not all of the random solutions generated here as actual tilings, so the prime distributions seem to be determined by number theoretical considerations, not topological ones.

The sources used are here; randorder32v4.cpp and bk2prime.cpp

There are 368 primes less than 2505, if I am picking 1 number at random from 1 to 2505, I would have expected the probability it being prime to be 368/2505. If I am picking 32 numbers, a rough estimate would be 32*368/2505 = 4.7 primes on average. The mean and mode of the order 32 SPSS data, and the sum of squares random numbers data is about 6.

How to account for the discrepancy?

What if I relax the condition that the sum of the elements squared (or just the sum of the elements) has to be a square (changing part of condition 3)?
I made this change to the program and generated 144161 solutions in seconds, then decomposed the solutions into primes and composites and counted the prime frequencies like before; 822, 4471, 12253, 21863, 28189, 27568, 21867, 14076, 7549, 3476, 1333, 511, 133, 44, 3, 3. The frequencies range from 0 primes to 15 primes. This time the 2 highest frequencies are 4 and 5, this is consistent with a 4.7 prime average. The requirement for a squared sum appears to raise the average number of primes from approx 5 to approx 6. Why is that so? Perhaps Fermat's theorem on the sum of two squares or the theory of quadratic forms has some relevance?