Build a Primitive Root Diffuser Introduction

When a room has bare, parallel walls, it can create unpleasant repeating acoustic reflections (echoes).  A diffuser is a device mounted on a wall which creates a blocky surface of many different heights.  This has the effect of breaking up the reflected sound into many smaller reflections (each with a different time delay).  The acoustic quality of the room is thought to be greatly enhanced.

Primitive Root Diffuser

The primitive root diffuser uses a grid of (typically wooden) posts, each with a different height (to obtain a different reflection delay time).  The heights of the posts are chosen according to successive powers of a primitive root G, modulo N (a prime number).

Here are some pictures of a primitive root diffuser.

Task

You must write a program or function which takes a prime number N > 2 as its only input.  It must then select a primitive root of N.  G is a primitive root of N if the successive powers of G modulo N generate all the integers from 1 to N-1.  A more formal definition of primitive roots can be found in the Wikipedia article.

Your code must then place these successive powers into an X-by-Y array (or matrix).  Because we want to hang this on the wall, the diffuser should be a convenient shape - as close to square as possible.  The number of entries in the array will be N-1, so its dimensions X and Y should be chosen such that X*Y = N-1, and X is near but not equal to Y (specifically, X and Y are coprime).  For convenience of display, the horizontal dimension should be greater than the vertical dimension.

The order of placement of successive powers of G (modulo N) should start at one corner of the array, then traverse diagonally through the array, wrapping around in both directions until all array positions are filled.  Then the function should return that array, or print it to stdout, whichever you prefer.

We won't concern ourselves with the measurement units used to actually cut the posts to length; we'll assume the units are appropriate for the range of audio frequencies we want to diffuse, and stick to integers.

There's an online calculator, in case you want to check your work.  However, the site appears to be down currently.

Example

Suppose we choose N = 41.  Our program should select one of the primitive roots of N, which are [6, 7, 11, 12, 13, 15, 17, 19, 22, 24, 26, 28, 29, 30, 34, 35].  Suppose our program selects the smallest primitive root, G = 6.

Our program then chooses the squarest dimensions X,Y for the array such that X*Y = N-1.  The clear choice is X=8 and Y=5.

Then our program populates the array with successive powers of G modulo N (starting with G^0), traversing diagonally.  The program then returns/prints the array:

1   14  32  38  40  27  9   3
18  6   2   28  23  35  39  13
37  26  36  12  4   15  5   29
10  17  33  11  31  24  8   30
16  19  20  34  25  22  21  7


Note how it begins with 6^0 mod 41 = 1, then 6^1 mod 41 = 6, then 6^2 mod 41 = 36, then 6^3 mod 41 = 11, etc.

To clarify, the sequence of powers of 6 modulo 41 is:

[1, 6, 36, 11, 25, 27, 39, 29, 10, 19, 32, 28, 4, 24, 21, 3, 18, 26, 33, 34, 40, 35, 5, 30, 16, 14, 2, 12, 31, 22, 9, 13, 37, 17, 20, 38, 23, 15, 8, 7]


Another Example

N = 31

Assuming the program selects G = 11, X = 6, Y = 5, the output is:

1   26  25  30  5   6
4   11  7   27  20  24
16  13  28  15  18  3
2   21  19  29  10  12
8   22  14  23  9   17


A Third Example

N = 37

Assuming the program selects G = 2, X = 9, Y = 4, the output is:

1   12  33  26  16  7   10  9   34
31  2   24  29  15  32  14  20  18
36  25  4   11  21  30  27  28  3
6   35  13  8   22  5   23  17  19


Scoring

Standard code golf rules - the fewest bytes of code producing correct output wins!
Haskell, 220 213 207 201 194 182 178
f n=[map(x%)[0..j-1]|x&lt;-[0..i-1]]where j=[x|x&lt;-l,x^2&gt;n,(n-1)&amp;x&lt;1,gcd(n`div`x)x&lt;2]!!0;i=n`div`j;0%0=1;a%b=([x|x&lt;-l,all((&gt;1).(&amp;n).(x^))l]!!0*(a-1)&amp;i%((b-1)&amp;j))&amp;n;l=[1..n-2]
(&amp;)=mod

this is a function which returns a matrix (a list of lists of Int).
it works by calculating the smallest possible g, the closet coprime j and k and defines a recursive function that given the indices of a slot in the grid returns the Integer in it, and maps this function over all slots.
example output:
*Main&gt; mapM_ print $ f 13 -- used for printing the matrix line-by-lie
[1,5,12,8]
[3,2,10,11]
[9,6,4,7]

*Main&gt; mapM_ print $ f 43
[1,21,11,16,35,4,41]
[37,3,20,33,5,19,12]
[36,25,9,17,13,15,14]
[42,22,32,27,8,39,2]
[6,40,23,10,38,24,31]
[7,18,34,26,30,28,29]

*Main&gt; mapM_ print $ f 37
[1,12,33,26,16,7,10,9,34]
[31,2,24,29,15,32,14,20,18]
[36,25,4,11,21,30,27,28,3]
[6,35,13,8,22,5,23,17,19]