Suresh asked me to assemble my comments above into an answer, so here it is. I'm not really sure it's an answer to the original question, though, since it's not obvious how to make it polynomial time when the dimension of the input Euclidean space is non-constant. It does at least have the advantage of avoiding any problem with large $1/\epsilon$ as the original question asks, because it doesn't involve any approximation, and it looks polynomial for constant $d$.
Anyway: from integral geometry, there is a standard measure on sets of hyperplanes in $d$-dimensional Euclidean space that is invariant under Euclidean congruences. It has the property that the length of any bounded-length curve $C$ is proportional to the measure of the hyperplanes that cross $C$ (with multiplicity, meaning that if a hyperplane crosses $C$ twice then it contributes twice to the total measure of the hyperplanes crossing $C$). In particular if $C$ is a line segment then the multiplicity complication doesn't arise and we can normalize the measure on the hyperplanes crossing $C$ to be exactly the length of $C$. (The hyperplanes containing $C$ have measure zero, so don't worry about infinite multiplicity.)
Now, given a set of n points in d-dimensional space, make an $\ell_1$ coordinate for each of the partitions of the points into two subsets induced by a hyperplane that doesn't pass through any of the points. Give the points on one side of the partition coordinate value zero and the points on the other side of the partition coordinate value equal to the measure of the set of hyperplanes inducing that partition.
If $p$ and $q$ are any two of the $n$ points, let $K$ be the set of hyperplanes crossing line segment $pq$, and let $K_i$ be the subsets of $K$ formed by each possible hyperplane partition that has $p$ on one side and $q$ on the other. Then $K$ is the disjoint union of the $K_i$, and the coordinate differences between $p$ and $q$ are just the measures of the subsets $K_i$. Therefore, the $\ell_1$ distance between the coordinatizations of $p$ and $q$ (the sum of the measures of the $K_i$) is the measure of $K$, which is just the original $\ell_2$ distance between $p$ and $q$.
For computational geometers, an alternative description of the same construction may be helpful: use projective duality to transform the $n$ input points into $n$ hyperplanes, and separating hyperplanes into points. The integral geometry measure on sets of hyperplanes is then transformed into a more standard measure on sets of points, the distance between $p$ and $q$ dualizes to the measure of the double wedge between two hyperplanes, and the hyperplane arrangement partitions this double wedge into smaller cells. The coordinate value for a point is either the measure of one of the cells in the arrangement (if the dual hyperplane is below that coordinate's cell) or zero (if the dual hyperplane is above the cell). Therefore, the $\ell_1$ distance between $p$ and $q$ is just the sum of the measures of the cells in the double wedge, which is the same as the measure of the whole double wedge. This dual point of view also makes it easy to calculate the dimension of the embedding found in this way: it's just the number of cells in the hyperplane arrangement, which is $O(n^d)$, or more precisely at most $\sum_{i=0}^d{n\choose i}$.
So far, this gives a completely deterministic and exact embedding in $\ell_1^{O(n^d)}$. But we wanted a smaller dimension, $\ell_1^{n\choose 2}$. Here's where Luca's comment about Carathéodory's theorem comes in. The set of $\ell_1$-representable metrics forms a polyhedral cone in the $n\choose 2$-dimensional space of all functions from unordered pairs of points to real numbers, and the geometric argument above says that the Euclidean metric belongs to this cone. The points on the extreme rays of the cone are one-dimensional $\ell_1$ pseudometrics (in which the points are split into two sets, all distances within a single set are zero, and all distances across the split are equal), and Carathéodory says that any point within the cone (including the one we care about) can be represented as a convex combination of points on extreme rays whose number is at most the dimension of the ambient space, $n\choose 2$. But a convex combination of at most $n\choose 2$ one-dimensional $\ell_1$ metrics is a $\ell_1^{n\choose 2}$ metric.
Finally, how can we go about actually computing the $n\choose 2$-dimensional embedding? At this point we have not just a point in the $n\choose 2$-dimensional convex cone of $\ell_1$ metrics (the distance metric we started with), but we also have a set of $O(n^d)$ extreme points of the cone (corresponding to the partitions of the input into two subsets induced by hyperplanes) such that our metric is a convex combination of these extreme points — for small $d$, this is a big improvement over the $2^n-2$ extreme rays that the cone has overall. Now all we need to do is apply a greedy algorithm that gets rid of extreme points from our set, one by one, until only $n\choose 2$ of them are left. At each step, we need to maintain as an invariant that our metric is still inside the convex hull of the remaining extreme points, which is just a linear programming feasibility problem, and if we do this Carathéodory will ensure that there always remains a set of $n\choose 2$ extreme points whose convex hull contains the input metric.
