You need to solve polynomial equation in the bottom of the picture:

Keep in mind that some of these values are vectors so you need to use dot product in some places.

Also:

1. If player is on a platform then it isnt affected by gravity, so gravity shouldnt be included in player acceleration term

2. It only accurate if player acceleration/velocity is constants

I am not sure if this equation can be solved analitically... If it cant than it is quiet a messy task to solve. You first need to find a first root of a polynomial numerically, than divide initial polynom by (x - first_root) using Horner's method and obtain cubic polynom. Then you need to solve cubic polynomial equation (either analitically or numerically).

I'd consider what Alberth has written before implementing all of this. Are you really need super accurate enemies? If you drop down accelerations (a = 0) you'll end up with simple quadratic equation that can be easily solved. Acquired solution can be used as initial guess (+- some coefficient * distance to player to take gravity into account). If an enemy undershots/overshots next shot is corrected appropriately.

I wrote an article about how to solve this type of problem: http://www.gamedev.net/page/resources/_/technical/math-and-physics/leading-the-target-r4223

I'll try to write a piece of code specific to your situation when I have a chance (probably late tonight).

Some bullets are quite slow and fast shooting interval, so correcting appropriately won't really work, plus it's gonna be a real mess. It would be cool to get the most accurate direction and then later I could just add some randomness based on bot's noobness level.

What I suggest you is to first check whether simplified solutions suit your needs. Either solution with no acceleration at all or acceleration is taken into account by heuristic:

1. find t from: (player_velocity^2 - bullet_speed^2) * t^2 + 2 * (player_pos - shooter_pos) * player_velocity * t + (player_pos - shooter_pos)^2 = 0

2. find direction from: direction = normalize((player_velocity * t + player_pos) - shooter_pos)

3. (optionally) bullet_velocity = direction * speed + gravity_heuristic(t), where gravity_heuristic: some_coefficient * t

The thing is that the player probably will be moving and trying to dodge bullets. If some bullets will go off target he may not notice this. It may save you development time that you'd spend on the feature the player wont notice. Even if you personally would not like the result I suggest you to test it on some user group.

If after all you'll decide to go to the end. Well, there is no simple solution for this problem. You'll have to use iterative method. You'll need to determine either you undershoot or overshoot and correct bullet_velocity accordingly. But if you need accurate solution you'll have to make all the calculations before shooting real bullet.

You have the following equation vector equation that should be rewritten as two separate equations for x and y (sorry... :) ):

bullet_velocity equals bullet_speed * normalize(K_dist), a is deflection angle. Initially, a = 0.

Then you need to solve this vector equations and find tx, ty. If tx == ty you have found the solution.

If not:

1. Either there is no solution at all.

2. Or you need to "rotate gun" which is done by changing value of a (decrementing or incrementing).

There is a way to determine whether there are solutions and whether you need to decrement or increment value of a, but it is really late night (actually, it is early morning) so I leave it for now. I suggest you to try the following: make a small increment, find positive tx, ty => dPos = |tx - ty|*, make small decrement, find positive tx, ty dNeg => |tx - ty|. If dPos < dNeg use increment operation and decrementing otherwise. Make small steps in chosen direction to minimize |tx - ty|.

* if there were no positive tx, ty or no roots at all, take bigger step. If it wont help repeat several times and if it still fail then there are might be no solutions at all.

P.S.: It is really late... maybe I am missing something... maybe there is simpler solution or errors in mine.

I wrote an article about how to solve this type of problem: http://www.gamedev.net/page/resources/_/technical/math-and-physics/leading-the-target-r4223

I'll try to write a piece of code specific to your situation when I have a chance (probably late tonight).

Acceleration mess everything

