How to approximate $PV\int_0^\infty \frac{\tan x}{x}\text{d}x?$

The integral is Pi/2, a proof of which may be found on math.SE. Here's a numerical check, integrating $(\tan z)/z$ over parallel paths $z = x \pm a i$ and subtracting $(2 a)/(a^2 + x^2)$ to get NIntegrate to converge. The OP's integral is equal to the sum of the two below.

a = 50;
1/2 NIntegrate[
  Tan[x + a I]/(x + a I) + Tan[x - a I]/(x - a I) -
     (2 a)/(a^2 + x^2) //
    ComplexExpand // Simplify,
  {x, 0, Infinity}, AccuracyGoal -> 16]
1/2 Integrate[(2 a)/(a^2 + x^2), {x, 0, Infinity}, 
  Assumptions -> a > 0]
(*
  -7.83867*10^-47
  Pi/2
*)

Update

It was late and I was tired. A shortcut (due to the symmetry of Tan[x]/x) occurred to me this morning:

a = 50;
NIntegrate[Re[Tan[x]/x], {x, 0 + a I, Infinity}, WorkingPrecision -> 50]

(*  1.5707963267948966192313216916397514399782555923723  *)

You can evaluate principal-valued integrals using NIntegrate using Method->"PrincipalValue" but you must specify the singular point(s). In your case, there are an infinite number of them. For starters, integrating over just the first one at $\pi/2$, we can write:

NIntegrate[Tan[x]/x, {x, 0, Pi/2, 3}, 
 Method -> "PrincipalValue"]

Out[42]= 1.36501

Notice how I inserted the first singular point $\pi/2$ between the limits of integration. If we wanted to integrate over two of them, then:

In[44]:= NIntegrate[Tan[x]/x, {x, 0, Pi/2, 3 Pi/2, 6}, 
 Method -> "PrincipalValue"]

Out[44]= 1.46884

And if I wanted to integrate over 25 singular points, I would sequence-in a table of 25 singular points into the limits of integration and notice in the plot below it's convergence to $\pi/2$:

totalN = 25;
myTable = Table[
   {index, 
    NIntegrate[
     Tan[x]/x, {x, 0, 
      Sequence @@ Table[Pi/2 + n Pi, {n, 0, index}], (index + 1) Pi}, 
     Method -> "PrincipalValue"]},
   {index, 0, totalN}];
ListPlot[myTable, Joined -> True]

The method "PrincipalValue" does not work on an integral with infinitely many poles, since each pole must be specified. Confining to a finite region of integration is problematic on integrals like $\int_0^\infty (\tan x/x)\,dx$ that converge slowly. NIntegrate has built in the "ExtrapolatingOscillatory" strategy, but even if it would work together with the "PrincipalValue" strategy, you would still have the problem of specifying infinitely many poles. Instead, we can manually call NSum on the principal value integral over each period of Tan[x]. NSum will do the extrapolation for us. (Something similar was done here and is shown in the documenation.) The series still converges slowly, and the results are not as good as my other answer.

ClearAll[cc];
cc[n_Integer] := cc[n] =
   NIntegrate[Tan[x]/x, {x, 0 + n Pi, Pi/2 + n Pi, Pi + n Pi}, 
    WorkingPrecision -> 50, 
    Method -> {"PrincipalValue", Method -> "GaussKronrodRule"}];

NSum[cc[n], {n, 0, Infinity}, NSumTerms -> 50, 
  WorkingPrecision -> 50] // AbsoluteTiming

(*  {1.76767, 1.57067}  *)

---

Alternative

We can speed up the principal value integral by subtracting a null integral that smashes the pole (something similar was done here). The slow convergence persists, but the following is twice as fast as with the "PrincipalValue" method.

ClearAll[cc];
cc[n_Integer] := cc[n] = NIntegrate[
    Tan[x]/x - Tan[x]/(Pi/2 + n Pi), 
    {x, 0 + n Pi, Pi + n Pi}, 
    WorkingPrecision -> 50, Method -> "GaussKronrodRule"];

NSum[cc[n], {n, 0, Infinity}, NSumTerms -> 1000, 
  WorkingPrecision -> 50] // AbsoluteTiming

(*  {10.2142, 1.57078947}  *)

Last@% - Pi/2   (* error *)
(*  -6.85*10^-6  *)