Linear combination of discrete variables $T_i$ with $P(T_i=1)=P(T_i=-1)=1/2$ Let $T_1,...,T_n$ be iid with a Rademacher distribution; i.e., $P(T_i=1)=P(T_i=-1)=1/2$; and let $w = (w_1,...,w_n) \in \mathbb{R}^n$ without further constraints on $w$.
Is there a way to compute $P(T_1\sum_{i=1}^n T_i w_i \geq 0)$ without stepping through all $2^n$ possible outcomes for the $T$'s?
 A: Yes.
Let's begin by simplifying the question.  The event
$$T_1\sum_{i=1}^n T_i w_i \geq 0$$
is the union of the disjoint events (determined by $T_1=\pm 1$)
$$w_1 + \sum_{i=2}^n T_i w_i \geq 0,\ -w_1 + \sum_{i=2}^n T_i w_i \geq 0.$$
After simple algebraic manipulation, and using the fact that each $T_i$ has the same distribution as $-T_i$ (it is "symmetric"), these events are both equivalent to those of the form
$$\sum_{i=2}^n T_i w_i \leq w$$
for $w=\pm w_1$.  Again exploiting the symmetries of the $T_i$ (to make all the $w_i$ non-negative) and re-ordering them, this can be written
$$\sum_{i=2}^n T_i |w_i| \leq w$$
where $|w_2| \ge |w_3| \ge \cdots \ge |w_n|$.
We need to compute the chances of two such events and wish to do so efficiently.  I will assume that the chance needs to be computed to high accuracy, so that approximations would be unacceptable.  (For many configurations of the $|w_i|$, Normal approximations to the distributions might otherwise work well.  Deeper analysis can be carried out using the characteristic function $\phi$ of $\sum_{i=2}^n T_i |w_i|$, which by definition is $$\phi(t) = \prod_{i=2}^n\left(\frac{\exp(-j |w_i|t)}{2} + \frac{\exp(j |w_i|t)}{2}\right)=\prod_{i=2}^n\cos(w_it);$$ $j = \sqrt{-1}$.)

From now on, reduce $n$ by $1$, start the indexing with $i=1$, and assume all the $w_i$ are positive.  (Obviously, any zero values can be ignored.)  Let $1 \le k \lt n$.  Consider the distribution of $X = \sum_{i=1}^n T_i w_i$ conditional on the first $k$ values of $T_i$:
$$\Pr(X \le w\,|\, T_1, \ldots, T_k) = {\Pr}_{(T_{k+1},\ldots,T_n)}\left(\sum_{i=k+1}^n T_i w_i \le w - \sum_{i=1}^k T_i w_i = w_0\right).$$
Since $|T_i| \le 1$, the left hand sum is bounded:
$$-u_{k+1} = -(w_{k+1} + \cdots + w_n) \le \sum_{i=k+1}^n T_i w_i \le w_{k+1} + \cdots + w_n = u_{k+1}.$$
Let $u_{k+1}$ be the right hand side: it is a cumulative sum of the $w_i$, accumulated from the rightmost (lowest) values.  Obviously now if $u_{k+1} \le w_0$ then $X \le w$ is certain; and if  $-u_{k+1} \gt w_0$, then $X \le w$ has zero probability.  This leads to a simple branch and bound algorithm, because we needn't search any further to assess the distribution.  Whereas exhaustive enumeration would have had to examine $2^{n-k}$ possible cases, we have made a determination of their contribution to the distribution in $O(1)$ time.
There is not much hope that this improvement will lead to a worst-case algorithm that is better than $O(2^n)$ in performance (although I think it can be reduced to $O(2^{n/2})$ which--although much better--is still non-polynomial).  However, the improvement is good enough to be worth considering, especially for moderate values of $n$ where exhaustive enumeration starts becoming impracticable (somewhere above $20$ and certainly above $40$).  Let us therefore turn to benchmarking the algorithm.  How efficient is it?
The worst case is when most of the $w_i$ have comparable sizes, for then the branch-and-bound heuristic rarely accomplishes anything.  Fortunately, this is exactly the situation where the Central Limit Theorem can supply excellent approximations!  It will therefore be perhaps even more interesting to explore situations where the sizes of the $w_i$ are highly spread out.
To this end, I created four datasets with $n=15$ from four distributions with radically different shapes: Beta$(1/5,1/5)$ (sharply bimodal), Exponential (high mode near zero), Uniform, and Gamma$(20)$ (nearly Normal, with values closely arranged near $20$).  Each dataset was normalized to a maximum of $1$.  Using the branch-and-bound method, $\Pr(X \le w)$ was computed for $19$ values of $w$ ranging from $0$ up to $2\sqrt{n}$.  (Negative values of $w$ need not be shown because the distributions of $X$ are nearly symmetric.)  The figure displays the results, graphing the probabilities (the empirical cumulative distribution function of $X$) in the top row and the efficiencies (on log-linear scales) in the bottom row.  The "efficiency" is the ratio of the number of configurations needed for exhaustive enumeration ($2^n$) to the number of configurations considered by the branch-and-bound algorithm.

