[Question]
  [
Write a program that will determine if a given matrix represents a quandle. A [quandle](https://en.wikipedia.org/wiki/Racks_and_quandles) is a set equipped with a single (non-commutative, non-associative) operation \$◃\$ which obeys the following axioms:
* The operation is closed, meaning that \$a◃b = c\$ is always an element of the set if \$a\$ and \$b\$ are elements of the set.
* The operation is right-self-distributive: \$(a◃b)◃c = (a◃c)◃(b◃c)\$.
* The operation is right-divisible: For any chosen pair of \$a\$ and \$b\$, there is a single unique \$c\$ such that \$c◃a = b\$
* The operation is idempotent: \$a◃a = a\$
A finite quandle can be represented as a square matrix. Below is an example of an order-\$5\$ quandle ([source](https://arxiv.org/pdf/math/0412417.pdf)).
```
0 0 1 1 1
1 1 0 0 0
3 4 2 4 3
4 2 4 3 2
2 3 3 2 4
```
The value located at the \$n\$-th row and \$m\$-th column (0-indexed) is the value of \$n◃m\$. For example, in this quandle, \$4◃1 = 3\$. Some of the quandle properties are easy to see from this matrix:
* It is closed because only values \$0-4\$ appear in this \$5\times5\$ matrix.
* It is idempotent because the matrix diagonal is \$0, 1, 2, 3, 4\$
* It is right-divisible because no column contains any duplicate values. (The rows can, and usually will.)
The property of right-self-distributivity is harder to test. There might be a shortcut, but the simplest method is to iterate over each possible combination of three indexes to verify that `m[m[a][b]][c] = m[m[a][c]][m[b][c]]`.
## Input
Input will be the list of rows of a square matrix, using either 0-indexing or 1-index (your choice). Each entry will be a single digit number from `0` to `8` (or `1` through `9`). I will be flexible on the input format. Some acceptable formats include:
* Your language's most natural formatting for matrices or lists, such as `[[0 0 0][2 1 1][1 2 2]]` or `(0,0,0,2,1,1,1,2,2)`.
* The list of values delimited by whitespace, newlines, commas, etc.
* A single string consisting of all values concatenated together, such as `000211122`.
You are also allowed to take the transpose of the matrix as input (swapping rows with columns). Just be sure to state this in your answer.
## Output
A single truthy/falsey value indicating the matrix's status as a quandle.
## Examples of quandles
```
0
0 0
1 1
0 0 0
2 1 1
1 2 2
0 0 1 1
1 1 0 0
3 3 2 2
2 2 3 3
0 3 4 1 2
2 1 0 4 3
3 4 2 0 1
4 2 1 3 0
1 0 3 2 4
```
## Examples of non-quandles
**not closed**
```
1
0 0 0
2 1 1
1 9 2
```
**not right-self-distributive**
```
0 0 1 0
1 1 0 1
2 3 2 2
3 2 3 3
(3◃1)◃2 = 2◃2 = 2
(3◃2)◃(1◃2) = 3◃0 = 3
```
**not right-divisible**
```
0 2 3 4 1
0 1 2 3 4
3 4 2 2 2
3 3 3 3 3
4 1 1 1 4
0 1 2 3
3 1 2 0
3 1 2 3
0 1 2 3
```
**not idempotent**
```
1 1 1 1
3 3 3 3
2 2 2 2
0 0 0 0
2 1 0 4 3
3 4 2 0 1
4 2 1 3 0
1 0 3 2 4
0 3 4 1 2
```
  
[Answer]
## Haskell, 100 bytes
This answer uses **transposed** input.
```
q m=and$(elem<$>v<*>m)++[a#a==a&&a#b#c==a#c#(b#c)|a<-v,b<-v,c<-v]where v=[0..length m-1];i#j=m!!j!!i
```
Seems like I can't use a pattern guard to bind an infix operator, so I'm using `where` in this case. 
(First version was 108 bytes but missed idempotency test, fixed version was 120 bytes, later versions had 108, 103 and 98 bytes but I had to realize thanks to @nimi that they were all wrong: of course I need to do the right-divisibility test (which implies closedness) before doing any dangerous `!!` operations, but I could still use most of my later golfing ideas and with one more, it was 102 bytes, now improved by changing the operand order of `#` (which is nicer anyway to compensate the transposition) to make better use of its associating to the left)
Use like this:
```
*Main> q [[0,1,2,3],[0,1,3,2],[1,0,2,3],[0,1,2,3]]
False
```
[Answer]
# [Python 2](https://docs.python.org/2/), ~~104~~ ~~103~~ 102 bytes
```
t=input();e=enumerate
[0%(a==A[a]in B>C[B[a]]==t[C[b]][C[a]])for(a,A)in e(t)for(b,B)in e(t)for C in t]
```
Input is transposed. Output is via exit code, so **0** (success) is truthy and **1** (failure) is falsy.
[Try it online!](https://tio.run/nexus/python2#jVFNb4MwDD2XX@HLBEg5QMNlnTKprbrrLrshDsBcDa0DlAZt@/XMEKehUzVNihI/fzw/O3X3iqAgDMPRqKbtBxPFD6iwHT5QlwaDPLmLSqW2eVk0Lewe9/mOzEIpk@/zqijoJhgfOx2VYhtTDkZmhpXYLSDsgYApRmoVBJ9vzQnhRQ@4CVZGf9O9wi@soSY9ZPe6aQ2E56Gu8XwOAwrW2Bs4PD8dtO70lF9pLN9dZOOLjmVzGjSGYx4lIi6C6YE0FmBfdqwXPkJL05UIkFcRh8i0BPIX8tSUmnEDOddmszuZHJkrSD2D9bkGlsASZzNrKm7rvr@t22n7S/cikwutgkk6i7DCpB9k7YZLeRB55bjMKZ3uC2XCcT7cgA@zT8kgJA/M3K7wP2gq9Mv3cV6@H3zt9pL5D5F@t37581p/AA "Python 2 – TIO Nexus")
### How it works
`e(t)` returns the enumerated rows of the input matrix **t** – which represents an operator **◃** – as (index,row) pairs. `for(a,A)in e(t)`, e.g., iterates over these, storing the index in **a** and the row itself in **A**, so `A` becomes a shortcut for `t[a]`.
Between the former, `for(b,B)in e(t)`, and `for C in t`, we iterate over all possible ordered tuples **(a, b, c)** in the Cartesian power **t3**.
For each of these tuples, we evaluate the expression
```
0%(a==A[a]in B>C[B[a]]==t[C[b]][C[a]])
```
The value of the parenthesized Boolean is *False* if and only if one or more of the following individual comparisons do the same.
* `a==A[a]` will fail (for some value of **a**) iff **◃** is not idempotent.
* `A[a]in B` will fail if **B** doesn't contain all indices of **A**. 
Since **A** has **n** indices and **B** has **n** elements, non-failure means that the elements of **B** match the indices of **A**, so **◃** is closed and right-divisible.
* `B>C[B[a]]` is a tautology, since Python 2 considered numbers "smaller" than iterables.
* `C[B[a]]==t[C[b]][C[a]]` will fail for some value if **◃** is not right-self-distributive.
If any one of the comparisons returns *False*, the expression `(0%...)` will throw a *ZeroDivisionError*. Also, if **◃** is not closed `A[a]` or `C[b]` may also throw an *IndexError*. In both cases, the program exits with status code **1** (failure).
If all test passed, the program will exit normally with status code **0** (success).
[Answer]
# [Python 2](https://docs.python.org/2/), 138 bytes
```
def f(m):R=range(len(m));return all(m[i][i]==i<set(zip(*m)[i])==set(R)>m[m[j][k]][i]==m[m[j][i]][m[k][i]]for i in R for j in R for k in R)
```
`m` is a list of lists of integers.
[Try it online!](https://tio.run/nexus/python2#jZLPboMwDMbvfgrfSKapgtDL/jBpL1CpvTIOCMKWQlIUqFTt5Ts7RGJdL1NA@fzZJj/ktLpDr@tWWPlMYj57h2Vp3CwuEruTxwsah34zjYOZhayC59mz0Us@XCIrgH1h61Hwtx4xSZIUIMUUMsyCIKmQgwwVqsVawiwkc8xDgl5SORfkuKWkCn0p6RzYUdwIvGdUwQekoXVLZ65EzCQBdn@g7lieVpY0smSgIku@sqiFBrgw6MiylMUFzMtrC/@FXn/zHv/a0mw6Hsyh8LX71GLQjkL5EudUD4OwpanoKQrzOulZfJtRPFhJjiwKNg7yzZa2PFZlXy2FMTQUWjJZ8EwNz/QQxntcZR@kvNbTpP0cTuxubsZewq@cO814m99RwejpPmHyPgzYfOmmn3A6N43WrW43yfUH "Python 2 – TIO Nexus")
[Answer]
## JavaScript (ES6), 150 bytes
```
a=>!(a.some((b,i)=>b[i]-i)|a.some(b=>[...new Set(b)].sort()+''!=[...b.keys()])||a.some((_,i)=>a.some((_,j)=>a.some((b,k)=>b[a[j][i]]-a[b[j]][b[i]]))))
```
Takes input as an array of column arrays of integers.
[Answer]
# Mathematica, 122 bytes
```
(n_±m_:=#[[m,n]];#&@@Union[Sort/@#]==Range@l==Array[#±#&,l=Length@#]&&And@@Flatten@Array[±##2±#==(#2±#)±(#3±#)&,{l,l,l}])&
```
Pure function taking a 2D array of integers (1-indexed) as input, with rows and columns reversed from the convention in the question, and returning `True` or `False`. The first line defines the binary infix operation `n_±m_` to be the quandle operation.
For an `l` x `l` array, closed and right-divisible is equivalent to every row (in my case) being some permutation of `{1, ..., l}`, and idempotent is equivalent to the main diagonal being exactly `{1, ..., l}`. So `#&@@Union[Sort/@#]==Range@l==Array[#±#&,l=Length@#]` detects for these three conditions. (The use of `Sort/@#` here is why I chose to swap rows and columns.)
For right-distributive, we literally check all possibilities using `Array[±##2±#==(#2±#)±(#3±#)&,{l,l,l}])`. (Note that `±##2±#` automatically expands to `(#2±#3)±#`, since `##2` represents the sequence of the second and third arguments to the three-variable pure function being arrayed over.) Then `&&And@@Flatten@` checks whether every single test was passed. For some non-closed quandles, errors might be thrown when it tries to access a part of a matrix that doesn't exist, but the correct answer is still returned.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~38~~ 36 bytes
```
ŒDḢJƑ;ZQƑ€$;Jṗ3Ṗịṛ/ịⱮ¥œịʋ=Ṫṭœị⁺ʋɗ€ƲẠ
```
[Try it online!](https://tio.run/##jVC7SkNBEO3nK7awdx@3kZAuVTobQXsbyQ/YmTQpbpUETMBGFBtRwSAhm4DFDfdC8hezP7KZmd0khY0sy5yZc86e3b277fXuY6zHHVy@dJtR6@ayGYXBx1mri37q0D/iukT/dE4lfH9Vb/WE0K5so39H/yld6K925XZKrmaOq@eI6x9czqm9oh0eZpthPUY/AxL719D/rRbVYjMk7jpGDaCVBqOMAIJWcWOUVTaNUmuEdMoJQZuQY4FTBZFWfJqwA55YNgJXQwoO0GItAP4kXZySdE4yYHOSOyXZlAUsFJyTkiwv4NvwKuAgJJKrztUd55CU5uCWZ/Fpcj8y/PtRx2/YAw "Jelly – Try It Online")
Takes input as a matrix. The Footer converts the test case formats into matrix formats for you.
Truly a monstrosity of a program. This is broken into 3 stages, each of which checks a different axiom (except closedness as this is implied by right-divisibility):
```
ŒDḢJƑ - Checks for idempotentness
ZQƑ€  - Checks for right-divisibility and closedness
Jṗ3Ṗịṛ/ịⱮ¥œịʋ=Ṫṭœị⁺ʋɗ€ - Checks for right-self-distributiveness
```
We group these 3 results together using `;` and check that `Ạ`ll are true.
## How each stage works
We use `œị` to calculate \$b ◃ a\$, when the left argument is `[a, b]` and the right argument is `M`. `œị` is the "multi-dimensional index-into" command, which retrieves the `b`th element of the `a`th row of `M`. As we're checking every triple of coordinates, and the arguments are reversed for both sides of the indexing equality, the fact that we're calculating \$b ◃ a\$ instead of \$a ◃ b\$ doesn't matter. That means that instead of checking that \$(a◃b)◃c=(a◃c)◃(b◃c)\$, we instead check that \$c◃(b◃a) = (b◃c)◃(a◃c)\$ for all triples \$a, b, c\$.
Each stage is run as a "monadic link", meaning that the only argument passed to it is the inputted matrix.
```
ŒDḢJƑ - Argument: an n×n matrix M
ŒD    - Yield the diagonals of M
  Ḣ   - Extract the first i.e. the main diagonal
   JƑ - Is it in the form [1, 2, 3, ..., n]?
ZQƑ€ - Argument: an n×n matrix M
Z    - Transpose M
   € - Over each row:
 QƑ  -   Is it unchanged when deduplicated?
Jṗ3Ṗịṛ/ịⱮ¥œịʋ=Ṫṭœị⁺ʋɗ€ - Argument: an n×n matrix M
J                      - Extract n and yield [1, 2, 3, ..., n]
 ṗ3                    - Yield all triples of elements of [1, 2, 3, ..., n]
                    ɗ€ - Over each triple [a,b,c], do the following, with M as the right argument:
              =        -   Are the two following results equal?
             ʋ         -     First result:
   Ṗ                   -       Yield [a,b]
         ¥             -       Group the previous 2 links into a dyad:
     ṛ/                -         Yield c
        Ɱ              -         Over each row in M:
       ị               -           Get the c'th element
                                  This returns a list [1◃c, 2◃c, 3◃c, ..., n◃c]
    ị                  -       Get the a'th and b'th elements
                                  This yields [a◃c, b◃c]
          œị           -       Calculate (b◃c)◃(a◃c); See above
                     ʋ   -     Second result:
               Ṫ         -       Yield c and replace the left argument with [a, b]
                 œị     -        Calculate (b◃a)
                ṭ        -       Tack on c; [b◃a, c]
                    ⁺    -       Repeat the previous link (œị), yielding c◃(b◃a)
              =          -     Does c◃(b◃a) equal (b◃c)◃(a◃c)?
```
* The first yields either `1` or `0`, depending on whether the matrix is idempotent or not.
* The second yields an array of \$n\$ `1`s or `0`s, for each row in `M`. `M` is only right-divisible and closed iff all of these are `1`
* The third yields a list of \$n^3\$ `1`s if `M` is right-self-distributive, or a list with at least one `0` if not.
Finally, we concatenate all of these results into a single list. If that list contains a `0`, then `M` is not a quandle, otherwise it is.
[Answer]
# C++14, 175 bytes
As unnamed lambda, assuming `n` to be like `std::vector<std::vector<int>>` and returning via reference parameter. 0 is false, everything else is true.
```
#define F(x);for(x=-1;++x<s;){
[](auto m,int&r){int s=r=m.size(),a,b,c F(a)auto A=m[a];r*=s==A.size()&&A[a]==a;int u=0 F(b)u|=1<<m[b][a];r*=A[b]<s F(c)r*=m[A[b]][c]==m[A[c]][m[b][c]];}}r*=!(u-(1<<s)+1);}}
```
Ungolfed and usage:
```
#include<vector>
#include<iostream>
auto f=
#define F(x);for(x=-1;++x<s;){
[](auto m,int&r){
 int s=r=m.size(),a,b,c
 F(a)
  auto A=m[a];               //shortcut for this row
  r*=s==A.size()&&A[a]==a;   //square and idempotet
  int u=0                    //bitset for uniqueness in col
  F(b)
   u|=1<<m[b][a];            //count this item
   r*=A[b]<s                 //closed
   F(c)
    r*=m[A[b]][c]==m[A[c]][m[b][c]];
   }
  }
  r*=!(u-(1<<s)+1);          //check right-divisibility
 }
}
;
int main(){
 int r;
 std::vector<std::vector<int>>
  A = {
   {0, 0, 1, 1},
   {1, 1, 0, 0},
   {3, 3, 2, 2},
   {2, 2, 3, 3},
  },
  B = {
   {0, 2, 3, 4, 1},
   {0, 1, 2, 3, 4},
   {3, 4, 2, 2, 2},
   {3, 3, 3, 3, 3},
   {4, 1, 1, 1, 4},
  };
 f(A,r);
 std::cout << r << "\n";
 f(B,r);
 std::cout << r << "\n";
}
```
]