SAKO, 93 91 bytes

CZYTAJ:N
*)GDYMOD(N,M*2)=0:1,INACZEJ2
2)POWTORZ:M=N(-1)1
1)DRUKUJ(9,0):M,N/M*2
STOP1
KONIEC

Perl 5 -n, 39 bytes

($c=sqrt$_/++$b)=~/\./&&redo;say"$c $b"

Try it online!

Brachylog, 7 bytes

~×Ċ√ᵗℕᵐ

Outputs [b,a]. Try it online!

Explanation

~×Ċ√ᵗℕᵐ
~×       Un-multiply the input number into a list of multiplicands
  Ċ      The list must have exactly 2 elements
   √ᵗ    Take the square root of the second element
     ℕᵐ  Both elements must be nonnegative integers

The order in which ~× tries options guarantees that the first element will be as small as possible.

APL(NARS), 39 chars

{⍵≤1:1 1⋄(√,⍵×÷)×/∊{0=2∣≢⍵:⍵⋄1↓⍵}¨⊂⍨π⍵}

factorize the input number ⍵, make partition of its factors, return the list of list of factors in the way each list has even lenght (exaple 8 has (2 2 2) get the list (2 2) because the lenght of (2 2 2) is odd) multiply them all in N, return the sqrt(N), ⍵/N.

Test:

  f←{⍵≤1:1 1⋄(√,⍵×÷)×/∊{0=2∣≢⍵:⍵⋄1↓⍵}¨⊂⍨π⍵}
  f¨(⍳8),36,4×9×25×29
┌10─────────────────────────────────────────────────────────────────────┐
│┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2───┐ ┌2─────┐│
││ 1 1│ │ 1 2│ │ 1 3│ │ 2 1│ │ 1 5│ │ 1 6│ │ 1 7│ │ 2 2│ │ 6 1│ │ 30 29││
│└~───┘ └~───┘ └~───┘ └~───┘ └~───┘ └~───┘ └~───┘ └~───┘ └~───┘ └~─────┘2
└∊──────────────────────────────────────────────────────────────────────┘
  f 8×9×7
┌2────┐
│ 6 14│
└~────┘
    

sed -E, 173 bytes

Takes unary input of xs. Outputs a and b, separated by a space.

s/^(x+)\1+?$/\1 &/
s/x+/&A&B/
:Z
s/^(x+)(Ax*)xB(x*) \1\3/\1\2B\3\1 \1\3/
tZ
s/(x+)AB(x+) (\2+$)/\1 C\2 \3/
tO
s/^(x+)x\w+( \1+$)/\1A\1B\2/
tZ
:O
s/C(x+) \1/xC\1 /
tO
s/C.+//

Try it online!

Ungolfed, 459 bytes

# find largest a such that N is a times m
s/^(x+)\1+?$/\1 &/

# calculate a times a
s/x+/&A&B/
## a A iterator B sum
## tZero is unnecessary
:Zero
## add a to sum and decrement iterator
## obtw stop iteration when a+sum > N
s/^(x+)(Ax*)xB(x*) \1\3/\1\2B\3\1 \1\3/
tZero
## go to next step if iterator is zero
s/(x+)AB(x+) (\2+$)/\1 C\2 \3/
tOne
## Reduce a to try again
s/^(x+)x\w+( \1+$)/\1A\1B\2/
tZero

# find b that N equals a times a times b
:One
s/C(x+) \1/xC\1 /
tOne

# clean up
s/C.+//

Try it online!

Jelly, 6 bytes

²×¥þœi

Try It Online!

TFW you outgolf Dennis by 33% :P

Outputs [b, a].

This does not use any prime factorization built-ins. Instead, it uses an array-oriented built-in œi (multi-dimensional index of), which is interestingly not found in any of APL/J/K/BQN.

How it works

