Find all integer pairs that produce a given Loeschian number

Question

Inspired by and drawns from Is this number Loeschian?

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 (see illustrations for \$k = 4\$ and for \$k = 7\$). (Because of this property, these numbers find application in mobile cellular communication networks.)

• See more characterizations in the OEIS page 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...

¹while (\$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?

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 1

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!

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 2

Haskell, 46 bytes

f x|l<-[0..x]=[(i,j)|i<-l,j<-l,i*i+j*j+i*j==x]

Try it online!

Answer 3

Husk, ¹⁵ 13 bytes

fo=¹§+Πṁ□π2…0

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 4

Jelly, 11 10 9 bytes

Żp`ḋÄ$=¥Ƈ

Try it online!

Outputs [] for non-Loeschian numbers

-1 byte thanks to Sisyphus

Not particularly efficient, but that can be fixed for an additional 2 bytes.

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 5

Wolfram Language (Mathematica), 46 bytes

Solve[i^2+j^2+i*j==#&&i>=j>=0,{i,j},Integers]&

Try it online!

Answer 6

Scala, 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!

Thanks to user for -2

Answer 7

Japt, 13 bytes

ô ï f@¶Xx²+X×

Try it

ô ï 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 8

Perl 5, 81 77 74 bytes

sub{map{$i=$_;grep{$k==$i**2+$i*$_+$_**2&&($_=[$i,$_])}$i..$k}0..($k=pop)}

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 9

Scala, 53 bytes

| =>for(i<-0 to|;j<-0 to|if| ==i*i+j*j+i*j)yield(i,j)

Try it online!

Answer 10

C (gcc), 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!

Answer 11

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!

---

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!

Answer 12

JavaScript (V8), 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!

Answer 13

Charcoal, 25 bytes

ＮθＦ⊕₂θＦ⊕ι¿⁼θ⁻Ｘ⁺ικ²×ικＩ⟦ικ

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 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! Very slow, so don't try anything over about 500.

Answer 14

Wolfram Language (Mathematica), 42 bytes

Array[(+##)^2-##&,{#,#}+1,0]~Position~#-1&

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 *)