# 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 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.

• I guess we have to specify how triplets are ordered and be consistent.. E.g. All (x,y,z) or all (z,x,y) and so on – AZTECCO Nov 4 '19 at 17:06
• @AZTECCO Good point, I will add that specification. – 79037662 Nov 4 '19 at 17:08
• I'm not convinced that solutions with z=0 are trivial. – Neil Nov 4 '19 at 21:22
• @Neil I didn't say they were, but they are invalid for the purposes of this question. – 79037662 Nov 4 '19 at 21:37
• @RosLuP Correct. – 79037662 Nov 8 '19 at 14:35

# Julia 1.0, 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!

• You can drop the space after if to save a byte. – Glen O Nov 8 '19 at 1:21

# 05AB1E, 15 10 bytes

-2 bytes by Kevin Cruijssen
-3 bytes by Grimmy

<L3ãʒImÆ¹Ö


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.

• The specification says the order of the numbers in the triple does not matter (as long as it's consistent). Is the í necessary? – 79037662 Nov 4 '19 at 18:38
• Ops! Yes, then it is not needed. Good catch, thanks. – Wisław Nov 4 '19 at 18:43
• 0Q can be _ to save a byte. – Kevin Cruijssen Nov 5 '19 at 12:37
• 13 bytes by outputting as triplets in the order [z,y,x]. – Kevin Cruijssen Nov 5 '19 at 12:46
• 10 bytes – Grimmy Nov 5 '19 at 15:50

# APL (Dyalog Unicode), 19 bytes

{⍸0=⍵|-/¨⍺*⍨⍳3⍴⍵-1}


The triples are returned in the format x z y. 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.


# Jelly, 13 bytes

*Œ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!

### 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


## GAP, 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 ] ]


# Octave, 74 72 bytes

2 bytes saved thanks to @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!

• Do you need the extra pair of parentheses in a((k=~mod(a.^p+b.^p-c.^p,n)))? – 79037662 Nov 4 '19 at 17:14
• @79037662 No, they are not needed. They were a leftover from some (unsuccessful) tests to try to get rid of the cell array. Thanks! – Luis Mendo Nov 4 '19 at 17:22

# 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!

# R, 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!

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\$$.

• nice, I tried a different approach and only got down to 68 – MickyT Nov 4 '19 at 17:58
• @MickyT Nice idea using which with 3d indices! Here is a 65 byte version of yours, which you should post separately as it is substantially different. – Robin Ryder Nov 4 '19 at 19:18

# R, 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!

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.

# Wolfram Language (Mathematica), 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!

# Ruby, 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!

Thanks Chas Brown for -1 byte.

• Save a byte with (a**p+b**p-c**p)%n<1 – Chas Brown Nov 5 '19 at 21:29
• I never knew you could pass multiple arrays to product! Very nice. – IMP1 Nov 6 '19 at 10:09

# Python 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!

Pretty basic.

# 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!

# Jelly, 13 bytes

*+_ƭ/³ḍ
Ṗṗ3çƇ


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


# 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
└∊─────────────────────────────────────────────────────────────┘


# K (oK), 30 23 bytes

-7 bytes thanks to ngn!

{(~x!-/*/y#)#+1+!3#x-1}


Try it online!

• (*/y#)'' -> */y#, – ngn Nov 10 '19 at 9:28
• using func#list as "filter": {(~x!-/*/y#)#+1+!3#x-1} – ngn Nov 10 '19 at 9:32
• @ngn Filter, of course! Thanks! – Galen Ivanov Nov 10 '19 at 9:36

# Zsh, 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!

Creates the list l=('1 ' '2 ' '3 ' ... '(n-1) ' ), then ${=${:-$^l$^l\$^l} gets the Cartesian product and splits on spaces.

# Pyth, 17 bytes

AQf!%-F^RHTG^StG3


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