²×¥þœi    Monadic link. Input: the positive integer n.
²×¥þ      Compute the table of (x^2 * y) for both x and y in 1..n
          For n = 3, this returns:
          × | 1  4  9
          --+--------
          1 | 1  4  9
          2 | 2  8 18
          3 | 3 12 27
    œi    Find the multi-dimensional index of n in the matrix above
          [3, 1]

Since it finds the first occurrence in the matrix, b (first index) is guaranteed to be minimized, and n always appears at least once in the matrix at the coordinates [n, 1].

Vyxal, 9 bytes

K'∆²;t₍/√

Try it Online!

Very cool and good. Outputs as [constant, surd]

Explained

K'∆²;t₍√/
K        # from the factors of the input, 
 '∆²;    # keep only square numbers
     t   # and keep the last.
      ₍√/ # Push [sqrt(^), input / ^]

J, 19 bytes

(0|:0 2#:])&.(_&q:)

Try it online!

Same as the Jelly solution.

Python, 74 Bytes

def e(k):a=filter(lambda x:k/x**2*x*x==k,range(k,0,-1))[0];return a,k/a**2

Straightforward enough.

JavaScript (ECMAScript 6), 35 bytes

f=(n,k=n)=>n/k%k?f(n,--k):[k,n/k/k]

JavaScript 1+, 37 B

for(k=n=prompt();n/k%k;--k);[k,n/k/k]

APL, 25 chars

 {(⊢,⍵÷×⍨)1+⍵-0⍳⍨⌽⍵|⍨×⍨⍳⍵}

In English:

• 0⍳⍨⌽⍵|⍨×⍨⍳⍵: index of the largest of the squares up to n that divides completely n;
• 1+⍵-: the index is in the reversed array, so adjust the index
• (⊢,⍵÷×⍨): produce the result: the index itself (a) and the quotient b (that is, n÷a*a)

Test:

     ↑{(⊢,⍵÷×⍨)⊃z/⍨0=⍵|⍨×⍨z←⌽⍳⍵}¨⍳36
1  1
1  2
1  3
2  1
1  5
1  6
1  7
2  2
3  1
1 10
1 11
2  3
1 13
1 14
1 15
4  1
1 17
3  2
1 19
2  5
1 21
1 22
1 23
2  6
5  1
1 26
3  3
2  7
1 29
1 30
1 31
4  2
1 33
1 34
1 35
6  1

05AB1E, 14 bytes

Lv¹ynÖi¹yn/y‚ï

Explained

Lv              # for each x in range(1,N) inclusive
  ¹ynÖi         # if N % x^2 == 0
       ¹yn/y‚ï  # create [N/x^2,x] pairs, N=12 -> [12,1] [3,2]
                # implicitly output last found pair

Try it online

Python 2.7 (ungolfed) - 181 Bytes

def e(n):   
 for x in range(1,n+1):
  r=(1,x)
  for i in range(1,x+1):
   l=i**.5
   for j in range(1,x+1): 
    if i*j==x and l%1==0 and j<r[1]:r=(int(l),j)                
  print x,r

Run as: e(number) eg. e(24)

Sample output:

>> e(24)
1 (1, 1)
2 (1, 2)
3 (1, 3)
4 (2, 1)
5 (1, 5)
6 (1, 6)
7 (1, 7)
8 (2, 2)
9 (3, 1)
10 (1, 10)
11 (1, 11)
12 (2, 3)
13 (1, 13)
14 (1, 14)
15 (1, 15)
16 (4, 1)
17 (1, 17)
18 (3, 2)
19 (1, 19)
20 (2, 5)
21 (1, 21)
22 (1, 22)
23 (1, 23)
24 (2, 6)

Haskell, 43> 42 bytes

Brute force solution.

Saved 1 byte thanks to Xnor

f n=[(x,y)|y<-[1..],x<-[1..n],x*x*y==n]!!0

JavaScript (ECMAScript 2016), 40 bytes

n=>{for(k=n;n%k**2;k--);return[k,n/k/k]}

Basically a JavaScript port of Dennis's Python 2 answer.

Try it on JSBin.

Note: it doesn't work in strict mode, because k is not initialized anywhere. To make it work in strict mode, k=n in the loop should be changed to let k=n.

Mathematica 34 bytes

#/.{a_ b_^_:>{a, b},_[b_,_]:>{1,b}}&

This says to replace all the input (#) according to the following rules: (1) a number, a, times the square root of b should be replaced by {a, b} and a function b to the power of whatever should be replaced by {1,b} . Note that the function assumes that the input will be of the form, Sqrt[n]. It will not work with other sorts of input.

This unnamed function is unusually cryptic for Mathematica. It can be clarified somewhat by showing its full form, followed by replacements of the original shorter forms.

Function[
   ReplaceAll[
      Slot[1],
      List[
         RuleDelayed[Times[Pattern[a,Blank[]],Power[Pattern[b,Blank[]],Blank[]]],List[a,b]],
         RuleDelayed[Blank[][Pattern[b,Blank[]],Blank[]],List[1,b]]]]]

which is the same as

   ReplaceAll[
      #,
      List[
         RuleDelayed[Times[Pattern[a,Blank[]],Power[Pattern[b,Blank[]],Blank[]]],List[a,b]],
         RuleDelayed[Blank[][Pattern[b,Blank[]],Blank[]],List[1,b]]]]&

and

ReplaceAll[#, 
  List[RuleDelayed[
    Times[Pattern[a, Blank[]], 
     Power[Pattern[b, Blank[]], Blank[]]], {a, b}], 
   RuleDelayed[Blank[][Pattern[b, Blank[]], Blank[]], {1, b}]]] &

and

ReplaceAll[#, 
  List[RuleDelayed[Times[a_, Power[b_, _]], {a, b}], 
   RuleDelayed[Blank[][b_, _], {1, b}]]] &

and

ReplaceAll[#, {RuleDelayed[a_*b^_, {a, b}], RuleDelayed[_[b_, _], {1, b}]}]&

and

ReplaceAll[#, {a_*b^_ :> {a, b}, _[b_, _] :> {1, b}}] &

Matlab, 51 bytes

x=input('');y=1:x;z=y(~rem(x,y.^2));a=z(end)
x/a^2

Explanation

x=input('')       -- takes input
y=1:x             -- numbers from 1 to x
z=y(~rem(x,y.^2)) -- numbers such that their squares divide x
a=z(end)          -- biggest such number (first part of output)
x/a^2             -- remaining part

PARI/GP, 12 bytes

n->core(n,1)

core returns the squarefree part of n by default, but setting the second argument flag to 1 makes it return both parts. Output order is (b, a), e.g. (n->core(n,1))(12) -> [3, 2].

Julia, 32 bytes

\(n,k=n)=n%k^2>0?n\~-k:[k,n/k^2]

Try it online!

MATL, 12 bytes

t:U\~f0)GyU/

Try it online!

Explanation

t     % Take input n implicitly. Duplicate
:U    % Push [1 4 9 ... n^2]
\~    % True for entries that divide the input
f0)   % Get (1-based) index of the last (i.e. largest) dividing number
G     % Push input again
y     % Duplicate index of largest dividing number
U     % Square to recover largest dividing number
/     % Divide input by that. Implicitly display stack

Python 2, 43 bytes

k=n=input()
while n%k**2:k-=1
print k,n/k/k

Test it on Ideone.

Pyth, 15 bytes

/Q^
ef!%Q^T2SQ2

Test suite.

Jelly, 9 bytes

ÆE;0d2ZÆẸ

Try it online! or verify all test cases.

How it works

ÆE;0d2ZÆẸ  Main link. Argument: n

ÆE         Exponents; generate the exponents of n's prime factorization.
  ;0       Append 0 since 1ÆE returns [].
    d2     Divmod by 2.
      Z    Zip/transpose to group quotients and remainders.
       ÆẸ  Unexponent; turn the exponents of prime factorizations into integers.