[Question]
  [
# Fermat's Last Theorem, mod n
It is a well known fact that for all integers \$p>2\$, there exist no integers \$x, y, z>0\$ such that \$x^p+y^p=z^p\$. However, this statement is not true in general if we consider the integers modulo \$n\$.
You will be given \$n\$ and \$p\$, which are two positive integers with \$n>1\$. Your task will be to write a function or program to compute all positive integers \$x, y, z<n\$ such that \$(x^p+y^p)\$ and \$z^p\$ give the same remainder when divided by \$n\$.
# Input
Any reasonable method of input is allowed. E.g. two separate user inputs, ordered pair, two function parameters, etc.
# Output
Any reasonable method of output is valid, it may be produced by a function or output to the screen. The order the triples are listed does not matter. Triples such as `(1, 2, 3)` and `(2, 1, 3)` are considered distinct, and all distinct triples should be listed exactly once. No invalid/trivial triples such as `(0, 0, 0)` should be output.
The numbers \$x, y, z\$ may have any ordering within each of the triples, but that order should be consistent. For example if \$2^p+2^p\$ and \$3^p\$ have the same remainder when divided by \$n\$, you may list this triple as `(2, 2, 3)` or `(3, 2, 2)`.
# Examples
```
n p -> Possible Output
----------------------------------------------------------------
2 3 -> []
3 3 -> [(1,1,2),(2,2,1)]
3 4 -> []
4 3 -> [(1,2,1),(1,3,2),(2,1,1),(2,2,2),(2,3,3),(3,1,2),(3,2,3)]
```
# Scoring
Shortest code in bytes with no standard loopholes wins.
  
[Answer]
# [Julia 1.0](http://julialang.org/), 61 bytes
```
f(n,p,q=1:n-1)=[(x,y,z) for x=q,y=q,z=q if(x^p+y^p)%n==z^p%n]
```
-1 thanks to Glen O
[Try it online!](https://tio.run/##yyrNyUw0rPj/P00jT6dAp9DW0CpP11DTNlqjQqdSp0pTIS2/SKHCtlCnEoirbAsVMtM0KuIKtCvjCjRV82xtq@IKVPNi/zsUZ@SXK6RpGOkYa3LBOMaoHBMExwQo8x8A "Julia 1.0 – Try It Online")
[Answer]
# [APL (Dyalog Unicode)](https://www.dyalog.com/), 19 bytes
```
{⍸0=⍵|-/¨⍺*⍨⍳3⍴⍵-1}
```
The triples are returned in the format `x z y`. [Try it online!](https://tio.run/##SyzI0U2pTMzJT///v/pR7w4D20e9W2t09Q@teNS7S@tRL5DabPyodwtQVNew9n/ao7YJhJVxPeqbqgBUqWCkkKZgDOcZo/BMMORM/gMA "APL (Dyalog Unicode) – Try It Online")
### Explanation
```
{⍸0=⍵|-/¨⍺*⍨⍳3⍴⍵-1} ⍝ n is our right argument, ⍵, and p is our left, ⍺.
             3⍴⍵-1  ⍝ First, we get a triple of (n-1, n-1, n-1).
            ⍳       ⍝ Then we get the Cartesian product of three range[1, n-1).
                    ⍝ As you will see, this gives a list of (x z y) triples.
         ⍺*⍨        ⍝ We take each element of each triple to the power of p.
      -/¨           ⍝ Here we subtract over each triple, using (element - rest)
                    ⍝ This gives us (x^p - (z^p - y^p)) = x^p + y^p - z^p.
    0=              ⍝ Then we check which (x^p + y^p - z^p) equals 0.
   ⍸                ⍝ And use these results, a Boolean array,
                    ⍝ to find get the indices of the correct triples,
                    ⍝ which are our (x z y)s.
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), ~~15~~ 10 bytes
-2 bytes by Kevin Cruijssen  
-3 bytes by Grimmy
```
<L3ãʒImÆ¹Ö
```
[Try it online!](https://tio.run/##yy9OTMpM/f/fxsf48OJTkzxzD7cd2nl42v//JlzGAA "05AB1E – Try It Online")
Returns a triples in the format `[z,y,x]`.
**Explanation**:
```
<L                | Push the interval [1,n)
  3ã              | Push the triple cartesian product [1,n) x [1,n) x [1,n)
                  | E.g. with n = 3, the list begins as follows:
                  | [[1, 1, 2], [1, 2, 1], [2, 1, 3], ...]
    ʒ             | Filter the list:
     Im           | Take the p-th power of each input
                  | E.g. [1, 1, 2] --> [1^p, 1^p, 2^p]
       Æ          | Reduce by subtraction
                  | E.g. [1^p, 1^p, 2^p] -->  1^p - 1^p - 2^p
        ¹Ö        | Check if 1^p - 1^p - 2^p = 0 (mod n), and filter out if not.
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 13 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
```
*ŒH§%³E
Ṗṗ3çƇ
```
A full program accepting command line arguments `n p` which prints a list representation of the results (Do note that the empty list is represented as no output).
**[Try it online!](https://tio.run/##ASMA3P9qZWxsef//KsWSSMKnJcKzRQrhuZbhuZczw6fGh////zT/Mw "Jelly – Try It Online")**
### How?
```
*ŒH§%³E - Link 1: evaluate a triple: [x, y, z], p
*       - exponentiate               [x^p, y^p, z^p]
 ŒH     - split into halves          [[x^p, y^p], [z^p]]
   §    - sums                       [x^p+y^p, z^p]
     ³  - program's 1st argument     n
    %   - modulo                     [(x^p+y^p)%n, (z^p)%n]
      E - all equal?                 (x^p+y^p)%n == (z^p)%n
Ṗṗ3çƇ - Main Link: n, p
Ṗ     - pop (implicit range of n)    [1,2,3,...,n-1]
  3   - literal three                3
 ṗ    - Cartesian power              [[1,1,1],[1,1,2],...,[n-1,n-1,n-1]]
    Ƈ - filter keep those for which:
   ç  -   last Link (1) as a dyad    Link 1([[1,1,1],[1,1,2],...,[n-1,n-1,n-1]], p)
      - implicit print
```
[Answer]
## [GAP](http://gap-system.org), 63 bytes
```
Filtered(Tuples([1..n-1],3),i->(i[1]^p+i[2]^p-i[3]^p) mod n=0);
```
For example:
```
gap> n:=4; p:=3;
4
3
gap> Filtered(Tuples([1..n-1],3),i->(i[1]^p+i[2]^p-i[3]^p) mod n=0);
[ [ 1, 2, 1 ], [ 1, 3, 2 ], [ 2, 1, 1 ], [ 2, 2, 2 ], [ 2, 3, 3 ], 
  [ 3, 1, 2 ], [ 3, 2, 3 ] ]
```
[Answer]
# [Octave](https://www.gnu.org/software/octave/), ~~74~~ 72 bytes
*2 bytes saved thanks to [@79037662](https://codegolf.stackexchange.com/users/89930/79037662)*
```
@(n,p){[c,b,a]=ndgrid(1:n-1) [a(k=~mod(a.^p+b.^p-c.^p,n)) b(k) c(k)]}{2}
```
[Try it online!](https://tio.run/##y08uSSxL/Z9mq6en999BI0@nQLM6OlknSScx1jYvJb0oM0XD0CpP11BTITpRI9u2Ljc/RSNRL65AOwlI6CYDCZ08TU2FJI1sTYVkIBFbW21U@z9Nw0jHWJMrTcMYRpmAKBMg7z8A)
[Answer]
# JavaScript (ES7), 104 bytes
Takes input as `(n)(p)`.
```
n=>p=>[...Array(n**3)].map((_,i)=>[i/n/n|0,i/n%n|0,i%n]).filter(([x,y,z])=>x*y*z&&(x**p+y**p)%n==z**p%n)
```
[Try it online!](https://tio.run/##dYpBCoMwFET3vUfk/zSNpbqN4DlESrBaUuxPiFIS6d3T4K6UbuYNb@ahX3oZvHHriextTJNKpBqnmk5K2XqvIxDnFfbyqR3AVRjMmymppPdZZLKdjHqUk5nX0QN0QUSx9fkYeORbUUDg3B1jDmSk1JYLI0yDpcXOo5ztHSa4IFSIh29Z/ZH1j6z3Z/oA "JavaScript (Node.js) – Try It Online")
[Answer]
# [R](https://www.r-project.org/), 67 bytes
```
function(n,p,m=t(expand.grid(a<-2:n-1,a,a)))m[,!c(1,1,-1)%*%m^p%%n]
```
[Try it online!](https://tio.run/##NcexCoAgEADQvb9oEO7iDEynqC@JAqkMBy8Rg/7eamh6vFTcWNzFa/YnA1OkMGbY72h5a4/kN7CD7HqWiixZRAwT1SsoUiQVikaEJQrBc3HQkcbKgf4xH@ZdeQA "R – Try It Online")
`expand.grid` creates all possible triplets of numbers. Annoyingly, the output of `expand.grid` is a data frame and not a matrix, so I transpose it to coerce to a matrix `m`. Take the matrix product of `[1,1,-1]` with `m^p` to get \$x^p+y^p-z^p\$, then keep only the rows for which we get \$ 0 \mod n\$.
[Answer]
# [R](https://www.r-project.org/), 65 bytes
3 btye's saved from my original idea by @robin-ryder by using subtraction, applying modulo the result.
```
function(n,p,q=(2:n-1)^p)which(!outer(outer(q,q,'+'),q,'-')%%n,T)
```
[Try it online!](https://tio.run/##K/qfZvs/rTQvuSQzP08jT6dAp9BWw8gqT9dQM65AszwjMzlDQzG/tCS1SANCFuoU6qhrq2uCKF11TVXVPJ0Qzf9pGkY6xppcaRrGMMoERJmAeP8B "R – Try It Online")
Makes use of `outer` to build a 3d array of the addition, then the subtraction. If the remainder of n is 0 output the array indices using `which`.
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 54 bytes
```
n_~f~p_:=Select[Range[n-1]~Tuples~3,n|{1,1,-1}.#^p&]
```
Named function taking the two inputs, like `f[4,3]`.
The `|` symbol typed above needs to be replaced by the 3-byte character representing the "divides" relation; `\[Divides]` works as well, and is used in the TIO example.
Works by testing every possible triple `{x,y,z}`. The only detail of note is that the dot product `{1,1,-1}.#^p&`, when applied to `{x,y,z}`, calculates `x^p+y^p-z^p`.
[Try it online!](https://tio.run/##y00syUjNTSzJTE78/z8vvi6triDeyjY4NSc1uSQ6KDEvPTU6T9cwti6ktCAntbjOWCcvJtolsywzJbU4ttpQx1BH17BWTzmuQC32f0BRZl5JtFJatJGOcayCrYKSDoQZywWXMUbIGKPLmCBkTJBlTBB6QMzY/wA "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# [Ruby](https://www.ruby-lang.org/), ~~67~~ 66 bytes
```
->n,p{w=*1...n;w.product(w,w).select{|a,b,c|(a**p+b**p-c**p)%n<1}}
```
[Try it online!](https://tio.run/##KypNqvyfZvtf1y5Pp6C63FbLUE9PL8@6XK@gKD@lNLlEo1ynXFOvODUnNbmkuiZRJ0knuUYjUUurQDsJSOgmAwlN1Twbw9ra/wUKadFGOsaxXCCGMYJhAmGYAEX@AwA "Ruby – Try It Online")
Thanks Chas Brown for -1 byte.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 13 bytes
```
*+_ƭ/³ḍ
Ṗṗ3çƇ
```
[Try it online!](https://tio.run/##ASQA2/9qZWxsef//Kitfxq0vwrPhuI0K4bmW4bmXM8Onxof///80/zM "Jelly – Try It Online")
-2 bytes thanks to caird coinheringaahing  
-1 byte borrowing a clever trick from Jonathan Allan
# Explanation
```
*+_ƭ/³ḍ  Helper Link - given (x, y, z), determine if x^p + y^p = z^p (mod n)
*        (x^p, y^p, z^p)
    /    Reduce over
 +_ƭ     Tied operator: cycle between + and - (x^p + y^p - z^p)
     ³ḍ  Divisibility check by n (if x^p + y^p - z^p is divisible by n then x^p + y^p is congruent to z^p mod n)
Ṗṗ3çƇ    Main Link
Ṗ        Pop (remove last element) - on an integer, implicit range (1 to n-1)
 ṗ3      Cartesian power ^ 3
   çƇ    Filter: keep triples that are valid by the helper link
```
[Answer]
# [Python 2](https://docs.python.org/2/), 94 bytes
```
n,p=input()
R=range(1,n)
print[(x,y,z)for x in R for y in R for z in R if(x**p+y**p-z**p)%n<1]
```
[Try it online!](https://tio.run/##RY49D4IwFEVn3q/oYihYBz4mIyOuJsTNOCAWaTSvzUuJLX@@Fh0c7slJ7nCv8XbSWIZB32WTpmlAYRqFZrY8g66hHh@SFwIzMKTQXrgTXizZqIk5ppB1bFX/1@WnauQuz83WR@yWiGyDh@Ia4gS8J/WS7Eyz3ENiyUcm0smBrScg6iCNZe3p2BJpWtsbyf4JoRQVVN/UUIvqAw "Python 2 – Try It Online")
Pretty basic.
[Answer]
# [Icon](https://github.com/gtownsend/icon), 95 bytes
```
procedure f(n,p)
x:=1to(t:=n-1)&y:=1to t&z:=1to t&(x^p+y^p-z^p)%n=0&write(x," ",y," ",z)&\q
end
```
[Try it online!](https://tio.run/##y0zOz/v/v6AoPzk1pbQoVSFNI0@nQJOrwsrWsCRfo8TKNk/XUFOtEsxVKFGrgjE0KuIKtCvjCnSr4go0VfNsDdTKizJLUjUqdJQUlHQqwWSVplpMIVdqXgqS@bmJmXkamlwKQADRkKZhpGOsiSyipKurq4SqxJgYJSaElJiATQG5BwA "Icon – Try It Online")
[Answer]
# APL(NARS), chars 129, bytes 258
```
r←f w;n;p;a;v;i;j;b
(n p)←w⋄a←n∣p*⍨v←⍳n-1⋄r←⍬
:for i :in v⋄:for j :in v
:if a∊⍨b←n∣a[i]+a[j]⋄r←r,⊂i,j,a⍳b⋄:endif
:endfor⋄:endfor
```
test:
```
  ⎕fmt f 2 3
┌0─┐
│ 0│
└~─┘
  ⎕fmt f 3 3
┌2────────────────┐
│┌3─────┐ ┌3─────┐│
││ 1 1 2│ │ 2 2 1││
│└~─────┘ └~─────┘2
└∊────────────────┘
  ⎕fmt f 3 4
┌0─┐
│ 0│
└~─┘
  ⎕fmt f 4 3
┌7─────────────────────────────────────────────────────────────┐
│┌3─────┐ ┌3─────┐ ┌3─────┐ ┌3─────┐ ┌3─────┐ ┌3─────┐ ┌3─────┐│
││ 1 2 1│ │ 1 3 2│ │ 2 1 1│ │ 2 2 2│ │ 2 3 3│ │ 3 1 2│ │ 3 2 3││
│└~─────┘ └~─────┘ └~─────┘ └~─────┘ └~─────┘ └~─────┘ └~─────┘2
└∊─────────────────────────────────────────────────────────────┘
```
[Answer]
# [K (oK)](https://github.com/JohnEarnest/ok), ~~30~~ 23 bytes
-7 bytes thanks to ngn!
```
{(~x!-/*/y#)#+1+!3#x-1}
```
[Try it online!](https://tio.run/##y9bNz/7/P82qWqOuQlFXX0u/UllTWdtQW9FYuULXsPZ/WrSJtXHsfwA "K (oK) – Try It Online")
[Answer]
# [Zsh](https://www.zsh.org/), 90 bytes
```
l=({1..$[$1-1]}\ )
for x y z (${=${:-$^l$^l$^l}})(((x**$2+y**$2-z**$2)%$1))||<<<"$x $y $z"
```
[Try it online!](https://tio.run/##LczLCoMwFIThfZ5isKd4jiVC1FXRJ@kFijRYEATbhYn67GlCO/x8y/HvIVgWDmPHqylLupDR5rZfIcpOMxY4eDCtHa1nTffx174LMy9FQdXJJbVPypGMyLa1bZvRAnIgn4VnP0zIu7gcwAHptx8f8@vjlChlUaGO1n@baIMmfAE "Zsh – Try It Online")
Creates the list `l=('1 ' '2 ' '3 ' ... '(n-1) ' )`, then `${=${:-$^l$^l$^l}` gets the Cartesian product and splits on spaces.
[Answer]
# [Pyth](https://github.com/isaacg1/pyth), 17 bytes
```
AQf!%-F^RHTG^StG3
```
[Try it online!](https://tio.run/##K6gsyfj/3zEwTVFV1y0uyCPEPS64xN34/38THWMA "Pyth – Try It Online")
Displays a list of triplets, each as \$(z, x, y)\$.
```
AQ                 # Q = eval(input()), G = Q[0] (=n), H = Q[1] (=p)
            ^StG3  # Get the following: range(1, n) ^ 3 -> List of (z,x,y) tuples
  f                # filter that on the following condition: lambda T:
       ^RHT        # map(lambda x: x^H, T)  --> (z^p, x^p, y^p)
     -F            # reduce by subtraction -> z^p - x^p - y^p = z^p - (x^p + y^p)
    %      G       # modulo n  --> (z^p - (x^p + y^p)) mod n
   !               # not (check if == 0) --> z^p == x^p + y^p mod n
```
]