# 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 Nov 4, 2019 at 17:06
• @AZTECCO Good point, I will add that specification. Nov 4, 2019 at 17:08
• I'm not convinced that solutions with z=0 are trivial.
– Neil
Nov 4, 2019 at 21:22
• @Neil I didn't say they were, but they are invalid for the purposes of this question. Nov 4, 2019 at 21:37
• @RosLuP Correct. Nov 8, 2019 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. Nov 8, 2019 at 1:21

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


# 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? Nov 4, 2019 at 18:38
• Ops! Yes, then it is not needed. Good catch, thanks. Nov 4, 2019 at 18:43
• 0Q can be _ to save a byte. Nov 5, 2019 at 12:37
• 13 bytes by outputting as triplets in the order [z,y,x]. Nov 5, 2019 at 12:46
• 10 bytes Nov 5, 2019 at 15:50

# 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:
- 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)))? Nov 4, 2019 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! Nov 4, 2019 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 Nov 4, 2019 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. Nov 4, 2019 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 Nov 5, 2019 at 21:29
• I never knew you could pass multiple arrays to product! Very nice.
– IMP1
Nov 6, 2019 at 10:09

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

• 14 bytes Nov 7, 2019 at 8:11
• @cairdcoinheringaahing nice, thanks! Nov 8, 2019 at 13:15

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

# 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, 2019 at 9:28
• using func#list as "filter": {(~x!-/*/y#)#+1+!3#x-1}
– ngn
Nov 10, 2019 at 9:32
• @ngn Filter, of course! Thanks! Nov 10, 2019 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