How to calculate the square root of a number? By searching I found few methods but all of them involve guessing which is not what I want.
I need to know how to calculate the square root using a formula or something. In other words how does the calculator calculate the square root?

The easiest way to find $\sqrt[n]{a}$ for integer $n$ and $a>0$ efficiently is to use the Newton-Raphson approximation to invert the function $f : x \mapsto x^n - a$. But one must be careful with choosing the right starting point, so that the iteration will converge quadratically. Quadratic convergence means that at each step the error becomes approximately a constant times its square, which is equivalent to the error being proportional to $c^{2^k}$ after $k$ steps, for some $c \in (0,1)$
Let $x_0$ be such that $x_0 \in \sqrt[n]{a}[1,1+\frac{1}{4n})$
For each natural $k$ from $0$ to $\infty$:
  Let $x_{k+1} = x_k - \dfrac{f(x_k)}{f'(x_k)} = x_k - \dfrac{{x_k}^n-a}{n{x_k}^{n-1}} = \dfrac{(n-1){x_k}^n-a}{n{x_k}^{n-1}}$
Then $( x_k : k \in \mathbb{N} )$ converges quadratically to $\sqrt[n]{a}$ uniformly for all $a>0$
General Case
For any real function $f$ such that $f(r) = 0$ and $f' \ne 0$ and $f''$ exists and $\left|\frac{f''}{2f'(r)}\right| \le m$ for some $m$:
  Let $a = f'(r) \ne 0$
  Then $f(r+d) = a d + g(d) d^2$ for any $d$ for some function $g$ such that:
    $g(d) \in a [-m,m]$ for any $d$
  Also $f'(r+d) = a + h(d) d$ for any $d$ for some function $h$ such that:
    $h(d) \in a [-m,m]$ for any $d$
  Let $( x_k : k \in \mathbb{N} )$ be such that:
    $x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$ for any natural $k$
    $|x_0-r| \le \frac{1}{6m}$
  For each natural $k$ from $0$ to $\infty$:
    $x_k = r + d_k$ for some $d_k$
    $|d_k| \le |d_0| \le \frac{1}{6m}$ by invariance
    $x_{k+1} = (r+d_k) - \dfrac{ad_k+g(d_k){d_k}^2}{a+h(d_k){d_k}} \in (r+d_k) - \dfrac{d_k+[-m,m]{d_k}^2}{1+[-m,m]{d_k}}$
    Thus $d_{k+1} \in d_k - (d_k+[-m,m]{d_k}^2) (1-[-m,m]{d_k}+[0,2]([-m,m]{d_k})^2)$ because:
      $\frac{1}{1+t} \in 1-t+[0,2]t^2$ for any $t \ge -\frac{1}{2}$
    Thus $d_{k+1} \in d_k - (d_k+[-m,m]{d_k}^2) (1+[-m,m]{d_k}+\frac{1}{3}[-m,m]d_k) \\ \quad \subseteq \frac{7}{3}[-m,m]{d_k}^2 + \frac{4}{3}[-m,m]^2{d_k}^3 \subseteq \frac{7}{3}[-m,m]{d_k}^2 + \frac{7}{18}[-m,m]{d_k}^2 \\ \quad \subset 3[-m,m]{d_k}^2 \subset [-1,1]d_k$
    Thus the invariance is preserved
    Also $3 m |d_{k+1}| < ( 3 m |d_k| )^2$
  Therefore $3 m |d_k| < ( 3 m |d_0| )^{2^k} \le 2^{-2^k}$ for any natural $k$
  Thus $x_k \to r$ quadratically as $k \to \infty$
Notes
In the case of finding $r = \sqrt[n]{a}$, the function $f : x \mapsto x^n - a$ has $\frac{f''}{2f'(r)}$ being $x \mapsto \frac{(n-1)x^{n-2}}{2r^{n-1}}$ which is bounded on $r[1,1+\frac{1}{4n})$ by $m = \frac{2n}{3r}$ because $\frac{n}{2r} (\frac{x}{r})^{n-2} \le \frac{n}{2r} (1+\frac{1}{4n})^n < \frac{n}{2r} e^{1/4} < m$. Thus $|x_0-r| < \frac{r}{4n} = \frac{1}{6m}$.
The procedure to find $x_0$ for efficient arbitrary precision arithmetic can be as follows:
  Find the minimum integer $d$ such that $(2^d)^n \ge a$
  Binary search on $[2^{d-1},2^d]$ to find $r$ until within an error of $\frac{2^{d-1}}{4n}$
  Return the upper bound when the upper and lower bounds are within the error margin
  The upper bound is between $r$ and $r+\frac{2^{d-1}}{4n} < r+\frac{r}{4n}$
