1. Concise problem statement
Given positive integers x1, x2, a target position P (can be negative), and a total number of jumps K, decide whether it is possible to make exactly K jumps—each jump is either +x1, –x1, +x2, or –x2—so that the net displacement is exactly P. If yes, output “YES” and one quadruple (P1, N1, P2, N2) of nonnegative integers satisfying
     P1 + N1 + P2 + N2 = K
     P1·x1 − N1·x1 + P2·x2 − N2·x2 = P
otherwise output “NO.”

2. Detailed editorial

Let a = P1 − N1 and b = P2 − N2. We then need

  (1)  a·x1 + b·x2 = P
  (2)  |a| + |b| ≤ K
  (3)  |a| + |b| ≡ K mod 2

Indeed, once we fix integers a,b satisfying (1), we can set P1,N1 so that P1−N1=a and P1+N1 contributes exactly |a| jumps, and similarly for x2. The remaining K−(|a|+|b|) jumps must be inserted in ± pairs that cancel out displacement; that is possible iff K−(|a|+|b|) is even.

Step A: Solve a·x1 + b·x2 = P
  Use the extended Euclidean algorithm to find one particular integer solution (a0,b0). Let g = gcd(x1,x2). If g does not divide P, there is no solution: answer “NO.”

Step B: Parametrize all solutions
  General solution is
    a = a0 + t·(x2/g),
    b = b0 − t·(x1/g),
  for integer t.

Step C: Choose t to minimize S(t) = |a| + |b|
  Try to shift t so that a and b come closer to zero. Since S(t) is piecewise-linear, you can approximate the “best” real t and check its neighbors (or simply do two while-loops pushing t ±1 as long as it decreases S).

Step D: Check feasibility
  Let S0 = |a| + |b|. If S0 > K, impossible: answer “NO.”
  Otherwise compute rem = K − S0.
  If rem is even, proceed.
  If rem is odd, we need one extra “unit step” in the solution family to flip parity of S: because moving t by +1 changes S by |a+dx|+|b−dy|−(|a|+|b|) where dx = x2/g, dy = x1/g. We check if (dx+dy) is odd (so parity of S flips). If not odd, no solution; else shift t by +1 or −1 in the direction that keeps S ≤ K. Recompute S0 and rem; now rem will be even.

Step E: Build P1,N1,P2,N2
  For x1: if a ≥ 0 then P1 = a, N1 = 0; else P1 = 0, N1 = −a. Similarly for b and x2.
  Then distribute rem/2 extra +x1 and rem/2 extra −x1 jumps to keep net a unchanged; i.e. add rem/2 to both P1 and N1. That uses up all K jumps.

Time complexity: O(log (min(x1,x2))) for gcd plus a few corrections.

