[Question]
  [
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\$); ... 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 so that corresponding cells are the same distance apart (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 challenge
Given a **positive integer**, output a truthy result **if it is a Loeschian number**, or a falsy result otherwise.
The program or function should handle (say in less than a minute) inputs up to \$1000\$, or up to data type limitations.
Code golf. Shortest wins.
## Test cases
The following numbers should output a truthy result:
```
1, 4, 7, 12, 13, 108, 109, 192, 516, 999
```
The following numbers should output a falsy result:
```
2, 5, 10, 42, 101, 102, 128, 150, 501, 1000
```
  
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~11~~ 9 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
```
ÆF‘%3,2ḄȦ
```
[Try it online!](http://jelly.tryitonline.net/#code=w4ZG4oCYJTMsMuG4hMim&input=&args=OTk5) or [verify all test cases](http://jelly.tryitonline.net/#code=w4ZG4oCYJTMsMuG4hMimCsOH4oKs4oKsRw&input=&args=WzEsIDQsIDcsIDEyLCAxMywgMTA4LCAxMDksIDE5MiwgNTE2LCA5OTldLCBbMiwgNSwgMTAsIDQyLCAxMDEsIDEwMiwgMTI4LCAxNTAsIDUwMSwgMTAwMF0).
### Background
In [Elementary results on the binary quadratic form *a² + ab + b²*](https://arxiv.org/abs/math/0408107), the author proves the following theorem about Löschian numbers.
> 
> **Theorem 16.** *The necessary and sufficient condition of any non-negative integer to be in the form* a² + ab + b² *is that, in its prime factorization, all primes other than* 3 *that are not in the form* (6k + 1) *have even exponents.*
> 
> 
> 
As noted on the relevant [OEIS page](https://oeis.org/A003136), since all integers are are congruent to **0**, **1** or **2** modulo **3**, the number **3** is the only prime that is congruent to **0**, and all numbers of the form **(6k + 1)** are congruent to **1**, the theorem can be stated alternatively as follows.
*A non-negative integer* n *is a Löschian number if and only if all prime factors of* n *that are congruent to* 2 *modulo* 3 *have even exponents.*
### How it works
```
ÆF‘%3,2ḄȦ  Main link. Argument: n (integer)
ÆF         Yield the prime factorization of n, as prime-exponent pairs.
  ‘        Increment all primes and exponents, turning primes of the form 3k - 2
           into multiples of 3 and odd exponents into multiples of 2.
   %3,2    Reduce all incremented primes/exponents modulo 3/2.
           n is Löschian if and only if this does not result in a [0, 0] pair.
           Due to Jelly's form of vectorization, this yields [3, 2] if n = 1.
       Ḅ   Unbinary; convert each pair from base 2 to integer.
           Note that [x, y] = [0, 0] if and only if 2x + y = 0.
        Ȧ  All; return 1 if the result contains no zeroes, 0 otherwise.
```
[Answer]
## [Retina](https://github.com/m-ender/retina), ~~66~~ ~~63~~ ~~45~~ ~~43~~ 36 bytes
```
^()(\1(?<1>.\1))+(\1(.(?(4).\4)))*$
```
Despite the title saying Retina, this is just a plain .NET regex which accepts [unary representations](http://meta.codegolf.stackexchange.com/questions/5343/can-numeric-input-output-be-in-unary) of Loeschian numbers.
Inputs 999 and 1000 take well under a second.
[Try it online!](http://retina.tryitonline.net/#code=JShHYAouKwokKgpeKCkoXDEoPzwxPi5cMSkpKyhcMSguKD8oNCkuXDQpKSkqJA&input=MQoyCjMKNAo1CjYKNwo4CjkKMTAKMTEKMTIKMTMKMTQKMTUKMTYKMTc) (The first line enables a linefeed-separated test suite, and the next two take care of the conversion to unary for convenience.)
### Explanation
The solution is based on the classification that the input can be written as `i*i + j*(i + j)` for positive `i` and non-negative `j` (since we don't have to handle input `0`), and that `n*n` is just the sum of the first `n` odd integers. Golfing this was an interesting exercise in forward references.
A "forward reference" is when you put a backreference inside the group it refers to. Of course that doesn't work when the group is used the first time, since there is nothing to be backreferenced yet, but if you put this in a loop, then the backreference gets the previous iteration's capture each time. This in turn, let's you build up a larger capture with each iteration. This can be used to craft very compact patterns for things like triangular numbers, squares and Fibonacci numbers.
As an example, using the fact that squares are just sums of the first `n` odd integers, we can match a square input like this:
```
(^.|..\1)+$
```
On the first iteration, `..\1` can't work, because `\1` doesn't have a value yet. So we start with `^.`, capturing a single character into group `1`. On subsequent iterations, `^.` no longer matches due to the anchor, but now `..\1` is valid. It matches two more characters than the previous iteration and updates the capture. This way we match increasing odd numbers, getting a square after each iteration.
Now unfortunately, we can't use this technique as is. After matching `i*i`, we need to get `i` as well, so that we can multiply it by `j`. A simple (but long) way to do this is to make use of the fact that matching `i*i` takes `i` iterations, so that we've captured `i` things in group `1`. We could now use [balancing groups](https://stackoverflow.com/a/17004406/1633117) to extract this `i`, but like I said that's expensive.
Instead, I figured out a different way to write this "sum of consecutive odd integers" that also yields `i` in a capturing group at the end. Of course the `i`th odd number is just `2i-1`. This gives us a way to increment the forward reference only by 1 on each iteration. That's this part:
```
^()(\1(?<1>.\1))+
```
This `()` just pushes an empty capture onto group `1` (initialising `i` to `0`). This is pretty much equivalent to the `^.|` in the simple solution above, but using `|` in this case would be a bit trickier.
Then we have the main loop `(\1(?<1>.\1))`. `\1` matches the previous `i`, `(?<1>.\1)` then updates group `1` with `i+1`. In terms of the *new* `i`, we've just matched `2i-1` characters. Exactly what we need.
When we're done, we've matched some square `i*i` and group `1` still holds `i` characters.
The second part is closer to the simple square matching I showed above. Let's ignore the backreference to `1` for now:
```
(.(?(4).\1))*
```
This is basically the same as `(^.|..\4)*`, except that we can't make use of `^` because we're not at the start of the string. Instead we make use of a conditional, to match the additional `.\1` only when we've already used group `4`. But in effect this is exactly the same. This gives us `j*j`.
The only thing that's missing is the `j*i` term. We combine this with the `j*j` by making use of the fact that the `j*j` computation still takes `j` iterations. So for each iteration we also advance the cursor by `i` with `\1`. We just need to make sure not to write that into group `4`, because that would mess with matching consecutive odd numbers. That's how we arrive at the:
```
(\1(.(?(4).\1)))*
```
[Answer]
## CJam (16 15 bytes)
```
{mF{~\3%2=&},!}
```
[Online demo](http://cjam.aditsu.net/#code=%5B1%204%207%20%2012%2013%20%20108%20109%20192%20516%20999%0A%202%205%2010%2042%20101%20102%20128%20150%20501%201000%5D%0A%0A%7BmF%7B~%5C3%252%3D%26%7D%2C!%7D%0A%0A%25%60)
This is a block (an "anonymous function") which takes input on the stack and leaves `0` or `1` on the stack. It uses the characterisation that a number is Loeschian iff it has no prime factor equal to 2 mod 3 with odd multiplicity.
Thanks to [Dennis](https://codegolf.stackexchange.com/users/12012/dennis) for a one-byte saving.
[Answer]
# [Python 2](https://docs.python.org/2/), 49 bytes
```
lambda n:0in[(n-3*i*i+0j)**.5%1for i in range(n)]
```
[Try it online!](https://tio.run/##RcmxDkAwFAXQ3Ve8RaIVUqRLE1@CoUJ5wiVNF19fTNZzrjusJ@ro2j7u9hgnSzCK0WUoGsmSc7UJKUudVu70xMQgb7HMGcQQP8JPWgmTEDty7xq6PCMQ4gM "Python 2 – Try It Online")
Uses the equivalent quadratic form given on OEIS of `n == 3*i*i+j*j`. Check whether `n-3*i*i` is a perfect square for any `i` by taking its square root and checking if it's an integer, i.e. equals 0 modulo 1. Note that Python computes square roots of perfect squares exactly, without floating point error. The `+0j` makes it a complex number to avoid an error on the square root of a negative.
[Answer]
# Python 2, 56 bytes
```
lambda n:any(n==i*i%n+i/n*(i/n+i%n)for i in range(2*n*n))
```
[Answer]
## Haskell, 42 bytes
```
f k=or[k==i*i+j*j+i*j|i<-[0..k],j<-[0..i]]
```
Usage example: `f 501` -> `False`.
Tries all combinations of `i` from `0` to `k` and `j` from `0` to `i` . `or` returns `True` if the equality `k==i*i+j*j+i*j` holds for at least one of the combinations.
@flawr found a slightly different version with the same byte count:
```
f k|v<-[0..k]=or[(i+j)^2==k+i*j|i<-v,j<-v]
```
[Answer]
## Java 8, 81 bytes
```
k->{for(int i=0,j;i<=k;i++)for(j=0;j<=k;)if(i*i+j*j+i*j++==k)return 1;return 0;};
```
simple, naïve implementation. coincidentally same code as C# but uses `->` rather than `=>`.
[Answer]
# [Jellyfish](https://github.com/iatorm/jellyfish), ~~56~~ ~~43~~ ~~41~~ ~~29~~ 28 bytes
**2 bytes thanks to Zgarb**
```
p
n    <
+`/
`1*
/
+
&*r&;>i
```
[Try it online!](http://jellyfish.tryitonline.net/#code=cApuICAgIDwKK2AvCmAxKgovCisKJipyJjs-aQ&input=MTM)
A fork of [my Jelly answer](https://codegolf.stackexchange.com/a/88733/48934).
[Answer]
# [MATL](https://github.com/lmendo/MATL), ~~14~~ 13 bytes
```
t:0hU&+HM&*+m
```
[Try it online!](http://matl.tryitonline.net/#code=dDowaFUmK0hNJiorbQ&input=OTk5) Or [verify all test cases](http://matl.tryitonline.net/#code=IkAKdDowaFUmK0hNJiorbQpddiE&input=WzEsIDQsIDcsIDEyLCAxMywgMTA4LCAxMDksIDE5MiwgNTE2LCA5OTksIDIsIDUsIDEwLCA0MiwgMTAxLCAxMDIsIDEyOCwgMTUwLCA1MDEsIDEwMDBd).
Outputs `1` or `0`.
### Explanation
```
t:    % Implicitly input number k. Duplicate. Generate vector [1 2 ...k]
0h    % Concatenate a 0. Gives [1 2 ... k 0]
U     % Square, element-wise. Gives [1 4 ... k^2 0]
&+    % Sum of all pairs from this vector. Gives a (k+1)×(k+1) matrix
HM    % Push [1 2 ... k 0] again
&*    % Product of all pairs from this vector. Gives a (k+1)×(k+1) matrix
+     % Add the two matrices
m     % True if k is a member of the resulting matrix. Implicitly display
```
[Answer]
## Python, 67 bytes
```
lambda k,r=range:any(i*i+j*j+i*j==k for i in r(k+1)for j in r(k+1))
```
<https://repl.it/Cj6x>
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~15~~ ~~14~~ ~~13~~ 12 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
**1 byte thanks to miles.**
```
²S+P
‘ṗ2’Ç€i
```
[Try it online!](http://jelly.tryitonline.net/#code=wrJTK1AK4oCY4bmXMuKAmcOH4oKsaQ&input=&args=MTI)
[Verify the smaller testcases](http://jelly.tryitonline.net/#code=wrJTK1AK4oCY4bmXMuKAmcOH4oKsaQrDh-KCrOKCrA&input=&args=W1sxLCA0LCA3LCAxMiwgMTNdLFsyLCA1LCAxMCwgNDJdXQ).
A word of advice when testing for large numbers (bigger than 50): don't.
Truthy is a positive number. Falsey is zero.
### Explanation
```
‘ṗ2’Ç€i   main chain, argument: z
‘ṗ2’      generate all pairs of numbers between 0 and z inclusive
    Ç€    apply the helper link to each pair
      i   find the index of z in the result
²S+P   helper link, argument: [x,y] (a pair of numbers)
²      compute [x*x, y*y]
 S     x*x+y*y
  +P   x*x+y*y+x*y
```
[Answer]
# [Brachylog](https://github.com/JCumin/Brachylog), 13 bytes
```
ḋḅ{h%₃<2|~j}ᵐ
```
[Try it online!](https://tio.run/##SypKTM6ozMlPN/r//@GO7oc7WqszVB81NdsY1dRl1T7cOuH/f0tLSwA "Brachylog – Try It Online")
### Explanation
Uses the prime factors approach.
```
ḋḅ{h%₃<2|~j}ᵐ
ḋ              Prime factorization: gives a list like [7, 3, 2, 2]
 ḅ             Analyze into blocks of the same value: [[7], [3], [2, 2]]
  {        }ᵐ  This predicate must succeed for each block:
   h            The head (first number in the list, i.e. the prime factor)
    %₃           mod 3
      <2         is less than 2
        |       Or
         ~j      the list can be considered the result of joining two identical sublists
                 (i.e. its length is even, i.e. that prime factor has an even exponent)
```
[Answer]
# C (gcc), 71 69 bytes
```
i,j,r;f(n){for(r=i=n+1;i--;)for(j=n;j--;)r*=n!=i*i+j*j+i*j;return!r;}
```
[Answer]
## VBA, ~~68~~ 67 bytes
```
Function L(N):For a=0To N:For b=0To a:L=L+(N=a^2+a*b+b^2):Next b,a
```
Naive search, starting to slow down slightly for n=1000. Excel recognizes zero return as falsy, all other returns as truthy.
Note that investigation of negative *i* and *j* is not needed, since given *i>j>=0* :
> 
> *(-i)2 + (-i)(-j) + (-j)2 = i2 + ij + j2*
> 
> 
> 
(the same result as for *i* and *j*)
> 
> *(-i)2 + (-i)j + j2 = i2 - ij + j2*
> 
> 
> *i2 + i(-j) + (-j)2 = i2 - ij + j2*
> 
> 
> 
(if one is negative, it doesn't matter which one), and then
> 
> *(i-j)2 + (i-j)j + j2 = (i2 - 2ij + j2) + (ij - j2) + j2 = i2 - ij + j2*
> 
> 
> 
And since both *(i-j)* and *j* are non-negative, any generation of Loeschian numbers involving a negative number can be achieved using non-negative numbers.
---
*Saved a byte, `Next:Next` -> `Next b,a` thanks to Taylor Scott.*
[Answer]
## C#, ~~84~~ ~~82~~ 81 bytes
```
k=>{for(int i=0,j;i<=k;++i)for(j=0;j<=k;)if(i*i+j*j+i*j++==k)return 1;return 0;};
```
A naïve solution. 1 = true, 0 = false
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E/wiki/Commands), 10 [bytes](https://github.com/Adriandmen/05AB1E/wiki/Codepage)
```
ÝãDnOsP+Iå
```
[Try it online](https://tio.run/##yy9OTMpM/f//8NzDi13y/IsDtD0PL/3/39LSEgA) or [verify all test cases](https://tio.run/##yy9OTMpM/V9Waa@k8KhtkoKSfeX/w3MPL3bJ8y8O0K48vPS/zv9oQx0THXMdQyMdQ2MdQwMLILbUMbQ00jE1NNOxtLTUUQAygYI6JkAVBoZADKSNgMpMDXRMwXwDg1gA).
**Explanation:**
```
Ý           # Push a list in the range [0, (implicit) input]
 ã          # Create all possible pairs by taking the cartesian product on itself
  D         # Duplicate this list of pairs
   n        # Square each inner value in the pairs: [i²,j²]
    O       # Sum each inner pair together: i²+j²
     s      # Swap so the original list of pairs is at the top of the stack again
      P     # Take the product of each inner pair: i*j
       +    # Sum the values in the list at the same positions: i²+j²+i*j
        Iå  # And check if the input-integer is in this list
            # (after which the result is output implicitly)
```
[Answer]
# [Desmos](http://desmos.com/calculator), ~~62~~ 60 bytes
```
f(k)=\prod_{n=0}^k\prod_{m=0}^k\left\{mm+nn+mn=k:0,1\right\}
```
[![](https://i.stack.imgur.com/2IjGA.png)](https://i.stack.imgur.com/2IjGA.png)
Outputs \$0\$ for truthy and \$1\$ for falsey
[Try It On Desmos!](https://www.desmos.com/calculator/sl7klzxsih)
Explanation:  
The products(\$\prod\$) iterate through \$0\$ through \$k\$ inclusive for both \$m\$ and \$n\$. It then tests if \$m^2+n^2+mn=k\$. If it does, then we know that \$k\$ is a Loeschian number, and we multiply by \$0\$. Else, we multiply by \$1\$. So if we don't encounter any pairs of \$m\$ and \$n\$ that satisfy \$m^2+n^2+mn=k\$, it will keep on multiplying \$1\$, and therefore it will stay at \$1\$ at the end. So the falsey value is \$1\$. But if a pair does satisfy \$m^2+n^2+mn=k\$, then it multiplies by \$0\$, and because \$0\$ times anything is \$0\$, and we are essentially just multiplying over and over again, it is guaranteed that once we multiply by \$0\$, it will stay at \$0\$. Hence, the truthy value is \$0\$.
[Answer]
# Javascript (using external library - Enumerable) (63 bytes)
```
k=>_.Range(0,k+1).Any(i=>_.Range(0,k+1).Any(j=>i*i+j*j+i*j==k))
```
Link to library: <https://github.com/mvegh1/Enumerable>
Code explanation: Create a range of integers from 0 to k (call this the "i" range), and test if any "i" satisfies a certain predicate. That predicate creates a range from 0 to k (call this the "j" range), and tests if any "j" satisfies a certain predicate. That predicate is the loeschian formula
[![enter image description here](https://i.stack.imgur.com/Me2Wu.png)](https://i.stack.imgur.com/Me2Wu.png)
[Answer]
# [Perl 6](http://perl6.org), ~~52 51~~ 50 bytes
```
->\k{?first ->(\i,\j){k==i*i+j*j+i*j},(0..k X 0..k)}
->\k{?grep ->(\i,\j){k==i*i+j*j+i*j},(0..k X 0..k)}
{?grep ->(\i,\j){$_==i*i+j*j+i*j},(0..$_ X 0..$_)}
```
### Explanation:
```
{
  # Turn the following into a Bool
  # ( Technically not necessary as a list of 1 or more values is truthy )
  ?
  # find all where the code block returns a truthy value
  grep
  # pointy block that takes one value (list of 2 values)
  # and gives each of the values in it a name
  ->
    $ ( \i, \j )
  {
    # return true if the definition matches
    $_ == i*i + j*j + i*j
  },
  # a list of 2 element lists (possible i and j values)
  ( 0..$_ X 0..$_ )
}
```
### Test:
```
use v6.c;
use Test;
my @true = 0, 1, 4, 7, 12, 13, 108, 109, 192, 516, 999;
my @false = 2, 5, 10, 42, 101, 102, 128, 150, 501, 1000;
plan (@true + @false) * 2;
my &is-loeschian = {?grep ->(\i,\j){$_==i*i+j*j+i*j},(0..$_ X 0..$_)}
for |(@true X True), |(@false X False) -> ( $input, $expected ) {
  my ($result,$seconds) = $input.&time-it;
  is $result, $expected, ~$input;
  cmp-ok $seconds, &[<], 60, "in $seconds seconds"
}
sub time-it ( $input ) {
  my $start = now;
  my $result = $input.&is-loeschian;
  my $finish = now;
  return ( $result, $finish - $start )
}
```
```
1..42
ok 1 - 0
ok 2 - in 0.00111763 seconds
ok 3 - 1
ok 4 - in 0.00076766 seconds
...
ok 19 - 516
ok 20 - in 0.19629727 seconds
ok 21 - 999
ok 22 - in 0.1126715 seconds
ok 23 - 2
ok 24 - in 0.0013301 seconds
ok 25 - 5
ok 26 - in 0.00186610 seconds
...
ok 37 - 150
ok 38 - in 0.83877554 seconds
ok 39 - 501
ok 40 - in 9.2968558 seconds
ok 41 - 1000
ok 42 - in 37.31434146 seconds
```
[Answer]
## PowerShell v2+, ~~63~~ ~~56~~ 55 bytes
```
param($k)(0..$k|%{0..($i=$_)|%{$i*($i+$_)+$_*$_}})-eq$k
```
Takes input `$k`, loops upwards twice (outer loop `$i = 0 to $k`, inner loop `$j = 0 to $i`), each iteration generates the result of `i*i + j*j + i*j` (shortened to `i*(i+j) + j*j`). Those results are encapsulated in parens, and passed as an array to `-eq$k`. This acts as a filter to select only elements that equal the input. Outputs a nonzero (the number back) for truthy, or nothing (empty) for falsey. Processes `1000` in about 15 seconds on my machine.
### Test Cases
```
PS C:\Tools\Scripts\golfing> (1,4,7,12,13,108,109,192,516,999|%{.\loeschian-numbers.ps1 $_})-join','
1,4,7,12,13,108,109,192,516,999
PS C:\Tools\Scripts\golfing> (2,5,10,42,101,102,128,150,501,1000|%{.\loeschian-numbers.ps1 $_})-join','
PS C:\Tools\Scripts\golfing>
```
[Answer]
# Perl, 54 + 1 (`-n` flag) = 55 bytes
```
for$i(0..$_){for$j(0..$_){$i*$i+$j*$j+$i*$j-$_?1:say}}
```
Needs `-n` and `-M5.010` flags to run :
```
perl -nE 'for$i(0..$_){for$j(0..$_){$i*$i+$j*$j+$i*$j-$_?1:say}}'
```
Outputs some stuffs if the number is a Loeschian number, and nothing otherwise.
This implementation is quite boring, so here is another one, for 87 bytes, regex-based, just for the eyes :
```
perl -pE '$_=(1 x$_)=~/^(.*)(??{$1x(-1+length$1)})(.*)(??{$2x(-1+length$2)})(??{$1x length$2})$/'
```
Carefull with this one, as the backtracking will use a lot of memory, so don't try to test numbers too big! (especially numbers that aren't Loeschians)
[Answer]
# [Dyalog APL](http://goo.gl/9KrKoM), 19 [bytes](https://codegolf.meta.stackexchange.com/a/9429/43319)
```
⊢∊(∘.(×-⍨2*⍨+)⍨0,⍳)
```
Checks if *k* ∊ (*i* + *j*)² – *ij*, for any 0 ≤ *i*, *j* ≤ *k*.
    `⊢` is *k*  
`∊` a member of  
    `∘.` all combinations of  
        `×` *i* times *j*  
        `-⍨` subtracted from  
        `2*⍨` the square of  
        `+` *i* plus *j*  
    `⍨` for all *i* **and** *j* in  
    `0,` zero prepended to  
    `⍳` the integers 1 through *k*
1000 takes 3.3 seconds on my M540 and even less on [TryAPL](http://tryapl.org/?a=f%u2190%u22A2%u220A%28%u2218.%28%D7-%u23682*%u2368+%29%u23680%2C%u2373%29%20%u22C4%20%u2395TS%20%u22C4%20f%20999%20%u22C4%20%u2395TS%20%u22C4%20f%201000%20%u22C4%20%u2395TS&run).
[Answer]
## Matlab, ~~53~~ 52 bytes
```
n=input('');[a b]=ndgrid(0:n);find((a+b).^2-a.*b==n)
```
Simple search over all possibilities.  
Outputs empty array as falsy and a non-empty vector as truthy value.
Considering all-zeros matrix as falsy and not-all-zeros matrix as truthy we can get rid of the `find` function resulting in **~~47~~ 46 bytes** solution:
```
n=input('');[a b]=ndgrid(0:n);(a+b).^2-a.*b==n
```
*One byte saved thanks to [@flawr](https://codegolf.stackexchange.com/users/24877/flawr)*
[Answer]
# C, 66 bytes
Call `f()` with the number to test. The function returns the number of solutions it found.
```
q,r;f(n){for(r=q=0;q++<n*n;r+=n==q%n*(q%n+q/n)+q/n*q/n);return r;}
```
[Try it on ideone](http://ideone.com/fDyGD8).
[Answer]
# Mathematica, 44 bytes
```
MemberQ[(+##)^2-##&@@@0~Range~#~Tuples~2,#]&
```
Unnamed function taking an integer as input and returning `True` or `False`. The command `0~Range~#~Tuples~2` creates all ordered pairs of integers both between `0` and the input `#`. The function `(+##)^2-##&` computes the square of the sum of its arguments minus the product of its arguments; when called on two arguments `i` and `j`, this is exactly `i^2+j^2+ij` as desired. So that function is called on all the tuples, and then `MemberQ[...,#]` checks whether the input is one of the resulting values.
[Answer]
# ASP, 39 + 4 = 43 bytes
```
o:-k=I*I+J*J+I*J;I=1..k;J=1..k.:-not o.
```
Output: the problem is satisfiable iff k is Loeschian.
[Answer Set Programming](https://en.wikipedia.org/wiki/Answer_set_programming) is a logical language, similar to prolog.
I use here the [Potassco implementation](https://potassco.org), *clingo*.
Input is taken from parameters (`-ck=` is 4 bytes long). Call example:
```
clingo -ck=999
```
Output sample:
```
SATISFIABLE
```
Tried with 1000:
```
clingo -ck=1000
```
Output sample:
```
UNSATISFIABLE
```
You can try it in your [browser](https://potassco.org/clingo/run/) ; unfortunately, this method doesn't handle call flags, so you need to add the line `#const k=999` in order to make it work.
---
Ungolfed & explained code:
```
v(1..k).  % predicate v(X) holds for any X in [1..k]
o:- k=I*I+J*J+I*J ; v(I) ; v(J).  % o holds if k is Loeschian.
:- not o.  % discard models where o doesn't holds (make problem unsatisfiable)
```
[Answer]
# [Add++](https://github.com/cairdcoinheringaahing/AddPlusPlus), ~~28~~ 26 bytes
```
D,g,@@,*aaz€b*Fs
L,0rd‽gAe
```
[Try it online!](https://tio.run/##S0xJKSj4/99FJ13HwUFHKzGx6lHTmiQtt2IuHx2DopRHDXvTHVP///9v/i@/oCQzP6/4v27mf13d4oLU5My0StucxNyklEQFQwA "Add++ – Try It Online") or [verify the test cases](https://tio.run/##S0xJKSj4/99FJ13HwUFHKzGx6lHTmiQtt2IuHx2DopRHDXvTHVP/l1hFGyqYKJgrGBopGBorGBpYALGlgqGlUSxXmlW0kYIpkK9gApQ0MARiIG0EVGFqEMvF5Vqiw8WZkgkilHTtlIC0Sk5iblJKooKhHUjUn4vLX0kJqC6NgLr/AA)
Takes a long time for large inputs. Times out on TIO for the `501` and larger test cases.
## How it works
```
D,g,@@,		; Define a helper function that takes 2 arguments, i and j
		; Example: i = 0, j = 2			STACK = [0 2]
	*	; Multiply				STACK = [0]
	aaz	; Push [i i] and [j j]			STACK = [0 [0 0] [2 2]]
	€b*	; Take their products			STACK = [0 [0 4]]
	Fs	; Flatten and return the sum		STACK = [0 0 4 4]
		; Returns i*i + i*j + j*j		Returns 4
L,		; Main function, takes 1 argument, k
		; Example: k = 4			STACK = [4]
	0r	; Range from 0				STACK = [[0 1 2 3 4]]
	d	; Duplicate				STACK = [[0 1 2 3 4] [0 1 2 3 4]]
	‽	; Across all pairs...			STACK = [[[0 0] [0 1] [0 2] ... [4 2] [4 3] [4 4]]]
	    g	;   run the helper function		STACK = [[0 1 4 ... 28 37 48]]
	Ae	; Is k in the array?			Returns 1
```
[Answer]
# PHP, 70 bytes
```
for(;$i++<$k=$argv[1];)for($j=$i+1;$j--;)$i*$i+$j*$j+$i*$j-$k?:die(1);
```
takes input from command line argument; exits with `1` for Loeschian number, with `0` else.  
Run with `-nr`.
**breakdown**
```
for(;$i++<$k=$argv[1];)     # loop $i from 1 to $k
    for($j=$i+1;$j--;)      # loop $j from $i to 0
        $i*$i+$j*$j+$i*$j-$k?   # if $i,$j,$k do not satisfy the equation, do nothing
        :die(1);                # else exit with return code 1
                            # implicit: exit with code 0
```
---
# revisited, 67 bytes
by merging the loops to one:
```
for($b=1+$a=$n=$argv[1];$b--||$b=--$a;)$a*$a+$b*$b+$a*$b-$n||die(1);
```
not tested, though
pushed it a little further ...
## 57 bytes:
```
while($a*$a+$b*$b+$a*$b-$argn)$b++<$a||$b=1/($a++<$argn);
```
* takes input from pipe; run with -nR
* no output if input is Loeschian, division by zero error if not
* not tested either
Haven´t been here for a while. I thought I remembered that deliberate errors were forbidden; but I couldn´t find it. What I found instead was a deleted post in the [loopholes thread](https://codegolf.meta.stackexchange.com/questions/1061) that requested to forbid output to STDERR. Could someone refresh my memory?
**breakdown**
```
# $a and $b are implicitly initialized with NULL, which increments to 1
while($a*$a+$b*$b+$a*$b-$argn)      # while solution not found
    $b++<$a                         # inner loop: $b from 1 to $a
    ||$b=1/($a++<$argn);            # outer loop: $a from 0 to $argn(=input)
    # $a<$argn: 1/true -> 1
    # else: 1/false -> division by zero
```
[Answer]
# GAP 4.7.9, 75 bytes
```
k->Filtered(Collected(FactorsInt(k)),d->RemInt(d[2],2)+RemInt(d[1],3)=3)=[]
```
Every positive integer \$k\$ has a unique representation as the product of powers of distinct primes, \$p\_1^{e\_1}\dots p\_s^{e\_s}\$. \$k\$ is Loeschian iff no prime \$p\_i\$ occurs to an odd power \$e\_i\$ if \$p\_i=2\$ modulo 3.
The above GAP code is for a lambda-expression which returns `true` if `k` is Loeschian and `false` if not.
[Try it online!](https://tio.run/##S08s@B9QlJlXouGWmVOSWpSaohFtpKdnYRir8z9b1w4u6Jyfk5OaXAJkuSUml@QXFXsCtWRrauqk6NoFpeaCeCnRRrE6RpracC7QDGNNWyCKjv2vqWn9HwA "GAP – Try It Online")
Its algorithm is nice and simple but the long names for some of those GAP functions make for poor golf, unfortunately.
[Answer]
# [Japt](https://github.com/ETHproductions/japt), 13 [bytes](https://en.wikipedia.org/wiki/ISO/IEC_8859-1)
```
ô ï d@¶Xx²+X×
```
[Try it](https://petershaggynoble.github.io/Japt-Interpreter/?v=1.4.6&code=9CDvIGRAtlh4sitY1w&input=MTIKLVE)
]