Build a Primitive Root Diffuser

Question

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!

Answer 1

Mathematica 158

Spaces not needed:

f@n_ := (m = Mod;
        {p, q} = Sort[{y-x, x, y} /. Solve[x*y==n - 1 && x~GCD~y == 1 < x < y, Integers]][[1, 2 ;;]];
        (g[#~m~p + 1, #~m~q + 1] = m[PrimitiveRoot@n^#, n]) & /@ (Range@n - 1);
         g~ Array~{p, q})

The last one looks like this ...

Answer 2

Haskell, 220 213 207 201 194 182 178

f n=[map(x%)[0..j-1]|x<-[0..i-1]]where j=[x|x<-l,x^2>n,(n-1)&x<1,gcd(n`div`x)x<2]!!0;i=n`div`j;0%0=1;a%b=([x|x<-l,all((>1).(&n).(x^))l]!!0*(a-1)&i%((b-1)&j))&n;l=[1..n-2]
(&)=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> 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> 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> 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]

Answer 3

Python 2: 232 217 212 208 201 198 194 192 187 bytes

N=input()
w=n=G=N-1
Z=range(n)
while sorted(G**x%N-1for x in Z)!=Z:G-=1
while w*w/n|n%w|1-all(w%i|n/w%i for i in Z[2:]):w-=1
n/=w
a=[]
for x in Z:a+=[n*[0]];a[x%w][x%n]=G**x%N
print a[:w]

Now it selects the largest root instead of the smallest, so output is different. Example:

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

Answer 4

MATLAB w/ Symbolic Toolbox - 199 198 bytes

Minified:

function C=F(N),M=N-1;y=1:N;p=y-1;t=@(x)eval(mod(sym(x).^p,N));l=M./y;L=fix(l);k=y(l==L.*gcd(y,L)&y>M^.5);k=k(1);l=M/k;C=eye(l,k);g=y(arrayfun(@(x)sum(t(x)<2)<3,y));C(mod(p,k)*l+mod(p,l)+1)=t(g(1));

Expanded:

function C = primroot( N )

M = N-1;
y = 1:N;
p = y-1;
t = @(x) eval( mod( sym(x).^p, N ) );
l = M./y;
L = fix( l );
k = y( l == L.*gcd(y,L) & y > M^.5 );
k = k(1);
l = M/k;
C = eye( l, k );
g = y(arrayfun( @(x) sum( t(x) < 2 ) < 3, y ));
C(mod(p,k)*l+mod(p,l)+1) = t(g(1));

Sample Outputs

primroot(37)

ans =

 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

primroot(41)

ans =

 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

primroot(31)

ans =

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

primroot(113)

ans =

  1    48    44    78    15    42    95    40   112    65    69    35    98    71    18    73
106     3    31    19     8    45    13    59     7   110    82    94   105    68   100    54
 49    92     9    93    57    24    22    39    64    21   104    20    56    89    91    74
109    34    50    27    53    58    72    66     4    79    63    86    60    55    41    47
 28   101   102    37    81    46    61   103    85    12    11    76    32    67    52    10
 30    84    77    80   111    17    25    70    83    29    36    33     2    96    88    43
 16    90    26     5    14   107    51    75    97    23    87   108    99     6    62    38

A noble effort for MATLAB (even if I do say so myself), the language that requires a1 = b(c); a = a1(d); because it doesn't understand a = b(c)(d). :P

A few notes:

• the routine exploits the fact that mod( j^k, N ) for k = 0..N-1 is strictly positive for positive j, and that j is only a primitive root of N if the number 1 appears exactly twice in the sequence
• the routine exploits the fact that a number x >= 2 can only equal floor(x)*gcd(floor(x),y) for positive y if floor(x) == x and gcd(floor(x),y) == 1. Also, the fact that x and y are coprime iff gcd(x,y) == 1.

Other than that, I couldn't find any problem-specific optimizations.

MATLAB's mod function is particularly painful since I use it 3 times, and it takes a minimum of 8 characters per invocation, as opposed to 3 (i.e. x%y) for virtually every other programming language in existence. The arrayfun (which in every other language is either map or -> is also painful.

In light of these handicaps, I also present the following solution:

Fantasy MATLAB - far fewer bytes

C = @(N):
M = N-1;
y = 1:N;
p = y-1;
t = @(x) eval( sym(x).^p % N );
L = fix( l = M./y );
k = y( l == L.*gcd(y,L) & y > M^.5 )(1);
C = ones( l = M/k, k );
C( p%k*l + p%l + 1 ) = t(y( x -> sum( t(x) < 2 ) < 3 )(1));