3. C++ Solution
```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

int64_t x1v, x2v, p, k;

void read() { cin >> x1v >> x2v >> p >> k; }

int64_t extended_gcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
    if(b == 0) {
        x = 1;
        y = 0;
        return a;
    }

    int64_t x1, y1;
    int64_t g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

void solve() {
    // We need P1*x1 - N1*x1 + P2*x2 - N2*x2 = P with P1+N1+P2+N2 = K and all
    // four counts non-negative. Let a = P1 - N1 and b = P2 - N2; then
    // a*x1 + b*x2 = P, which is solvable via the extended Euclid only when
    // gcd(x1,x2) divides P. Among all solutions (shifting a by x2/g and b by
    // x1/g) we slide to the one minimising |a| + |b|, the smallest number of
    // jumps the net displacement forces.
    //
    // The leftover jumps last = K - (|a| + |b|) must be spent in
    // cancelling pairs (one positive plus one negative of the same length),
    // so last has to be non-negative and even. If it is odd we try shifting
    // by one (a +/- x2/g, b -/+ x1/g): this changes the parity of |a| + |b|
    // only when x1/g + x2/g is odd, otherwise no solution exists. Finally we
    // split each net count into positive and negative jumps and add half of
    // the leftover to a matching positive/negative pair.

    int64_t p1, p2;
    int64_t g = extended_gcd(x1v, x2v, p1, p2);

    if(p % g != 0) {
        cout << "NO\n";
        return;
    }

    int64_t dx = x2v / g;
    int64_t dy = x1v / g;
    p1 *= p / g;
    p2 *= p / g;

    while(llabs(p1 + dx) + llabs(p2 - dy) < llabs(p1) + llabs(p2)) {
        p1 += dx;
        p2 -= dy;
    }

    while(llabs(p1 - dx) + llabs(p2 + dy) < llabs(p1) + llabs(p2)) {
        p1 -= dx;
        p2 += dy;
    }

    if(llabs(p1) + llabs(p2) > k) {
        cout << "NO\n";
        return;
    }

    int64_t n1 = 0, n2 = 0;
    int64_t last = k - llabs(p1) - llabs(p2);
    if(last % 2 == 0) {
        if(p1 < 0) {
            n1 = -p1;
            p1 = 0;
        }

        if(p2 < 0) {
            n2 = -p2;
            p2 = 0;
        }

        p1 += last / 2;
        n1 += last / 2;
    } else {
        if((dx + dy) % 2 == 0) {
            cout << "NO\n";
            return;
        }

        if(llabs(p1 + dx) + llabs(p2 - dy) < llabs(p1 - dx) + llabs(p2 + dy)) {
            p1 += dx;
            p2 -= dy;
        } else {
            p1 -= dx;
            p2 += dy;
        }

        if(llabs(p1) + llabs(p2) > k) {
            cout << "NO\n";
            return;
        }

        last = k - llabs(p1) - llabs(p2);
        if(p1 < 0) {
            n1 = -p1;
            p1 = 0;
        }

        if(p2 < 0) {
            n2 = -p2;
            p2 = 0;
        }

        p1 += last / 2;
        n1 += last / 2;
    }

    cout << "YES\n";
    cout << p1 << ' ' << n1 << ' ' << p2 << ' ' << n2 << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

4. Python solution with detailed comments
```python
def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    """
    Return (g, x, y) such that g = gcd(a, b) and a*x + b*y = g.
    """
    if b == 0:
        return (a, 1, 0)      # Base: gcd(a,0)=a => a*1 + 0*0 = a
    # Recurse on (b, a mod b)
    g, x1, y1 = extended_gcd(b, a % b)
    # Now b*x1 + (a%b)*y1 = g
    # a%b = a - (a//b)*b
    # => a*y1 + b*(x1 - (a//b)*y1) = g
    x = y1
    y = x1 - (a // b) * y1
    return (g, x, y)

def solve_problem(x1: int, x2: int, P: int, K: int) -> None:
    # Step A: find any solution (a0,b0) to a0*x1 + b0*x2 = P
    g, a0, b0 = extended_gcd(x1, x2)
    if P % g != 0:
        print("NO")
        return
    # Scale to get a valid particular solution
    a0 *= P // g
    b0 *= P // g

    dx = x2 // g  # shift in 'a' when changing parameter t by +1
    dy = x1 // g  # shift in 'b' when changing parameter t by +1

    # S(a,b) = |a| + |b|
    def S(a, b): return abs(a) + abs(b)

    # Step C: adjust t to minimize |a|+|b|
    a, b = a0, b0
    # push t upward while it strictly improves S
    while True:
        a_up, b_up = a + dx, b - dy
        if S(a_up, b_up) < S(a, b):
            a, b = a_up, b_up
        else:
            break
    # push t downward similarly
    while True:
        a_dn, b_dn = a - dx, b + dy
        if S(a_dn, b_dn) < S(a, b):
            a, b = a_dn, b_dn
        else:
            break

    sum_abs = S(a, b)
    if sum_abs > K:
        print("NO")
        return
    rem = K - sum_abs

    # Step D: fix parity if needed
    if rem % 2 != 0:
        # We need dx+dy odd to flip parity of S
        if (dx + dy) % 2 == 0:
            print("NO")
            return
        # Try shifting t by +1; if that overshoots, shift by -1
        if S(a+dx, b-dy) <= K:
            a, b = a+dx, b-dy
        else:
            a, b = a-dx, b+dy
        sum_abs = S(a, b)
        rem = K - sum_abs
        if rem < 0 or rem % 2 != 0:
            print("NO")
            return

    # Step E: build P1,N1,P2,N2
    # a = P1 - N1  =>  P1 = max(a,0), N1 = max(-a,0)
    if a >= 0:
        P1, N1 = a, 0
    else:
        P1, N1 = 0, -a
    if b >= 0:
        P2, N2 = b, 0
    else:
        P2, N2 = 0, -b

    # Distribute remaining rem jumps in cancelling +x1/-x1 pairs:
    add = rem // 2
    P1 += add
    N1 += add

    print("YES")
    print(P1, N1, P2, N2)

def main():
    x1, x2, P, K = map(int, input().split())
    solve_problem(x1, x2, P, K)

if __name__ == "__main__":
    main()
```

5. Compressed editorial
Use extended gcd to find one integer solution (a0,b0) to a·x1 + b·x2 = P. General solutions are a=a0+t·(x2/g), b=b0−t·(x1/g). Adjust t greedily to minimize |a|+|b|. If this minimum exceeds K, or if leftover jumps K–(|a|+|b|) has wrong parity (and cannot be flipped because (x2/g + x1/g) is even), answer NO. Otherwise distribute leftover as cancelling ± jumps to fill exactly K moves and print the resulting nonnegative counts.
