Magical Modulo Squares

Jelly, 13 10 bytes

-3 thanks to Nick Kennedy

^{Feels like the repeated code should be is golf-able, but I have did not managed it...}

*€Ṗ%µQƑƇḢị

Try it online! (footer pretty-formats as a grid)

How?

*€Ṗ%µQƑƇḢị - Link: integer, p
 €         - for each n in [1..p]
*          -   exponentiate with:
  Ṗ        -     pop = [1..p-1]
           - ...i.e [[1^1,1^2,...,1^(p-1)],[2^1,2^2,...,2^(p-1)],...,[....,p^(p-1)]]
   %       - modulo p
    µ      - start a new monadic chain (call that list of lists X)
       Ƈ   - keep those which:
      Ƒ    -   are invariant under:
     Q     -     de-duplicate
        Ḣ  - head
         ị - index into the list of lists X

Charcoal, 36 bytes

≔Ｅ…¹θ﹪Ｘι…¹θＩθηＥ⊟Φη⁼¹№ι¹⪫Ｅ§η⊖ι◧ＩλＬ⊖θ 

Try it online! Link is to verbose version of code. Note: Trailing space. Explanation:

≔Ｅ…¹θ﹪Ｘι…¹θＩθη

Create a p-1 by p-1 array of powers of 1..p-1 to indices 1..p-1 (modulo p).

Ｅ⊟Φη⁼¹№ι¹

Map over one of the rows which has exactly one 1.

⪫Ｅ§η⊖ι◧ＩλＬ⊖θ 

Rearrange the rows into the order given by the selected row and format the output.

JavaScript (ES7),  91  86 bytes

This version attempts to compute the powers before applying the modulo and will fail for \$p\ge11\$ because of loss of precision. It is otherwise using the same logic as the commented version below.

f=(p,k)=>(g=k=>[...Array(i=p-1)].map(_=>k**++i%p))(k).sort()[1]>1?g(k).map(g):f(p,-~k)

Try it online!

---

JavaScript (ES6),  92  87 bytes

This version uses modular exponentiation to support (much) higher input values.

f=(p,k)=>(g=k=>[...Array(p-1)].map(_=>n=n*k%p,n=1))(k).sort()[1]>1?g(k).map(g):f(p,-~k)

Try it online!

How?

Finding the first row

Given \$1\le k<p\$, we use the helper function \$g\$ to compute \$a_k(n)=k^n\bmod p\$ for \$1\le n<p\$.

g = k =>              // k = input
  [...Array(p - 1)]   // generate an array of size p - 1
  .map(_ =>           // for each entry in there:
    n = n * k % p,    //   update n to (n * k) mod p
    n = 1             //   starting with n = 1
  )                   // end of map()

We look for \$k\$ such that there's only one value \$n\$ such that \$a_k(n)=1\$. We do that by sorting the array and testing if the 2^nd element is greater than \$1\$.

g(k).sort()[1] > 1

This works even in lexicographical order -- which is the default behavior of sort() -- because:

• if there are several \$1\$'s, they will all be moved to the front just like they would in numeric order
• if there's only a single \$1\$, the 2^nd value will be greater than \$1\$, no matter if it's really the 2^nd value in numeric order or not

Example:

For \$p=17\$:

• for \$k=1\$, we get:
  • \$a_1=[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\$
  • sorted as \$[ 1, \color{red}1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]\$
• for \$k=2\$, we get:
  • \$a_2=[ 2, 4, 8, 16, 15, 13, 9, 1, 2, 4, 8, 16, 15, 13, 9, 1 ]\$
  • sorted as \$[ 1, \color{red}1, 13, 13, 15, 15, 16, 16, 2, 2, 4, 4, 8, 8, 9, 9 ]\$
• for \$k=3\$, we get:
  • \$a_3=[ 3, 9, 10, 13, 5, 15, 11, 16, 14, 8, 7, 4, 12, 2, 6, 1 ]\$
  • sorted as \$[ 1, \color{red}{10}, 11, 12, 13, 14, 15, 16, 2, 3, 4, 5, 6, 7, 8, 9 ]\$

Building the matrix

Once we've found \$k\$, we invoke \$g(k)\$ again (to retrieve the unsorted version of the array) and invoke \$g\$ on each element of \$g(k)\$ to build the rows of the matrix.

This part can be simply written as:

g(k).map(g)