As others have indicated one has to distinguish a metric $(x,y)\mapsto d(x,y)$ which measures distances between points $x$, $y$ in a space $X$ and a Riemannian metric which tells us how lengths of curves $\gamma$ in a manifold $X$ should be computed.
In the complex plane we have the usual euclidean metric $d(z_1,z_2):=|z_2-z_1|$ and at the infinitesimal scale the Riemannian metric
$$ds^2:=|dz|^2=dx^2+dy^2\ .$$
The latter formula says that the length of an arbitrary curve
$$\gamma: \quad t\mapsto\bigl(x(t),y(t)\bigr)\qquad(a\leq t\leq b)$$
should be computed as
$$L(\gamma)=\int_\gamma ds=\int_a^b\sqrt{x'^2(t)+y'^2(t)}\ dt = \int_a^b|z'(t)|\ dt\ .$$
This formula implies that the length of a segment $\sigma$ connecting two points $z_1$ and $z_2$ is just $|z_2-z_1|$.
On the "Poincare disc" $P$ we are given a priori only a Riemannian metric
$$ds:= {|dz|\over 1-|z|^2}$$
(resp., $ds^2=\ldots$). This metric allows us to compute the lengths of arbitrary curves in $P$:
$$L(\gamma)=\int_a^b{|z'(t)|\over1-|z(t)|^2}\ dt\ .$$ The particular definition of $ds$ is chosen such that this hyperbolic length is invariant under arbitrary conformal movements of $P$ and the curves therein.
A posteriori one can define a metric $d(\cdot,\cdot)$ on $P$ by letting the distance $d(z_1,z_2)$ between two points $z_1$, $z_2\in P$ be the hyperbolic length of the shortest curve $\gamma$ connecting $z_1$ and $z_2$. The actual carrying out of this idea shows that $d(z_1,z_2)$ can be written as an elementary function (using ${\rm artanh}$, etc.) in terms of $z_1$ and $z_2$.
The problem is that things like scalar multiplication and vector addition/subtraction are not well defined in hyperbolic geometry -- you could invent some definition, but then they will not have the same properties as in Euclidean, so it is probably better to avoid such operations, or call them something else. What do you need "scalar multiplication" for?
I sometimes use scalar multiplication and vector addition/subtraction which act in the whole Minkowski space (not restricted to the hyperboloid), i.e., $(x_a,y_a,z_a)+(x_b,y_b,z_b) = (x_a+x_b, y_a+y_b,z_a+z_b)$. These operations are well defined and they have some uses, for example, to find a midpoint of $(A,B)$, you add the coordinates and project the result back to the hyperboloid.
I think the rounding errors are unavoidable if you are working with hyperbolic geometry extensively. The main reason is that the hyperbolic plane grows exponentially -- a circle of radius 200 will have area of roughly $e^{200}$, so even if you are using the Minkowski hyperboloid model with three 64-bit coordinates, you have only $2^{192}$ representable points. If you have, say, two creatures starting in a small distance $\epsilon$ and going in a roughly parallel direction, the distance between them will grow to roughly $\epsilon e^d$ after each of them moves $d$ units. Likewise, the $\epsilon$ could be caused by floating point errors, so floating point errors will blow up in cases similar to this. In HyperRogue, the matrices such as the view transformation naturally accumulate floating point errors, eventually the errors build up so much that the matrix does not make any sense anymore. To fix this, the transformations are "fixed" from time to time -- this operation transforms the given matrix to a nearby isometry, losing any built up errors.
You can see our implementation of hyperbolic geometry here: GitHub (it is a bit more complex because it works with Euclidean and spherical geometry too; also, points outside of the hyperboloid are used to represent the 3D hyperbolic space in a rather weird way). Not everything is optimal in terms of computational cost, but it works good enough for HyperRogue. (The main ideas are roughly described here.) For things you would use vector addition/subtraction in Euclidean geometry, HyperRogue generally uses functions "find a matrix which rotates the given point to the x axis" or "find a matrix which translates along the x axis to move the given point to 0", which are indeed rather inefficient. But IIRC the basic things like tiling generation do not need such operations, these are more useful for complex things like movement animations. No operation similar to scalar multiplication is used -- instead there is xpush(d) which translates $d$ unit along the $X$ axis, so instead of trying to multiply a vector by a scalar, you would say what direction and distance you want, and apply the rotation and xpush(d).
Best Answer
I assume you have read this: http://www.roguetemple.com/z/hyper/dev.php
You say you have computed the edge length -- you could base your implementation on this, but it is easier when you compute the distance between two centers of heptagons. You should get the following result: $d=1.090550...$
The center is $h_0 = \left(\begin{array}{c}0\\0\\1\end{array}\right)$ and it is also the center of the central heptagon.
To compute the coordinate of $h$ translated by $x$ units, compute $T_xh$, where $T_d$ is the following matrix: $\left(\begin{array}{c}\cosh x&0&\sinh x\\0&1&0\\\sinh x&0&\cosh x\end{array}\right)$ (note the similarity to the rotation of the sphere). So you can compute the center of the heptagon to the right from the central heptagon by $T_dh_0$.
To find the centers of other adjacent heptagons, you just need to multiply by the rotation matrix. Rotation by $\alpha$ is just like in Euclidean geometry: $R_\alpha = \left(\begin{array}{c}\cos \alpha&\sin \alpha&0\\-\sin\alpha&\cos\alpha&0\\0&0&1\end{array}\right)$. Thus the centers of seven adjacent heptagons are $R_{k\alpha}T_dh_0$, where $\alpha = 2\pi/7$ and $k=0..6$.
To get to the center of a heptagon in distance 2, you need to: rotate by angle $\alpha k$, move by $d$, rotate by $\pi$, rotate by angle $\alpha l$, move by $d$. Thus the formula for this point is: $R_{\alpha k} T_d R_\pi R_{\alpha l} T_d h_0$ By taking $k,l=0..6$ you get all heptagons in distance at most 2, some multiple times. (Since you have the graph, I assume you know how to clean this up.) Similarly you find the centers of heptagons in large distances.
The easiest way to compute the vertices is to take the sum $s$ of Minkowski hyperboloid of the three adjacent heptagon centers, and divide it by $k$ such that $s/k$ lies on the Minkowski hyperboloid.
In general the Minkowski hyperboloid is good because the formulas for it are directly analogous to the formulas for computations on the sphere (in three-dimensional coordinates, i.e., $x,y,z$ in the $\mathbb{R}^3$), which are rather easy to find out. They just have some signs changed, and cos/sin replaced with cosh/sinh if they correspond to distances. You could think about how to implement a dodecahedron -- the formulas for spherical geometry (in x,y,z coordinates) are in most cases the same as for hyperbolic geometry (in Minkowski hyperboloid model),