TI-Basic, 16 bytes

sum(not(fPart(seq(I²/Ans,I,1,Ans-1

Port of this answer. Takes input in Ans. Output is stored in Ans and is displayed.

Vyxal rs, 3 bytes

ʁ²Ḋ

Try it Online!

-4 thanks to Aaron Miller

Arn, 12 10 bytes

·£æ9Š3nòy├

Try it!

Explained

Unpacked: +v{!(v^2%}~

Uses the formula from the OEIS page: the number of solutions to $x^2≡0 (mod n)$

              ~     1-range (inclusive) to
                _   variable initialized to STDIN; implied
            +      folded with addition after
v{                  mapping with block (key of v)
  !                 Boolean NOT
    (               Begin expression
          v         
        ^           exponentiated by
          2         two 
      %             mod
        _           implied
    )               End expression; implied
}                   End block

J, 12 bytes

1#.0=[|2^~i.

Try it online!

Ruby, 29 bytes

->n,x=n{x-=1while n%x**2>0;x}

Try it online!

APL(NARS), 23 chars, 46 bytes

{√⍵÷×/(∪a)/⍨2∣≢¨⊂⍨a←π⍵}

test:

  f←{√⍵÷×/(∪a)/⍨2∣≢¨⊂⍨a←π⍵}
  f 4
2
  f 9
3
  f 12
2
  f 13
1
  f 108
6
  f 2×2×2×2×2×3×3
12

comment:

{√⍵÷×/(∪a)/⍨2∣≢¨⊂⍨a←π⍵}
                    π⍵  factor argument
                  a←    save that in a list "a" of prime factors
                ⊂⍨      partition "a" in a list of list each element is ugual factors found
            2∣≢¨        to each element of list of list find if number of elements is odd
    ×/(∪a)/⍨           so choice in ∪a the elements appear in list of list as odd and multiple them
  ⍵÷                   divide the argument for the number of factor contained odd times
 √                     make sqrt of that
  

C (gcc), ^{44 33} 32 bytes

i;f(n){for(i=n;n%(--i*i););n=i;}

Try it online!

Wolfram Language (Mathematica), 13 bytes

√#/._^_:>1&

Try it online!

For an integer argument, √ (Sqrt) returns in the desired a√b form (unless the argument was a perfect square).

Then, /._^_:>1 matches Power expressions and replaces them with 1. As a√b expands to Times[a,Power[b,1/2]], it becomes Times[a,1]=a.

Japt -mx,  4  3 bytes

^{Crossed out &nbsp4;  is no longer 4 :)}

²vN

My first Japt answer. :)
Port of my first 5-byter 05AB1E answer, but with smart use of Japt's flags for the range and sum.

-1 byte thanks to @Shaggy thanks to the shortcuts list: p)/p␠ to ²

Try it online.

Explanation:

-m       # Convert the (implicit) input-integer to a ranged list [0, input)
  ²      # Square each value in the list, and implicitly close the function
   vN    # Check which values are divisible by the input (1 if truthy; 0 if falsey)
     -x  # After which the sum is calculated of the resulting list
         # (before the result is output implicitly)

MathGolf, 5 bytes

╒²k÷Σ

Port of my 5-byter 05AB1E answer.

Try it online.

Explanation:

╒      # Push a list in the range [1, (implicit) input]
       # (could alternatively be `r` for a range [0, input) )
 ²     # Square each value in this list
  k÷   # Check which values are evenly divisible by the input (1 if truthy; 0 if falsey)
    Σ  # And sum those checks
       # (after which the entire stack joined together is output implicitly as result)

Jelly, 5 bytes

Ḷ²%ċ0

Uses the formula from the OEIS: the number of solutions to $$x^2 equiv 0 (mathrm{mod} n)$$ Explanation:

• implicit input
• range 0..n-1,
• square each
• modulo input (I got this part to work via trial and error)
• count zeroes
• implicit print

Try it online!

Retina, 123 72 bytes

.+                        We convert the input into unary
$&*_ $&*_                 and create a copy for factor checking
{` (_+)                   START LOOP: We square the input by multiplying
$& $.1*$1                 its string representation by its length
(?=^.* (_+) (_+))2+ .+   We check if the square is a factor of the input
$.1                       if so we replace the whole text with the current counter
 (_*)_.*                  Otherwise we decrement the counter by one
 $1                       ---
-- IMPLICIT LOOP END --
-- IMPLICIT OUTPUT --

Try it online!

---

This approach is essentially a port of Kevin Cruijssen's 05AB1E answer.
It checks all the numbers from the input downwards until it finds a number whose square divides the original.

---

-51 bytes: Many thanks to Neil, whose suggestion helped me to cut off a ridiculous amount of code.

I've also switched from separating with newlines to separating with spaces because . is anti-newline.

Charcoal, 15 12 bytes

Now based on @someone's formula.

ＮθＩＬΦθ¬﹪×ιιθ

Try it online! Link is to verbose version of code. For each number from 0 to the input, calculates whether its square is divisible by the input, and takes the number of matches.

Alternative version, also 12 bytes:

ＮθＩΣＥθ¬﹪×ιιθ

Try it online! Link is to verbose version of code. For each number from 0 to the input, calculates whether its square is divisible by the input, and takes the sum of the results.

Alternative version, also 12 bytes:

ＮθＩ№Ｅθ﹪×ιιθ⁰

Try it online! Link is to verbose version of code. For each number from 0 to the input, calculates the remainder when its square is divisible by the input, and counts the number of zeros.

Retina 0.8.2, 26 bytes

.+
$*
((^1|112)+)1*$
$#2

Try it online! Link includes test cases. Explanation:

.+
$*

Convert to unary.

((^1|112)+)

Find the largest square number...

1*$

... that divides the input...

$#2

... and output its root.

Bonus 63-byte version that for an input of √1, √2, √3, √4, √5, √6, √7, √8, √9... outputs 1, √2, √3, 2, √5, √6, √7, 2√2, 3... etc. (Previous bonus version did not handle √1 correctly.)

d+
$*
r`(?=^.(3)+)(.)3*((1$|114)+)
$#4$2$#1
D1$

^1(D)
$1

Try it online!

Gaia, 5 bytes

^{(this was produced by trying a bunch of languages from https://github.com/ETHproductions/golfing-langs until I found one that had the most useful built-ins for this problem)}

dụ⁇)u

Explanation:

d     divisors
 ụ⁇   keep only squares
   )  take last
    u square root

Try it online! (example stolen from the Jelly answer)

bc, 52 bytes

define f(n){for(i=n;i--;){if(!(n%(i*i))){return i}}}

Try it online!

Jelly, 6 bytes

ÆE:2ÆẸ

A monadic Link accepting a positive integer which yields a positive integer.

Try it online! Or see the first 100.

How?

ÆE:2ÆẸ - Link: integer, X          e.g. 9587193
ÆE     - factorisation vector (X)       [0,1,0,4,3] (since 2°×3¹×5°×7⁴×11³=9587193)
  :2   - integer divide by two          [0,0,0,2,1]
    ÆẸ - evaluate factorisation vector  539         (since 2°×3°×5°×7²×11¹=539)

Haskell, 30 bytes

f n=sum[0^mod(x^2)n|x<-[1..n]]

Try it online!

Based on the solution of my pronoun is monicareinstate, counting the number of solutions to $x^2 equiv 0 (mathbb{mod} n)$ using the range from 1 to n.

Haskell, 32 bytes

f n=until((<1).mod n.(^2))pred n

Try it online!

Start with n and repeatedly take the predecessor, until it satiafies this condition: when we square it and take the original n modulo it, the result is less than 1, that is equal to 0.

CJam, 15 bytes

q~_{_*1$%!},,;

Try it online!

Uses the new method in Kevin Cruijssen's 05AB1E answer.

---

CJam, 21 bytes
This is my old approach to the problem, pre-Kevin Cruijssen's new method which I don't understand

q~mF{[~2/]}%{~#}%{*}*

Try it online!

---

Explanation

q~                     Translate input into a CJam object (allows for easier testing)
  mF                   Factorise with exponents
    {     }%           For each factor
      ~2/              Halve the exponent [and round down]
     [   ]             Capture the base & exponent in an array
            {  }%      For each transformed factor
             ~#        Expand the base and exponent into b^e
                 {*}*  Multiply all the transformed factors together

This approach removes all single factors (those that would make up the radical part) whilst halving the paired factors (equivalent to square rooting the integer part).

---

CJam, 23 bytes
A CJam port of Kevin Cruijssen's 05AB1E answer

q~_,(;{_*1$%0>!},;);

Try it online!

CP-1610 machine code,  20  17 DECLEs¹ ≈ 22 bytes

As per the exception described in this meta answer, the exact score is 21.25 bytes (170 bits)

A routine expecting the input number in R0 and returning the result in R3.

1D2     |         CLRR    R2
1C9     |         CLRR    R1
0D1     | @@loop  ADDR    R2,     R1
00A     |         INCR    R2
084     |         MOVR    R0,     R4
10C     | @@sub   SUBR    R1,     R4
10C     |         SUBR    R1,     R4
114     |         SUBR    R2,     R4
22E 004 |         BGT     @@sub
20C 001 |         BNEQ    @@next
093     |         MOVR    R2,     R3
141     | @@next  CMPR    R0,     R1
226 00D |         BLE     @@loop
0AF     |         JR      R5

How?

The CP-1610 has no multiplication, no division, no modulo. We want to implement an algorithm that relies on additions and subtractions exclusively.

We start with $k=0$. At each iteration, we update $j$ in such a way that:

$$j = frac{k(k-1)}{2}$$

The good thing about this formula is that it's very easy to compute iteratively: we just need to add $k$ to $j$ and increment $k$ afterwards.

In order to test whether $n$ is divisible by $k^2$, we initialize a variable $x$ to $n$ and subtract $k^2$ until $xle 0$.

We do not explicitly store $k^2$, but it can be easily obtained with:

$$2j+k=k(k-1)+k=k^2$$

Each time we end up with $x=0$, we update the final answer to $k$.

We stop when $j$ is greater than $n$.

Here is a link to an implementation of the algorithm in low-level JS.

Full commented test code

        ROMW    10                ; use 10-bit ROM width
        ORG     $4800             ; map this program at $4800

PNUM    QEQU    $18C5             ; EXEC routine: print a number
MULT    QEQU    $1DDC             ; EXEC routine: signed multiplication

        ;; ------------------------------------------------------------- ;;
        ;;  main code                                                    ;;
        ;; ------------------------------------------------------------- ;;
main    PROC

        SDBD                      ; set up an interrupt service routine
        MVII    #isr,   R0        ; to do some minimal STIC initialization
        MVO     R0,     $100
        SWAP    R0
        MVO     R0,     $101

        EIS                       ; enable interrupts

        MVII    #$200,  R3        ; R3 = backtab pointer
        SDBD                      ; R4 = pointer to test cases
        MVII    #@@tc,  R4

@@loop  MVI@    R4,     R0        ; R0 = next test case
        TSTR    R0                ; stop if it's 0
        BEQ     @@done

        PSHR    R4                ; save R4
        PSHR    R3                ; save R3
        CALL    pSquare           ; invoke our routine
        MOVR    R3,     R0        ; copy the result into R0
        PULR    R3                ; restore R3
        CALL    print             ; print the result
        PULR    R4                ; restore R4
        B       @@loop            ; go on with the next test case

@@done  DECR    R7                ; done: loop forever

        ;; test cases
@@tc    DECLE   4, 9, 12, 13, 108, 300, 800, 900
        DECLE   0

        ENDP

        ;; ------------------------------------------------------------- ;;
        ;;  prints the result of a test case                             ;;
        ;; ------------------------------------------------------------- ;;
print   PROC

        PSHR    R5                ; save the return address on the stack

        MVII    #4,     R1        ; R1 = number of digits
        MOVR    R3,     R4        ; R4 = backtab pointer
        ADDI    #5,     R3        ; advance by 5 characters for the next one
        PSHR    R3                ; save R3
        CLRR    R3                ; R3 = attributes (black)
        CALL    PNUM              ; invoke the EXEC routine
        PULR    R3                ; restore R3

        PULR    R7                ; return

        ENDP

        ;; ------------------------------------------------------------- ;;
        ;;  ISR                                                          ;;
        ;; ------------------------------------------------------------- ;;
isr     PROC

        MVO     R0,     $0020     ; enable display
        MVI     $0021,  R0        ; color-stack mode

        CLRR    R0
        MVO     R0,     $0030     ; no horizontal delay
        MVO     R0,     $0031     ; no vertical delay
        MVO     R0,     $0032     ; no border extension
        MVII    #$D,    R0
        MVO     R0,     $0028     ; light-blue background
        MVO     R0,     $002C     ; light-blue border
        MVO     R0,     $002C     ; light-blue border

        JR      R5                ; return from ISR

        ENDP

        ;; ------------------------------------------------------------- ;;
        ;;  our routine                                                  ;;
        ;; ------------------------------------------------------------- ;;
pSquare PROC

        CLRR    R2                ; R2 = k
        CLRR    R1                ; R1 = k(k - 1) / 2

@@loop  ADDR    R2,     R1        ; add R2 to R1
        INCR    R2                ; k++

        MOVR    R0,     R4        ; start with R4 = n

@@sub   SUBR    R1,     R4        ; subtract 2 * (k(k - 1) / 2) = k² - k
        SUBR    R1,     R4        ; from R4
        SUBR    R2,     R4        ; subtract k from R4
        BGT     @@sub             ; until R4 is less than or equal to 0

        BNEQ    @@next            ; did we reach exactly 0? ...

        MOVR    R2,     R3        ; ... yes: update R3

@@next  CMPR    R0,     R1        ; go on while R1 is less than or
        BLE     @@loop            ; equal to R0

        JR      R5                ; return

        ENDP

Output

This is the output for the following test cases:

4, 9, 12, 13, 108, 300, 800, 900

^{screenshot from jzIntv}

---

_{1. A CP-1610 opcode is encoded with a 10-bit value (0x000 to 0x3FF), known as a 'DECLE'.}

Whitespace, 71 bytes

[S S S N
_Push_0][S N
S _Duplicate_0][T   N
T   T   _STDIN_as_integer][T    T   T   _Retrieve_input][S N
S _n=Duplicate_input][N
S S N
_Create_Label_LOOP][S T S S T   N
_Copy_0-based_1st_input][S T    S S T   N
_Copy_0-based_1st_n][S N
S _Duplicate_n][T   S S N
_Multiply][T    S T T   _Modulo][N
T   S S N
_If_0_Jump_to_Label_PRINT_RESULT][S S S T   N
_Push_1][T  S S T   _Subtract][N
S N
N
_Jump_to_Label_LOOP][N
S S S N
_Create_Label_PRINT_RESULT][T   N
S T _Print_as_integer]

Letters S (space), T (tab), and N (new-line) added as highlighting only.
[..._some_action] added as explanation only.

Try it online (with raw spaces, tabs and new-lines only).

Port of @Sok's Pyth answer, so make sure to upvote him! Whitespace doesn't have decimals, so his/her approach is ideal for Whitespace since it doesn't use square-root nor regular division, but only integers.

Explanation in pseudo-code:

Integer n = STDIN as integer
Integer r = n
Start LOOP:
  Integer s = r * r
  If(n % s == 0):
    Jump to Label PRINT
  r = r - 1
  Go to next iteration of LOOP

Label PRINT:
  Print r as integer to STDOUT
  (implicitly stop the program with an error: no exit defined)

perl -MList::Util=max -pl, 34 bytes

$n=$_;$_=max grep!($n%$_**2),1..$n

Try it online!

This finds the largest square which properly divides the input number. Very inefficient as it tries all numbers from 1 up to the input.

Python 2, 37 bytes

n=i=input()
while n%i**2:i-=1
print i

Try it online!

38 bytes

lambda n:sum(x*x%n<1for x in range(n))

Try it online!

Based on the solution of my pronoun is monicareinstate, counting the number of solutions to $x^2 equiv 0 (mathbb{mod} n)$ for $x$ from $0$ to $n-1$.

Husk, 9 bytes

-1 byte thanks to Razetime!

▲fo¬%1m√Ḋ

Try it online! or Try all cases!

Commented

▲          -- the maximum of ...
 fo        -- ... filter on the next two functions composed
   ¬       --     - boolean negate / (== 0)
    %1     --     - modulo 1
      m√   -- ... map square root on ...
        Ḋ  -- ... the list of divisors

Pyth, 9 bytes

ef!%Q^T2S

Try it online!

ef!%Q^T2S   
        S   Create range from 1 to (implicit) input
 f          Filter keep from the above, as T, where:
     ^T2      Square T
   %Q         Mod the input with the above
  !           Logical NOT
e           Take the last (largest) element of the filtered list, implicit print

Java 10, 43 40 bytes

n->{for(var c=n++;c/--n%n>0;);return n;}

Inspired by @Arnauld's JavaScript answer, so make sure to upvote him!

Try it online.

Explanation:

n->{                // Method with double as both parameter and return-type
  for(var c=n       //  Create a copy `c` of the input `n`
             ++     //  Then increase `n` by 1
      ;             //  Continue looping as long as:
       c/--n        //   (decrease `n` by 1 first before every iteration with `--n`)
                    //   `c` divided by `n`
            %n>0;)  //   is NOT a multiply of `n` nor 0
  ;return n;}       //  After the loop: return the modified `n` as result

Pari/GP, 15 bytes

n->core(n,1)[2]

Yes, there is a build-in.

core(n,{flag=0}): unique squarefree integer d dividing n such that n/d is a square. If (optional) flag is non-null, output the two-component row vector [d,f], where d is the unique squarefree integer dividing n such that n/d=f^2 is a square.

Try it online!

R, 44 (crossed-out) 36 33 32 30 bytes

((n=scan()):1)[!n%%(n:1)^2][1]

Try it online!

Or, a completely-different 25 byte approach based on equivalence to 'number of solutions to x^2==0(mod n)' (as pointed-out by my pronoun is monicareinstate), but that wasn't my own idea and hence seems to me to be cheating: sum(!(1:(n=scan()))^2%%n)

Brachylog, ^{8 6} 5 bytes

-1 thanks to Unrelated String

f↔∋√ℕ

Try it online!

How it works

f↔∋√ℕ
    ℕ output is a natural number (≥0) that is
   √  the root of … (Brachylog gives the negative root first)
  ∋   an element …
f↔    in the reverse factors list (so search starts with bigger values)

Alternative version with prime factors, 10 bytes

{ḋp⊇~j×}ᵘ⌉

Try it online! or verify all test cases.

How it works

{ḋp⊇~j×}ᵘ⌉
         ⌉  take the maximum of …
{      }ᵘ   all unique …
      ×       multiplications of …   10
    ~j        halves of …            [2,5]
   ⊇          ordered subsets from … [2,5,2,5]
  p           the permutations of …  [2,5,2,5,3]
 ḋ            the prime factors      [2,2,3,5,5]

05AB1E, 9 6 5 bytes

LnIÖO

Try it online or verify all test cases.

Previous 9 6 bytes approach:

LR.ΔnÖ

-3 bytes thanks to @ovs.

Try it online or verify all test cases.

Explanation:

L       # Push a list in the range [1, (implicit) input]
 n      # Take the square of each value in the list
  IÖ    # Check which squares are divisible by the input (1 if truthy; 0 if falsey)
    O   # And sum those checks
        # (after which this sum is output implicitly as result)

L       # Push a list in the range [1, (implicit) input]
 R      # Reverse it to [input, 1]
  .Δ    # Find the first value in this list which is truthy for:
    n   #  Square the current value
     Ö  #  Check if the (implicit) input is evenly divisible by this square
        # (after which the found value is output implicitly as result)

JavaScript (ES6), 27 bytes

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

Try it online!

How?

We recursively look for the greatest $kle n$ such that $dfrac{n}{k}equiv 0pmod k$, which is guaranteed to be satisfied for $k=1$ in the worst case.

This is a more golf-friendly way of testing $dfrac{n}{k^2}equiv 0pmod 1$.