JavaScript (Node.js), 144 138 125 74 73 70 bytes

f=(x,n=2,c=0)=>x%n?x-!c?f(x,n+1)/(n%4>2?n/=~c&1:n%4)**c:1:f(x/n,n,c+1)

Try it online!

-4 byte thanks @Arnauld!

Original approach, 125 bytes

a=>(F=(x,n=2)=>n*n>x?[x,0]:x%n?F(x,n+1):[n,...F(x/n,n)])(a).map(y=>r-y?(z*=[,1,.5,p%2?0:1/r][r%4]**p,r=y,p=1):p++,z=r=p=1)&&z

Try it online!

Inspired by the video Pi hiding in prime regularities by 3Blue1Brown.

For each prime factor $p^n$ in the factorization of the number, calculate $f(p^n)$:

• If $n$ is odd and $pequiv 3text{ (mod 4)}$ - $f(p^n)=0$. Because there is no place to go.


• If $n$ is even and $pequiv 3text{ (mod 4)}$ - $f(p^n)=frac{1}{p^n}$.


• If $p=2$ - $f(2^n)=frac{1}{2^n}$.


• If $pequiv 1text{ (mod 4)}$ - $f(p^n)=1$.

Multiply all those function values, there we are.

Update

Thanks to the effort of contributors from Math.SE, the algorithm is now backed by a proof

Clean, 189 185 172 171 bytes

import StdEnv

$n#r=[~n..n]

#p=[[x,y]\x<-r,y<-r|x^2+y^2==n]

=sum[1.0\_<-iter n(q=removeDup[k\[a,b]<-[[0,0]:p],[u,v]<-q,k<-[[a+u,b+v]]|all(e=n>=e&&e>0)k])p]/toReal(n^2)

Try it online!

Finds every position reachable in the n-side-length square cornered on the origin in the first quadrant, then divides by n^2 to get the portion of all cells reachable.

This works because:

• The entire reachable plane can be considered as overlapping copies of this n-side-length square, each cornered on a reachable point from the origin as if it were the origin.


• All movements come in groups of four with signs ++ +- -+ --, allowing the overlapping tiling to be extended through the other three quadrants by mirroring and rotation.

Mathematica, 80 bytes

d[n_]:=If[#=={},0,1/Det@LatticeReduce@#]&@Select[Tuples[Range[-n,n],2],#.#==n&];

This code is mostly reliant on a mathematical theorem. The basic idea is that the code asks for the density of a lattice given some generating set.

More precisely, we are given some collection of vectors - namely, those whose length squared is N - and asked to compute the density of the set of possible sums of these vectors, compared to all integer vectors. The math at play is that we can always find two vectors (and their opposite) that "generate" (i.e. whose sums are) the same set as the original collection. LatticeReduce does exactly that.

If you have just two vectors, you can imagine drawing an identical parallelogram centered at each reachable point, but whose edge lengths are the given vectors, such that the plane is completely tiled by these parallelograms. (Imagine, for instance, a lattice of "diamond" shapes for n=2). The area of each parallelogram is the determinant of the two generating vectors. The desired proportion of the plane is the reciprocal of this area, since each parallelogram has just one reachable point in it.

The code is a fairly straightforward implementation: Generate the vectors, use LatticeReduce, take the determinant, then take the reciprocal. (It can probably be better golfed, though)

Jelly,  25  24 bytes

ÆFµ%4,2CḄ:3+2Ịị,*/ʋ÷*/)P

A monadic link using the prime factor route.

Try it online!

How?

ÆFµ%4,2CḄ:3+2Ịị,*/ʋ÷*/)P - Link: integer, n               e.g. 11250

ÆF                       - prime factor, exponent pairs        [[2,1], [3,2], [5,4]]

  µ                   )  - for each pair [F,E]:

    4,2                  -   literal list [4,2]

   %                     -   modulo (vectorises)                [2,1]  [3,0]  [1,0]

       C                 -   complement (1-x)                  [-1,0] [-2,1]  [0,1]

        Ḅ                -   from base 2                         -2     -3      1      

         :3              -   integer divide by three             -1     -1      0

           +2            -   add two (call this v)                1      1      3

                  ʋ      -   last four links as a dyad, f(v, [F,E])

             Ị           -     insignificant? (abs(x)<=1 ? 1 : 0)   1      1      0

                */       -     reduce by exponentiation (i.e. F^E)  2      9     625

               ,         -     pair v with that                   [1,2]  [1,9]  [3,625]

              ị          -     left (Ị) index into right (that)     1      1     625

                    */   -   reduce by exponentiation (i.e. F^E)    2      9     625

                   ÷     -   divide                                1/2    1/9  625/625

                       P - product                                 1/18 = 0.05555555555555555

