Probability of winning in a Die-throw game

Given that 2 players are playing a die-throw game. The game is one player throws a die and he got a point, he can move forward accordingly to the point. All occurrences to get the point have equal probability. Let player1 start a point x and player2 start a point y. Player 1 can receive a point up to A, and Player 2 can receive a point up to B.  When a player in the position p and point got in the die q, he can move forward a minimum of (p+q) and W. The first player to reach W wins the game. Find the probability modulo that player 1 wins if he goes first modulo 998244353.

To find a probability modulo of 998244353. We can represent the probability as an irreducible fraction y/x, then x is indivisible by 998244353. So, there is a unique integer z between 0 and 998244352 such that xz ≡ y(mod998244353). Find this z.

Examples:

Input W = 6, x = 3, y =2, A = 5, B = 4
Output 586363698

Input W = 5, x = 3, y = 2, A = 1, B = 1
Output 0

Approach: This can be solved with the following idea:

This problem can be solved using dynamic programming. Let dp_{i, j, f} be the probability that player1 wins when player1 and player2  are at points i and j respectively, where f=0 if the turn is player1’s and f=1 if player2’s. The sought value is dp_{A, B, 0.} We first consider the value dp_{i, j, 0.} For k=1, 2, …A, player1 moves forward to point min(i+k, W) with equal probability, so  dp_{i, j, 0} = 1/A∑_A^k=1dp_{min(i+k, W), j, 1}. Similarly, we consider dp_{i, j, 1}  for k = 1, 2, …B, player2 moves forward to point min(j+k, W) with equal probability, so  dp_{i, j, 0} = 1/A∑_B^k=1dp_{i, min(j+k, W), 0.}​
By the relations of the DP array above, we can evaluate the DP array in descending order of i and j. We could find the answer using dp_{N, i, f} =1 and dp_{i, N, f} = 0 for i<W, we could find the answer.

Define a 3D DP array dpi, j, f where dpi, j, f is the probability that player 1 wins when player 1 is at position i and player 2 is at position j, and f is 0 if it’s player 1’s turn and 1 if it’s player 2’s turn.

Below are the steps involved in the implementation of the code:

• Initialize the DP array as follows:
  • Set dp_{i, j, 0} = 1 if i = N and dp_{i, j, 1} = 0 if i = N.
  • Set dp_{i, j, 0} = 0 if j = N and dp_{i, j, 1} = 1 if j = N.
  • For all other values of i, j, and f, set dp_{i, j, f} = -1 (to indicate that it hasn’t been computed yet).
• Evaluate the DP array in descending order of i and j, and for each f value (0 for player 1’s turn, 1 for player 2’s turn):
  • If dpi, j, and f have already been computed, skip to the next value.
  • If f = 0 (player 1’s turn), set dpi, j, 0 as follows: 
    •  For k = 1 to A, compute the new position i’ = min(i + k, W).
    • Compute the average of dp_{i’, j, 1} overall valid i’ values (i.e., i’ ≤ W).
    • Set dp_{i, j, 0} = that average.
  • If f = 1 (player 2’s turn), set dp_{i, j, 1} as follows:
    • For k = 1 to B, compute the new position j’ = min(j + k, W).
    • Compute the average of dp_{i, j’, 0} overall valid j’ values (i.e., j’ ≤ W).
    • Set dp_{i, j, 1} = that average.
  • The answer to the problem is dp_{i, j, 0}, where i = 0 and j = 0.

Below are the steps involved in the implementation of the code:

C++

// C++ code for the above approach:
#include <iostream>
using namespace std;

const int MOD = 998244353;

int dp[110][110][2];

// Function to find probability
int mdp(int W, int x, int y, int A, int B)
{

    // Initialize base cases
    for (int i = 0; i <= W; i++) {
        for (int f = 0; f < 2; f++) {

            // Base case 1: if W is
            // reached on one axis, set
            // the value to 1
            dp[W][i][f] = 1;

            // Base case 2: if W is reached
            // on the other axis, set the
            // value to 0
            dp[i][W][f] = 0;
        }
    }

    // Fill in the DP table
    for (int i = W - 1; i >= 0; i--) {
        for (int j = W - 1; j >= 0; j--) {
            for (int k = 1; k <= A; k++) {

                // If going in the x-axis,
                // add the value from the
                // previous cell in
                // the y-axis
                dp[i][j][0] += dp[min(W, i + k)][j][1];
            }

            // Take modulo
            dp[i][j][0] %= MOD;

            // Multiply by modular
            // inverse of A
            dp[i][j][0] *= 1LL * (MOD - 1) / A;
            dp[i][j][0] %= MOD;
            for (int k = 1; k <= B; k++) {
                dp[i][j][1] += dp[i][min(W, j + k)][0];
            }

            // Take modulo
            dp[i][j][1] %= MOD;
            dp[i][j][1] *= 1LL * (MOD - 1) / B;
            dp[i][j][1] %= MOD;
        }
    }

    // Output the answer
    return dp[x][y][0];
}