[EDIT: I added a comment indicating where you should start reading the code. Everything before is library-style code that doesn't add much to understanding how I solve the problem.]

Well, here's some code that seems to do the job. I wrote the whole thing today and I apologize for its many shortcomings, but it should serve the purpose of illustrating how your problem can be solved.

The program is rather long, because it contains a set of functions to handle polynomials, a set of functions to handle vectors and a few extra things to be able to handle "vectors" of polynomials (not really vectors, because polynomials don't form a field, so instead of a vector space we are dealing with a module; but that probably only bothers people who know entirely too much math). By using Vector<Polynomial> I can describe the equation that gives me the time to interception symbolically, which should make the code more understandable.

Enjoy!

#include <iostream>
#include <vector>
#include <initializer_list>
#include <cmath>
#include <cassert>
#include <limits>

struct Polynomial {
  std::vector<float> coefs;
  
  float operator[](size_t i) const {
    return i < coefs.size() ? coefs[i] : 0.0f;
  }
  
  template <typename It>
  Polynomial(It begin, It end) : coefs(begin, end) {
  }
  
  explicit Polynomial(size_t size) : coefs(size) {
  }
  
  Polynomial(std::initializer_list<float> const &list) : coefs(list) {
  }
  
  float eval(float x) const {
    float result = 0.0f;
    float x_pow = 1.0f;
    for (float c : coefs) {
      result += c * x_pow;
      x_pow *= x;
    }
    return result;
  }
  
  float eval_derivative(float x) const {
    float result = 0.0f;
    float x_pow = 1.0f;
    for (size_t i = 1; i < coefs.size(); ++i) {
      result += i * coefs[i] * x_pow;
      x_pow *= x;
    }
    return result;
  }
  
  float newton_raphson(float x) const {
    for (int i = 0; i < 10; ++i)
      x -= eval(x) / eval_derivative(x);
    return x;
  }
  
  bool find_root(float &x) const {
    for (int i = 0; i < 10; ++i) {
      float initial_guess = i;
      x = newton_raphson(initial_guess);
      if (std::abs(eval(x)) < 1.0e-4f) // I hate having a magic constant here, but I want to keep things simple
	return true;
    }
    return false;
  }
  
  Polynomial divide_by_x_minus(float r) const {
    assert(coefs.size() > 0);
    Polynomial result(++coefs.begin(), coefs.end());
    float carry = 0.0f;
    for (size_t i = coefs.size() - 1; i > 0; --i) {
      result.coefs[i - 1] = coefs[i] + carry;
      carry = result.coefs[i - 1] * r;
    }
    return result;
  }    
};

Polynomial operator*(Polynomial const &p, float scalar) {
  Polynomial result = p;
  for (size_t i = 0; i < p.coefs.size(); ++i)
    result.coefs[i] = p.coefs[i] * scalar;
  return result;
}

Polynomial operator*(float scalar, Polynomial const &p) {
  return p * scalar;
}

Polynomial operator+(Polynomial const &p1, Polynomial const &p2) {
  size_t max_size = std::max(p1.coefs.size(), p2.coefs.size());
  Polynomial result(max_size);
  
  for (size_t i = 0; i < max_size; ++i)
    result.coefs[i] = p1[i] + p2[i];
  
  return result;
}

Polynomial operator+(Polynomial const &p, float c) {
  return p + Polynomial({c});
}

Polynomial operator-(Polynomial const &p1, Polynomial const &p2) {
  Polynomial result(std::max(p1.coefs.size(), p2.coefs.size()));
  for (size_t i = 0; i < result.coefs.size(); ++i) {
    float a = (i < p1.coefs.size()) ? p1.coefs[i] : 0.0f;
    float b = (i < p2.coefs.size()) ? p2.coefs[i] : 0.0f;
    result.coefs[i] = a - b;
  }
  return result;
}

Polynomial operator-(Polynomial const &p, float c) {
  return p + Polynomial({-c});
}

Polynomial operator*(Polynomial const &p1, Polynomial const &p2) {
  size_t size = p1.coefs.size() + p2.coefs.size() - 1;
  Polynomial result(size);
  
  for (size_t i = 0; i < size; ++i) {
    float s = 0.0f;
    for (size_t j = std::max(0, int(i) - int(p2.coefs.size()) + 1); j <= i && j < p1.coefs.size(); ++j)
      s += p1.coefs[j] * p2.coefs[i - j];
    result.coefs[i] = s;
  }
  return result;
}

float smallest_positive_root(Polynomial p) {
  float root;
  float result = std::numeric_limits<float>::infinity();
  while (p.find_root(root)) {
    if (root > 0.0f && root < result)
      result = root;
    p = p.divide_by_x_minus(root);
  }
  return result;
}

template <typename S>
struct Vector {
  S x, y, z;
  
  Vector(S x, S y, S z) : x(x), y(y), z(z) {
  }
};

template <typename S>
Vector<S> operator*(S scale, Vector<S> const &v) {
  return Vector<S>(scale * v.x, scale * v.y, scale * v.z);
}

template <typename S>
Vector<S> operator*(Vector<S> const &v, S scale) {
  return scale * v;
}

template <typename S>
Vector<S> operator+(Vector<S> const &v1, Vector<S> const &v2) {
  return Vector<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
}

template <typename S>
Vector<S> operator-(Vector<S> const &v1, Vector<S> const &v2) {
  return Vector<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
}

template <typename S>
S dot_product(Vector<S> const &v1, Vector<S> const &v2) {
  return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

template <typename S>
S length_squared(Vector<S> const &v) {
  return dot_product(v, v);
}

Vector<Polynomial> operator*(Vector<float> const &v, Polynomial scalar) {
  return Vector<Polynomial>(v.x * scalar, v.y * scalar, v.z * scalar);
}

Vector<Polynomial> operator*(Polynomial scalar, Vector<float> const &v) {
  return v * scalar;
}

Vector<Polynomial> operator+(Vector<float> const &v1, Vector<Polynomial> const &v2) {
  return Vector<Polynomial>(Polynomial({v1.x}) + v2.x,
			    Polynomial({v1.y}) + v2.y,
			    Polynomial({v1.z}) + v2.z);
}

Vector<Polynomial> operator+(Vector<Polynomial> const &v1, Vector<float> const &v2) {
  return v2 + v1;
}

// --- START READING HERE ---

// This code assumes we are using coordinates where the shooter is at the origin and has velocity 0
float time_to_interception(Vector<float> target_position,
			   Vector<float> target_velocity,
			   Vector<float> gravity,
			   float initial_bullet_speed) {
  Polynomial t({0.0f, 1.0f});
  Polynomial equation = length_squared(0.5f * gravity * t * t - (target_position + target_velocity * t))
    - initial_bullet_speed * initial_bullet_speed * t * t;
  
  return smallest_positive_root(equation);
}

Vector<float> initial_velocity_to_reach_point_with_gravity(Vector<float> point, float t, Vector<float> gravity) {
  return point * (1.0f / t) - 0.5f * gravity * t;
}

std::ostream &operator<<(std::ostream &os, Vector<float> const &v) {
  return os << '(' << v.x << ',' << v.y << ',' << v.z << ')';
}

int main()  { 
  Vector<float> y0(3.0f, 4.0f, 0.0f);
  Vector<float> yv(1.0f, 0.0f, 0.0f);
  Vector<float> g(0.0f, 0.0f, -9.80665f);
  float s0 = 20.0f;
  
  std::cout << "Initial target position: " << y0 << '\n';
  std::cout << "Target velocity: " << yv << '\n';
  std::cout << "Gravity: " << g << '\n';
  std::cout << "Initial bullet speed: " << s0 << '\n';
  std::cout << '\n';
  
  float t = time_to_interception(y0, yv, g, s0);
  Vector<float> interception_position = y0 + yv * t;
  Vector<float> v0 = initial_velocity_to_reach_point_with_gravity(interception_position, t, g);
  
  std::cout << "Time to interception: " << t << '\n';
  std::cout << "Target will be at (" << interception_position.x << ',' << interception_position.y << ',' << interception_position.z << ")\n";
  std::cout << "We should shoot the bullet with the initial velocity (" << v0.x << ',' << v0.y << ',' << v0.z << ")\n";
  
  std::cout << "The length of that initial velocity is " << std::sqrt(length_squared(v0)) << " (sanity check) \n";
}