Vyxal, 3 bytes

εġ‹

Try it Online!

Port of Jonathan Allan's Jelly answer.

ε   # Absolute difference - [|x1-x2|,|y1-y2|]
 ġ  # gcd of that list
  ‹ # decrement

Jelly, 4 bytes

ạg/’

A dyadic Link that accepts a point as a pair on each side and yields the count.

Try it online!

How?

Implements the provided greatest common divisor formula.

ạg/’ - Link: pair of integers [x1, x2]; pair of integers [y1, y2]
ạ    - absolute difference (vectorises) -> [|x1-y1|, |x2-y2|]
  /  - reduce by:
 g   -   greatest common divisor -> GCD(|x1-y1|, |x2-y2|)
   ’ - decrement -> GCD(|x1-y1|, |x2-y2|) - 1

Nekomata, 4 bytes

≈đG←

Attempt This Online!

Takes input as two pairs of numbers, e.g. [-8,5] [0,5].

≈đG←
≈       Absolute difference (vectorized)
 đ      Unpair; get the two elements of a pair
  G     GCD
   ←    Decrement

I wonder how an approach that defines the solution to be an integer point on the line, and using -n for the final result. (I haven't really tried to learn Nekomata yet and somehow doubt it would be shorter for this problem, but seems interesting) – noodle man

The shortest I can get using this approach is 10 bytes:

Nekomata + -n, 10 bytes

≈:Ṁᵉ{~Z*}¦

Attempt This Online!

≈:Ṁᵉ{~Z*}¦      Take [-8,5] [0,5] as an example
≈           Absolute difference (vectorized)
                [-8,5] [0,5] -> [8,0]
 :          Duplicate
                [8,0] -> [8,0] [8,0]
  Ṁ         Maximum
                [8,0] [8,0] -> [8,0] 8
   ᵉ{       Apply the following block and then push the original top of stack
     ~Z       Choose any integer in [1,n)
                [8,0] 8 -> [8,0] 1 or [8,0] 2 or ... or [8,0] 7
       *      Multiply
                [8,0] 1 -> [8,0]
                [8,0] 2 -> [16,0]
                ...
                [8,0] 7 -> [56,0]
        }   End the block (ᵉ pushes the original top of stack)
                [8,0] -> [8,0] 8
                [16,0] -> [16,0] 8
                ...
                [56,0] -> [56,0] 8
         ¦  Divide and check if the result is an integer
                [8,0] 8 -> [1,0]
                [16,0] 8 -> [2,0]
                ...
                [56,0] 8 -> [7,0]

-n counts the number of solutions, so the output is 7.

APL(Dyalog Unicode), 7 bytes ^SBCS

¯1+∨/⍤-

Try it on APLgolf!

A tacit function which takes vectors on the left and right and returns an integer. I believe this solution also works in any dimension, not just 2D. It takes the difference of the points, GCD all of the differences, and adds -1.

Charcoal, 28 bytes

ＮθＮη≧⁻Ｎθ≔∨⁻ηＮθηＩＬΦ↔η∧ι¬﹪×ιθη

Try it online! Link is to verbose version of code. Explanation: No vectorised difference or gcd, so I have to do things the hard way.

ＮθＮη

Input the first co-ordinate.

≧⁻Ｎθ≔∨⁻ηＮθη

Subtract the second co-ordinate, but if it is horizontal then make it diagonal, which has the same number of intermediate points.

ＩＬΦ↔θ∧ι¬﹪×ιηθ

Generate a list of intermediate y-coordinates and see how many x-coordinates are integers.

Uiua, 12 bytes ^SBCS

-1;⍢⊃◿∘±°⊟⌵-

Try it!

-1;⍢⊃◿∘±°⊟⌵-
           -  # difference
          ⌵   # absolute value
        °⊟    # uncouple pair to stack
  ;⍢⊃◿∘±      # GCD
-1            # decrement

Retina 0.8.2, 97 bytes

,(.+),(.+),
,$2¶$1,
%O`[^,]+
\d+
$*
-(1+),(1+)
$1$2
(1+),-?\1|[^1¶]

mO^`^.*
^(1(1*))\1*¶\1*$
$.2

Try it online! Link includes test cases. Explanation:

,(.+),(.+),
,$2¶$1,

Exchange the first y-coordinate with the second x-coordinate and split the coordinates onto separate lines.

%O`[^,]+

Ensure that if either coordinate is negative then the first one is.

\d+
$*

Convert to unary.

-(1+),(1+)
$1$2

If the co-ordinates have different signs then add their absolute values together.

(1+),-?\1|[^1¶]

Otherwise take the absolute difference.

mO^`^.*

Sort them descending in case the first difference is zero.

^(1(1*))\1*¶\1*$
$.2

Calculate the decremented GCD.

R, 62 bytes

\(x,y)max((y=1:max(z<-abs(x-y)))[!z[1]%%y&!z[2]%%y|!all(z)])-1

Attempt This Online!

Base R has no built in for GCD.

JavaScript (Node.js), 50 bytes

(p,q,r,s)=>(g=r=>s?g(s,s=r%s):r*r)(r-p,s-=q)**.5-1

Try it online!

Basically based on the formula in the question. Calculate

$$ \sqrt{\left(\text{gcd}(x_2-x_1,y_2-y_1)\right)^2}-1 $$

where \$\text{gcd}\$ is defined as

$$ \text{gcd}(x,y)=\begin{cases} x & y=0 \\ \text{gcd}(y,x \bmod y) & \text{otherwise} \end{cases} $$

APL(NARS), 9 chars

¯1+∨/∣⎕-⎕

test&how use:

      ¯1+∨/∣⎕-⎕
⎕:
      ¯8 5
⎕:
      0 5
7
~

      ¯1+∨/∣⎕-⎕
⎕:
      5 10
⎕:
      10 5
4
~

05AB1E, 4 bytes

-¿Ä<

Try it online!

Uiua, 6 bytes

-1/∨⌵-

Try it! (Takes input as x1_y1 x2_y2)
Take the ⌵- absolute difference of the points, / reduce the difference with ∨ gcd, and -1 decrement.

C (gcc), 104 bytes

-5 bytes, thanks to @ceilingcat

g(a,b){return b?g(b,a%b):abs(a);}main(x,y,a,b){scanf("%d%d%d%d",&x,&y,&a,&b);printf("%u",g(x-a,y-b)-1);}

Try it online!