// Drivers Code
int main()
{
    int W = 6, x = 3, y = 2, A = 5, B = 4;

    // Function call
    cout << mdp(W, x, y, A, B);
    return 0;
}

Java

// JAVA code for the above approach:

class GFG {
    static final int MOD = 998244353;
    static int[][][] dp;

    // Function to find probability
    static int mdp(int W, int x, int y, int A, int B)
    {

        // Initialize the DP table
        dp = new int[W + 1][W + 1][2];

        // Initialize base cases
        for (int i = 0; i <= W; i++) {
            for (int f = 0; f < 2; f++) {

                // Base case 1: if W is reached on one axis,
                // set the value to 1
                dp[W][i][f] = 1;

                // Base case 2: if W is reached on the other
                // axis, set the value to 0
                dp[i][W][f] = 0;
            }
        }

        // Fill in the DP table
        for (int i = W - 1; i >= 0; i--) {
            for (int j = W - 1; j >= 0; j--) {
                for (int k = 1; k <= A; k++) {

                    // If going in the x-axis, add the value
                    // from the previous cell in the y-axis
                    dp[i][j][0]
                        += dp[Math.min(W, i + k)][j][1];
                }

                // Take modulo
                dp[i][j][0] %= MOD;

                // Multiply by modular inverse of A
                dp[i][j][0] *= (MOD - 1) / A;
                dp[i][j][0] %= MOD;

                for (int k = 1; k <= B; k++) {
                    dp[i][j][1]
                        += dp[i][Math.min(W, j + k)][0];
                }

                // Take modulo
                dp[i][j][1] %= MOD;
                dp[i][j][1] *= (MOD - 1) / B;
                dp[i][j][1] %= MOD;
            }
        }

        // Output the answer
        return dp[x][y][0];
    }

    // Driver code
    public static void main(String[] args)
    {
        int W = 6, x = 3, y = 2, A = 5, B = 4;

        // Function call
        System.out.println(mdp(W, x, y, A, B));
    }
}

// This code is contributed by shivamgupta310570

C#

// C# code for the above approach:

using System;
using System.Collections.Generic;

class Program {
    static int MOD = 998244353;
    static int[, , ] dp;

    // Function to find probability
    static int mdp(int W, int x, int y, int A, int B)
    {

        // Initialize the DP table
        dp = new int[W + 1, W + 1, 2];

        // Initialize base cases
        for (int i = 0; i <= W; i++) {
            for (int f = 0; f < 2; f++) {

                // Base case 1: if W is reached on one axis,
                // set the value to 1
                dp[W, i, f] = 1;

                // Base case 2: if W is reached on the other
                // axis, set the value to 0
                dp[i, W, f] = 0;
            }
        }

        // Fill in the DP table
        for (int i = W - 1; i >= 0; i--) {
            for (int j = W - 1; j >= 0; j--) {
                for (int k = 1; k <= A; k++) {

                    // If going in the x-axis, add the value
                    // from the previous cell in the y-axis
                    dp[i, j, 0]
                        += dp[Math.Min(W, i + k), j, 1];
                }

                // Take modulo
                dp[i, j, 0] %= MOD;

                // Multiply by modular inverse of A
                dp[i, j, 0] *= (MOD - 1) / A;
                dp[i, j, 0] %= MOD;

                for (int k = 1; k <= B; k++) {
                    dp[i, j, 1]
                        += dp[i, Math.Min(W, j + k), 0];
                }

                // Take modulo
                dp[i, j, 1] %= MOD;
                dp[i, j, 1] *= (MOD - 1) / B;
                dp[i, j, 1] %= MOD;
            }
        }

        // Output the answer
        return dp[x, y, 0];
    }

    // Driver code
    static void Main(string[] args)
    {
        int W = 6, x = 3, y = 2, A = 5, B = 4;

        // Function call
        Console.WriteLine(mdp(W, x, y, A, B));
    }
}

// This code is contributed by Tapesh(tapeshdua420)

Output

586363698

Time Complexity: O(W² * A * B)
Auxiliary Space: O(W²)