[Question]
  [
Inspired by and drawns from [Is this number Loeschian?](https://codegolf.stackexchange.com/q/88732/85527)
> 
> A positive integer \$k\$ is a **Loeschian number** if
> 
> 
> * \$k\$ can be expressed as \$i^2 + j^2 + i\times j\$ for \$i\$, \$j\$ integers.
> 
> 
> For example, the first positive Loeschian numbers are: \$1\$ (\$i=1, j=0\$); \$3\$ (\$i=j=1\$); \$4\$ (\$i=2, j=0\$); \$7\$ (\$i=2, j=1\$); \$9\$ (\$i=-3, j=3\$)1; ... Note that \$i, j\$ for a given \$k\$ are not unique. For example, \$9\$ can also be generated with \$i=3, j=0\$.
> 
> 
> Other equivalent characterizations of these numbers are:
> 
> 
> * \$k\$ can be expressed as \$i^2 + j^2 + i\times j\$ for \$i, j\$ non-negative integers. (For each pair of integers \$i, j\$ there's a pair of non-negative integers that gives the same \$k\$)
> * There is a set of \$k\$ contiguous hexagons that forms a tesselation on a hexagonal grid (see illustrations for [\$k = 4\$](https://i.stack.imgur.com/keJqMl.png) and for [\$k = 7\$](https://i.stack.imgur.com/Dj0lKl.png)). (Because of this property, these numbers find application in [mobile cellular communication networks](http://www.wirelesscommunication.nl/reference/chaptr04/cellplan/reuse.htm).)
> * See more characterizations in the [OEIS page](https://oeis.org/A003136) of the sequence.
> 
> 
> 
The first few Loeschian numbers are
```
0, 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 27, 28, 31, 36, 37, 39, 43, 48, 49, 52, 57, 61, 63, 64, 67, 73, 75, 76, 79, 81, 84, 91, 93, 97, 100, 103, 108, 109, 111, 112, 117, 121, 124, 127, 129, 133, 139, 144, 147, 148, 151, 156, 157, 163, 169, 171, 172, 175, 181, 183, 189, 192...
```
1while (\$i=-3, j=3\$) produces 9, stick to non-negative integers, so (\$i=0, j=3\$).
Loeschian numbers also appear in [determining if a coincident point in a pair of rotated hexagonal lattices is closest to the origin?](https://math.stackexchange.com/q/3705547/284619)
## The challenge
Given a **non-negative integer** \$k\$, output all pairs of non-negative integers \$i, j\$ such that \$i^2 + j^2 + i\times j=k\$. If none are found (i.e. \$k\$ is not Loeschian) then return nothing or some suitable flag *other than* \$(0, 0)\$ since that produces the first Loeschian number, \$0\$.
For reversed order pairs like \$(0, 4)\$ and \$(4, 0)\$ either include both, or one member of the pair, but it should be the same for all cases (i.e. not sometimes one and other times both).
The program or function should handle (say in less than a minute) inputs up to \$100,000\$, or up to data type limitations.
This is code golf so shortest code wins.
## Test cases
```
 in       out
 0      (0, 0)
 1      (0, 1), (1, 0)
 3      (1, 1)
 4      (0, 2), (2, 0)
 9      (0, 3), (3, 0)
 12     (2, 2)
 16     (0, 4), (4, 0)
 27     (3, 3)
 49     (0, 7), (3, 5), (5, 3), (7, 0)
 147    (2, 11), (7, 7), (11, 2)
 169    (0, 13), (7, 8), (8, 7), (13, 0)
 196    (0, 14), (6, 10), (10, 6), (14, 0)
 361    (0, 19), (5, 16), (16, 5), (19, 0)
 507    (1, 22), (13, 13), (22, 1)
 2028   (2, 44), (26, 26), (44, 2)
 8281   (0, 91), (11, 85), (19, 80), (39, 65), (49, 56), (56, 49), (65, 39), (80, 19), (85, 11), (91, 0)
 12103  (2, 109), (21, 98), (27, 94), (34, 89), (49, 77), (61, 66), (66, 61), (77, 49), (89, 34), (94, 27), (98, 21), (109, 2)
```
  
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), 11 9 bytes
*Thanks to @ovs for -2!*
```
ÝãʒãÀ¦POQ
```
Oh my, I am getting in the hang of code golfing!
Prints a list of all the valid pairs (e.g `[[1, 0], [0, 1]]`). If there are none, the list is empty (`[]`). Also outputs both of any reverse integer pairs.
[Try it online!](https://tio.run/##ARoA5f9vc2FiaWX//8Odw6PKksOjw4DCplBPUf//OQ)
## How?
You may count this as a port of the other answers, but I took a look only at the Husk answer before writing the program!
```
Ý             # Push a list of all numbers from 0 to the input.
 ã            # Push the cartesian power of lists. (Basically, finding all possible pairs)
  ʒ           # For each pair...
   ãÀ¦        # Find all other permutations of the pair.
      P       # Multiply each permutation.
       O      # Add the products.
        Q     # If the result is not equal to the input, yeet (throw) them from the list.
              # Automatically print the pairs not yeeted.
```
[Answer]
# [Haskell](https://www.haskell.org/), 46 bytes
```
f x|l<-[0..x]=[(i,j)|i<-l,j<-l,i*i+j*j+i*j==x]
```
[Try it online!](https://tio.run/##JY2xDoIwFEV/5cU4AH2QtmiBhPoFOjk2jWHQ0FqQCMYO/HutYTn35iyn7@bn3bkQHuBX1@aKFoXXUiUGbbqaNndo/zCZITazxGRWSq/D0JkRJAzddLnB9Fmuy/s8gkrm/vUFnwIhsAN5iohvs3uIjXT1W4RxTYhiAnmFhwbZoUIm4jYCS8HwSCvklNdY85oh44yWWocf "Haskell – Try It Online")
[Answer]
# [Husk](https://github.com/barbuz/Husk), ~~15~~ 13 bytes
```
fo=¹§+Πṁ□π2…0
```
[Try it online!](https://tio.run/##ASMA3P9odXNr//9mbz3CucKnK86g4bmB4pahz4Ay4oCmMP///zE5Ng "Husk – Try It Online")
-2 bytes from Zgarb.
Outputs `[]` for non-Loeschians.
## Explanation
```
fo=¹§+Πṁ□π2…0
           …0 range from 0..n
         π2   create all possible pairs using 0..n
fo            filter by the following two functions:
    §         f: fork: § f g h x = f (g x) (h x)
     +           add
       ṁ□        sum of squares
      Π          and fold by multiplication
  =¹          g: is that equal to 1?
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~11~~ ~~10~~ 9 bytes
```
Żp`ḋÄ$=¥Ƈ
```
[Try it online!](https://tio.run/##y0rNyan8///o7oKEhzu6D7eo2B5aeqz9/@H2R01rjk56uHPG//8GOgqGOgpGOgrGOgomOgqmQNISAA "Jelly – Try It Online")
Outputs `[]` for non-Loeschian numbers
-1 byte [thanks to](https://codegolf.stackexchange.com/questions/212300/find-all-integer-pairs-that-produce-a-given-loeschian-number/212310?noredirect=1#comment500552_212310) [Sisyphus](https://codegolf.stackexchange.com/users/48931/sisyphus)
Not particularly efficient, but that can be fixed for an additional [2 bytes](https://tio.run/##ASIA3f9qZWxsef//w4bCvcW7cGDCslMrUMayPcKlxof///84Mjgx).
Uses the fact that a Loeschian number can be expressed as \$i\times i + j\times(i+j)\$ by using Jelly’s vectorisation and cumulative sum.
## How it works
```
Żp`ḋÄ$=¥Ƈ - Main link. Takes n on the left
Ż         - Yield [0, 1, ..., n]
 p`       - Cartesian product with itself, yielding [[0, 0], [0, 1], ..., [n, n]]
       ¥Ƈ - Filter the pairs, keeping those where the following is true:
     $=   -   The pair equals n after the following is done:
    Ä     -     Cumulative sum. Yield [i, i+j]
   ḋ      -     Dot product with [i, j]; Yields i×i + j×(i+j)
```
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 46 bytes
```
Solve[i^2+j^2+i*j==#&&i>=j>=0,{i,j},Integers]&
```
[Try it online!](https://tio.run/##DYq9CsIwGEVfRRAy6IXmzyQdUrq6CY6iECTEL9gKNbiUPnvMcDnnwJ1CecUpFHqGmny9ft6/eKOHPOY2OmTv94zR4PPgOVZC3nCeS0xx@d5ZvSw0l103pm5cOQQUNHoICWEgLXRzbVs09gbKCJy4heTSwUkn2lNwtdX6Bw "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# [Scala](http://www.scala-lang.org/), ~~64~~ 62 bytes
```
k=>0.to(k)flatMap(i=>0.to(k)filter(j=>i*i+j*j+i*j==k)map(i->))
```
[Try it online!](https://tio.run/##VcoxC8IwEAXgvb/ixlyrUlERhQQcHTo5ikOsKVwa09geIoi/PbZTyfa@995Qa6djd7emZqg0eTAfNv4xwCkE@GZv7aA5wtkzSAUX87qKMS@mAm8gYytVueJOtNg4zZUOguaGHJteWKkop8LmtqDcStnic7otFWIMPXl2XjSiRMxmrRNtEm0THdIt5a7cj/7FPw "Scala – Try It Online")
Thanks to [user](https://codegolf.stackexchange.com/users/95792/user) for -2
[Answer]
# [Japt](https://github.com/ETHproductions/japt), 13 [bytes](https://en.wikipedia.org/wiki/ISO/IEC_8859-1)
```
ô ï f@¶Xx²+X×
```
[Try it](https://petershaggynoble.github.io/Japt-Interpreter/?v=1.4.6&code=9CDvIGZAtlh4sitY1w&input=MTIKLVE)
```
ô ï f@¶Xx²+X×     :Implicit input of integer U
ô                 :Range [0,U]
  ï               :Cartesian product
    f             :Filter by
     @            :Passing each X through the following function
      ¶           :  Is U equal to
       Xx         :  X reduced by addition
         ²        :  After squaring each
          +X×     :  Plus X reduced by multiplication
```
[Answer]
# [Perl 5](https://www.perl.org/), ~~81~~ ~~77~~ 74 bytes
```
sub{map{$i=$_;grep{$k==$i**2+$i*$_+$_**2&&($_=[$i,$_])}$i..$k}0..($k=pop)}
```
[Try it online!](https://tio.run/##ZVTbbqMwEH3PV1iRlULrRrYBY5e6yg/s275RhHZVEtE2ASWptFGUb8@OB0i4EAkP4zNnzvEE6mL/HV23J0LX9nr4@Xve/qnPtLQ0Tzb7AsIva2n5@Cif4E7zJ5pDvFh4NLcpLRnNM/9Cy@WSfl34cukBvK5q/3JNZj@HgvwuDseXl1/Vvkhm0GN1hGfrzQhcKeFuIfYNwpQzwjMCP9ZuiuGmyBhJxRgU9EDCgfqb4ZBBOgY5ZjBDUOBAwUSLvIOkYxpsqgFD6BjCMYOM76DAtRkINQOGuNUQuTVqNcUTTWHc1yREi8JqIaYqzeA4O1LtVt2VTZ0bNShDdwoCjnhIKQwmhgMlBoWmNSMavGrtCTMujHjcH6mUnbBGspTjMUsude8gQpQooYHEVmE4PgottehJM6I7MX2TpNFeAJHCXAhRhHQREIdoRrnRYKRvBnXUTcKI6Z9I8OCmU3AskAAzOAMJwzCoPQDF2nRtY5yMApxCAQoEqGbYcSdFAy7AWuPcYoWBocrGGjesfwJ@MltXew/fRf@Mqe3Jo@WuZrT4V/t2Be9@myYryNjmQ0DzlGev1q0iu6wc9g7bVEe7oGuk8W9peFpXdu6S5Pn5jcyXn1W58x7IA4PPzDxtOFlDmc2Zo2mry0P@URT198l7d1n27pQwJATEZfZR7YrcWSh3m@T6Hw "Perl 5 – Try It Online")
A bit ungolfed:
```
sub f {
  my $k=pop;                      #gangnam style, k=pop from input
  grep { $k==pop@$_ }             #pop last of three elems
                                  #...in the candidate array
                                  #...and return as result
                                  #...if last = i*i+i*j+j*j = k
  map  {                          #two loops from 0 to sqrt $k
    my $i=$_;                     #outer loop var
    map {
      my $j=$_;                   #inner loop var
      [$i, $j, $i*$i+$i*$j+$j*$j] #result candidate
    }
    0..sqrt$k                     #or  $i..sqrt$k  to return only i<=j
  }
  0..sqrt$k
}
```
Note: Saving bytes by removing the two sqrt makes it run A LOT slower, but it will still return the correct result.
[Answer]
# [Scala](http://www.scala-lang.org/), 53 bytes
```
| =>for(i<-0 to|;j<-0 to|if| ==i*i+j*j+i*j)yield(i,j)
```
[Try it online!](https://tio.run/##VcnNCgIhGIXh/VzFt1RnAqMi@jFo2aJVy2hhMwqfiNqMRNF07WYQhLvznHdopZXJX41qIxwlOlCPqFw3wD4EeFV3aUGv4eAiiB2c1O1M8m6@B72ASGO@te8Jbiccoh835jdQ5ySQYW2YqZEZ@kRlO4KNoSn06KJ1RBNOafXXtNCs0LzQqmwlF3yZ/U4f "Scala – Try It Online")
[Answer]
# [C (gcc)](https://gcc.gnu.org/), 75 bytes
```
i,j;f(x){for(i=j=x;~j;i-=!i?j--,-x:1)i*i+j*j+i*j-x||printf("(%d,%d)",i,j);}
```
[Try it online!](https://tio.run/##NY3LCsIwFET3/YprpJA0iVjEjWnqt0hi9AZtSx8QaOuvx6i4mmE4nDHyZkyMKLxyNLDZtT1F7XVQL69Q6g2evZRChlPJsEDuC8@x8DIsS9djMzpKaG5FbhkRScLUGrfYmMdkr1ANo8V2d6@zBMLzgg1lMGfw/UgLgoYSVMoKjvtUOEdgCYC/Oj9Y0DUkNVNpdvSX3TQOlJBPX7M1vgE "C (gcc) – Try It Online")
[Answer]
## Python 3 - ~~79~~ 70 Bytes
Just the trivial solution to get started with Python 3
~~```
k=9  # not counted
```
```
def f(k):
    a=range(k);return [(i,j) for i in a for j in a if k==i*i+j*j+i*j]
```
```
print(f(k))   # not counted
```
[Try it online!](https://tio.run/##VYy7DoAgDAB3vqIjDzcnTfgS42AiaEtSSIODX4/o5o2Xy5W7npnH1pKf1B4iRJ3MrKAjoV7CsGgcyEDMAgjIIBsfoUefoZ/BCMl7tOjIkkNLqyqCXPU7Na09 "Python 3 – Try It Online")~~
---
**Thanks to @Jakque**
lambda, testing framework in TIO, spaces, `k+1`, factorization
```
lambda k:[(i,j)for i in range(k+1)for j in range(k+1)if(i+j)*i+j*j==k]
```
[Try it online!](https://tio.run/##XVJBcoMwDLznFT6axAdkjGNnxi@hHNJpaIEUMoQc@npqCQw4HLBYaaVd4cff@NN32VS5j@l@/f38urL2UvBaNEnVD6xmdceGa/d94@0JCGpiqK54fWqSo38dG@fachpvz/HJXHFg/ilSwegpuI/SpCzFjMMeh0QwDlE@W/OA@RVXe55Enox4dp/PMJ/Fc6VY8hL5G67FxlPIUxFPnkM@w76bHrvjnZd5OZ75Mv8cz1dzI5oPsBQQEeBNkRXbhkIrg6cJjDdvVu8YZMIDkFKphzQFsa9Mw45jF@Ewl@rFCtiIk6fBAyqWQcmsUcrof8lUGhH8KtIkfVtJA5SKHBtpQAQtFsJOzKrBkJXMR5owXH5OnXLfU5F6jXunyKyOTB52bddbVh7wNreCDbfn6z7iraabeyEt/Wt8vEbmWMXbhJDHUHcjx0/nZopA9lyYHKZ/ "Python 3 – Try It Online")
[Answer]
# [JavaScript (V8)](https://v8.dev/), 57 bytes
Prints the pairs \$(x,y),\:x\le y\$.
```
n=>{for(y=n+1;x=y--;)for(;x--;)x*x+y*y+x*y-n||print(x,y)}
```
[Try it online!](https://tio.run/##TY3BCoMwDIbve4ocdqhaReumlVJvu@4FxkCRuTmkFRVpUZ/dpbcdEr7/S0K@9VJPzdgNc7jwo5WHkuXa6pFYqYJEGGnDUHhOCOPI@Cawvg2Mb0O1bcPYqZkYar39EI@YQkIhpXChUCAzrIwCy9G4fMmdcFSgTjNcvsboWMw4Bc544o6SOH2eInx5q5sPUSBLWKHRatL9K@r1m1R3kHBe1V55AlqisP@PMe7e8QM "JavaScript (V8) – Try It Online")
[Answer]
# [Charcoal](https://github.com/somebody1234/Charcoal), 25 bytes
```
ＮθＦ⊕₂θＦ⊕ι¿⁼θ⁻Ｘ⁺ικ²×ικＩ⟦ικ
```
[Try it online!](https://tio.run/##ZYy7CgIxFER7vyLlDURw13JLsdjCJaidWMSYxWAe5ibRz493wc5iYDhzGP1QqKNyrY3hVctU/c0gJD6s5ogMxqDReBOKucMpVYXmGGOhnXP2J1iCdmawJ9FlSIIdbKgZZPzQp3RUrWBPLlhPOVtvfmB5k2hDgZ3KBS4LvHI@tNb13Wbb1m/3BQ "Charcoal – Try It Online") Link is to verbose version of code. Only outputs those pairs where `i>=j`. `₂` speeds the code up so that the larger test cases complete within a minute, but it is not needed for smaller test cases. Explanation:
```
Ｎθ
```
Input `k`.
```
Ｆ⊕₂θ
```
Loop `i` from `0` to `√k` inclusive.
```
Ｆ⊕ι
```
Loop `j` from `0` to `i` inclusive.
```
¿⁼θ⁻Ｘ⁺ικ²×ικ
```
If `k=(i+j)²-ij`, then...
```
Ｉ⟦ικ
```
Output `i` and `j` on separate lines.
Just for fun, here's a 73-byte [Retina](https://github.com/m-ender/retina/wiki/The-Language) 1.0 answer that only finds nontrivial solutions (i.e. neither `i` nor `j` is zero):
```
.+
*
L$w`^((_)+)(?=(?<-2>\1)+(?(2)$.)(_(_)*)(?<-4>\1\3)*$(?(4).))
$.1 $.3
```
[Try it online!](https://tio.run/##K0otycxLNPz/X0@bS4vLR6U8IU5DI15TW1PD3lbD3kbXyC7GUFNbw17DSFNFT1MjHiinpQmSMAFKxBhraqkA5Uw09TQ1uVT0DBVU9Iz//zc1MAcA "Retina – Try It Online") Very slow, so don't try anything over about 500.
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 42 bytes
```
Array[(+##)^2-##&,{#,#}+1,0]~Position~#-1&
```
[Try it online!](https://tio.run/##FcvLCoMwEIXh/bzGgLQ4gVysl0VL@wZCl2IhSKVZqJBmU8S8ejrZzf9xZrHh815scJNN8zU9vLe/4VQinl9aIBa0I@FRKpJj7LevC25bIwpVpN67NQwobvMdxyI@J7vGHSSBItAEhqAi6Dg5VM3WsOSumgz56phNzfuLZNNStwStblV@UtLAkf4 "Wolfram Language (Mathematica) – Try It Online")
Gets slow on larger inputs.
```
Array[          (* Create a table of *)
(+##)^2-##&,    (* (i+j)^2-i j *)
{#,#}+1,0]      (* for i,j = 0...k *)
~Position~#-1   (* and find where that expression equals k *)
```
]