Previous 25 was:

ŒRp`²S⁼ɗƇ⁸+€`Ẏ;Ɗ%³QƊÐLL÷²

Full program brute forcer; possibly longer code than the prime factor route (I might attempt that later).

Try it online!

Starts by creating single moves as coordinates then repeatedly moves from all reached locations accumulating the results, taking modulo n of each coordinate (to restrict to an n by n quadrant) and keeping those which are distinct until a fixed point is reached; then finally divides the count by n^2

05AB1E, 27 26 25 bytes

ÓεNØ©<iozë®4%D≠iyÈ®ymz*]P

Port of @ShieruAsakoto's JavaScript answer, so make sure to upvote him as well!

Try it online or verify all test cases.

Explanation:

Ó                   # Get all prime exponent's of the (implicit) input's prime factorization

                    #  i.e. 6 → [1,1]      (6 → 2**1 * 3**1)

                    #  i.e. 18 → [1,2]     (18 → 2**1 * 3**2)

                    #  i.e. 20 → [2,0,1]   (20 → 2**2 * 3**0 * 5**1)

                    #  i.e. 25 → [0,0,2]   (25 → 2**0 * 3**0 * 5**2)

 ε                  # Map each value `n` to:

  NØ                #  Get the prime `p` at the map-index

                    #   i.e. map-index=0,1,2,3,4,5 → 2,3,5,7,11,13

    ©               #  Store it in the register (without popping)

     <i             #  If `p` is exactly 2:

       oz           #   Calculate 1/(2**`n`)

                    #    i.e. `n`=0,1,2 → 1,0.5,0.25

      ë             #  Else:

       ®4%          #   Calculate `p` modulo-4

                    #    i.e. `p`=3,5,7,11,13 → 3,1,3,3,1

          D         #   Duplicate the result (the 1 if the following check is falsey)

           ≠i       #   If `p` modulo-4 is NOT 1 (in which case it is 3):

             yÈ     #    Check if `n` is even (1 if truthy; 0 if falsey)

                    #     i.e. `n`=0,1,2,3,4 → 1,0,1,0,1

             ®ymz   #    Calculate 1/(`p`**`n`)

                    #     i.e. `p`=3 & `n`=2 → 0.1111111111111111 (1/9)

                    #     i.e. `p`=7 & `n`=1 → 0.14285714285714285 (1/7)

              *     #    Multiply both with each other

                    #     i.e. 1 * 0.1111111111111111 → 0.1111111111111111

                    #     i.e. 0 * 0.14285714285714285 → 0

 ]                  # Close both if-statements and the map

                    #  i.e. [1,1] → [0.5,0.0]

                    #  i.e. [1,2] → [0.5,0.1111111111111111]

                    #  i.e. [2,0,1] → [0.25,1.0,1]

                    #  i.e. [0,0,2] → [1.0,1.0,1]

  P                 # Take the product of all mapped values

                    #  i.e. [0.5,0.0] → 0.0

                    #  i.e. [0.5,0.1111111111111111] → 0.05555555555555555

                    #  i.e. [0.25,1.0,1] → 0.25

                    #  i.e. [1.0,1.0,1] → 1.0

                    # (and output implicitly as result)

Retina 0.8.2, 126 bytes

.+

$*

+`^(11+)(1)+b

$1;1$#2$*

b1(1111)+b

1

b(1(11)+);1b

$1_$1

.*b(1(11)+)b.*

0

_

;

+`G1(?=.*;(1+))|;1+$

$1

11+

1/$.&

Try it online! Link includes test cases. Uses the prime factorisation. Explanation:

.+

$*

Convert to unary.

+`^(11+)(1)+b

$1;1$#2$*

Get the prime factors.

b1(1111)+b

1

Delete factors of the form 4k+1.

b(1(11)+);1b

$1_$1

Temporarily mark any duplicate pairs of remaining odd factors (i.e. of the form 4k+3).

.*b(1(11)+)b.*

0

If there are any unmatched odd factors remaining then the result is zero.

_

;

Remove the markings.

+`G1(?=.*;(1+))|;1+$

$1

Multiply any remaining factors back together. (Taken from the Unary arithmetic tutorial.)

11+

1/$.&

Compute the reciprocal as a decimal fraction.