For the Beta distribution, whose strong bimodality essentially limits the calculations to the higher half of the data, efficiencies are greatest.  As expected, they are least for the Gamma distribution, whose values are all comparable.  Nevertheless, even in this difficult case, the efficiencies all exceed $3$.
The smallest efficiencies are usually for $w=0$.  Eventually efficiency will increase as $|w|$ increases.  Collectively these results are strong evidence that the approach described here not only is computationally more efficient than exhaustive enumeration, it tends to be much more efficient.

The implementation of the algorithm is straightforward and simple. R code follows.  It includes a section that tests the algorithm (by comparing its output to an exhaustive enumeration) and another section to reproduce the figure.
#
# The algorithm.
#
f <- function(w0, w) {
  w <- sort(abs(w), decreasing=TRUE)
  w.sum <- c(rev(cumsum(rev(w)))[-1], 0)
  count <- 0                              # Counts calls to f()
  f <- function(u, u.sum, w0) {
    count <<- count + 1
    y <- sapply(w0 + c(-1,1)*u[1], function(w1) {
      if (w1 < -u.sum[1]) return(0)
      if (w1 >= u.sum[1]) return(1)
      return(f(u[-1], u.sum[-1], w1))
    })
    return(mean(y)) # (This could easily be changed to accommodate 
                    # other probabilities for the T[i].)
  }
  list(Value=f(w, w.sum, w0), Count=count)
}
#
# Test with complete enumeration.  The plot pairs should exactly coincide.
#
binary <- function(a, b, zero=-1, one=1) rep(c(rep(zero, 2^a), rep(one, 2^a)), 2^b)
n <- 9
b <- sapply(1:n, function(a) binary(n-a, a-1))
par(mfrow=c(2,2))
set.seed(17)
for (i in 1:4) {
  w <- rexp(n)
  x <- b %*% w
  y <- sapply(x, function(w0) f(w0, w)$Value) #$
  plot(ecdf(x))
  points(x, y, pch=16, cex=1/2, col="Red")
}
#
# Explore the efficiencies actually achieved.
#
n <- 15
qb <- function(q) qbeta(q, 1/5, 1/5); sigma <- sqrt(1/(1+2*1/5))/2
qg <- function(q) qgamma(q, 20)
dist <- list(Gamma=qg, Normal=qnorm, Uniform=qunif, Exponential=qexp, Beta=qb)
par(mfcol=c(2,4))
for (s in c("Beta", "Exponential", "Uniform", "Gamma")) {
  d <- dist[[s]]
  w <- d(((1:n)-1/2)/n)
  w <- w / max(abs(w))
  print(c(system.time({
    y <- sapply(x <- seq(0, 2*sqrt(n), length.out=19), 
                function(w0) {x <- f(w0, w); c(x$Value, x$Count)})
  })["elapsed"], count=sum(y[2, ])))

  plot(x, y[1, ], ylim=c(1/2, 1), type="b", ylab="Probability", main=s)
  plot(x, 2^n/y[2, ], ylim=c(1, 2^n), type="b", log="y", ylab="Efficiency")
}

