Jump to content

Dice game probabilities

From Rosetta Code
Task
Dice game probabilities
You are encouraged to solve this task according to the task description, using any language you may know.

Two players have a set of dice each. The first player has nine dice with four faces each, with numbers one to four. The second player has six normal dice with six faces each, each face has the usual numbers from one to six.

They roll their dice and sum the totals of the faces. The player with the highest total wins (it's a draw if the totals are the same). What's the probability of the first player beating the second player?

Later the two players use a different set of dice each. Now the first player has five dice with ten faces each, and the second player has six dice with seven faces each. Now what's the probability of the first player beating the second player?

This task was adapted from the Project Euler Problem n.205: https://projecteuler.net/problem=205

11l

Translation of: C
F throw_die(n_sides, n_dice, s, [Int] &counts)
   I n_dice == 0
      counts[s]++
      R
   L(i) 1..n_sides
      throw_die(n_sides, n_dice - 1, s + i, counts)

F beating_probability(n_sides1, n_dice1,
                      n_sides2, n_dice2)
   V len1 = (n_sides1 + 1) * n_dice1
   V C1 = [0] * len1
   throw_die(n_sides1, n_dice1, 0, &C1)

   V len2 = (n_sides2 + 1) * n_dice2
   V C2 = [0] * len2
   throw_die(n_sides2, n_dice2, 0, &C2)

   Float p12 = (n_sides1 ^ n_dice1) * (n_sides2 ^ n_dice2)

   V tot = 0.0
   L(i) 0 .< len1
      L(j) 0 .< min(i, len2)
         tot += Float(C1[i]) * C2[j] / p12
   R tot

print(‘#.16’.format(beating_probability(4, 9, 6, 6)))
print(‘#.16’.format(beating_probability(10, 5, 7, 6)))
Output:
0.5731440767829814
0.6427886287176272


Action!

Library: Action! Tool Kit
INCLUDE "D2:REAL.ACT" ;from the Action! Tool Kit

BYTE FUNC Roll(BYTE sides,dices)
  BYTE i,sum

  sum=0
  FOR i=1 TO dices
  DO
    sum==+Rand(sides)+1
  OD
RETURN (sum)

PROC Test(BYTE sides1,dices1,sides2,dices2)
  INT i,count=[10000],sum1,sum2,wins1,wins2,draws
  REAL r1,r2,p

  wins1=0 wins2=0 draws=0
  FOR i=1 TO count
  DO
    sum1=Roll(sides1,dices1)
    sum2=Roll(sides2,dices2)
    IF sum1>sum2 THEN
      wins1==+1
    ELSEIF sum1<sum2 THEN
      wins2==+1
    ELSE
      draws==+1
    FI
  OD

  IntToReal(wins1,r1)
  IntToReal(wins2,r2)
  RealDiv(r2,r1,p)

  PrintF("Tested %I times%E",count)
  PrintF("Player 1 with %B dice of %B sides%E",dices1,sides1)
  PrintF("Player 2 with %B dice of %B sides%E",dices2,sides2)
  PrintF("Player 1 wins %I times%E",wins1)
  PrintF("Player 2 wins %I times%E",wins2)
  PrintF("Draw %I times%E",draws)
  Print("Player 1 beating Player 2 probability:")
  PrintRE(p)
  PutE()
RETURN

PROC Main()
  Put(125) PutE() ;clear the screen

  Test(4,9,6,6)
  Test(10,5,7,6)
RETURN
Output:

Screenshot from Atari 8-bit computer 

Tested 10000 times
Player 1 with 9 dice of 4 sides
Player 2 with 6 dice of 6 sides
Player 1 wins 5770 times
Player 2 wins 3493 times
Draw 737 times
Player 1 beating Player 2 probability:
.6053726169

Tested 10000 times
Player 1 with 5 dice of 10 sides
Player 2 with 6 dice of 7 sides
Player 1 wins 6417 times
Player 2 wins 3140 times
Draw 443 times
Player 1 beating Player 2 probability:
.4893252298

Ada 

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Discrete_Random;

procedure Main is
   package real_io is new Float_IO (Long_Float);
   use real_io;

   type Dice is record
      Faces    : Positive;
      Num_Dice : Positive;
   end record;

   procedure Roll_Dice (The_Dice : in Dice; Count : out Natural) is
      subtype Faces is Integer range 1 .. The_Dice.Faces;
      package Die_Random is new Ada.Numerics.Discrete_Random (Faces);
      use Die_Random;
      Seed : Generator;
   begin
      Reset (Seed);
      Count := 0;
      for I in 1 .. The_Dice.Num_Dice loop
         Count := Count + Random (Seed);
      end loop;
   end Roll_Dice;

   function Win_Prob
     (Dice_1 : Dice; Dice_2 : Dice; Tries : Positive) return Long_Float
   is
      Count_1      : Natural := 0;
      Count_2      : Natural := 0;
      Count_1_Wins : Natural := 0;
   begin
      for I in 1 .. Tries loop
         Roll_Dice (Dice_1, Count_1);
         Roll_Dice (Dice_2, Count_2);
         if Count_1 > Count_2 then
            Count_1_Wins := Count_1_Wins + 1;
         end if;
      end loop;
      return Long_Float (Count_1_Wins) / Long_Float (Tries);
   end Win_Prob;

   D1 : Dice := (Faces => 4, Num_Dice => 9);
   D2 : Dice := (Faces => 6, Num_Dice => 6);
   D3 : Dice := (Faces => 10, Num_Dice => 5);
   D4 : Dice := (Faces => 7, Num_Dice => 6);

   P1 : Long_Float := Win_Prob (D1, D2, 1_000_000);
   P2 : Long_Float := Win_Prob (D3, D4, 1_000_000);
begin
   Put ("Dice D1 wins = ");
   Put (Item => P1, Fore => 1, Aft => 7, Exp => 0);
   New_Line;
   Put ("Dice D2 wins = ");
   Put (Item => P2, Fore => 1, Aft => 7, Exp => 0);
   New_Line;
end Main;
Output:
Dice D1 wins = 0.5727800
Dice D2 wins = 0.6427660

BASIC

BASIC256

Translation of: Yabasic
dado1 = 9: lado1 = 4
dado2 = 6: lado2 = 6
total1 = 0: total2 = 0

for i = 0 to 1
    for cont = 1 to 100000
        jugador1 = lanzamiento(dado1, lado1)
        jugador2 = lanzamiento(dado2, lado2)
        if jugador1 > jugador2 then
            total1 = total1 + 1
        else
            if jugador1 <> jugador2 then total2 = total2 + 1
        endif
    next cont

    print "Lanzado el dado "; (cont - 1); " veces"
    print "jugador1 con "; dado1; " dados de "; lado1; " lados"
    print "jugador2 con "; dado2; " dados de "; lado2; " lados"
    print "Total victorias jugador1 = "; total1; " => "; left(string(total2 / total1), 9)
    print "Total victorias jugador2 = "; total2
    print (cont - 1) - (total1 + total2); " empates" + chr(10)

    dado1 = 5: lado1 = 10
    dado2 = 6: lado2 = 7
    total1 = 0: total2 = 0
next i
end

function  lanzamiento(dado, lado)
    total = 0

    for lanza = 1 to dado
        total = total + int(rand * lado) + 1
    next lanza
    return total
end function
Output:
Igual que la entrada de Yabasic.

FreeBASIC

Translation of: Gambas
Dim As Integer lado, jugador1, jugador2, total1, total2, cont, i
Dim As Integer dado1 = 9, lado1 = 4
Dim As Integer dado2 = 6, lado2 = 6

Randomize Timer

Function Lanzamiento(dado As Integer, lado As Integer) As Integer
    Dim As Short lanza, total
    
    For lanza = 1 To dado
        total += Int(Rnd * lado) + 1
    Next lanza
    Return total
End Function 

For i = 0 To 1
    For cont = 1 To 100000
        jugador1 = Lanzamiento(dado1, lado1)
        jugador2 = Lanzamiento(dado2, lado2)
        If jugador1 > jugador2 Then
            total1 += 1
        Elseif jugador1 <> jugador2 Then
            total2 += 1
        End If
    Next cont
    
    Print Using "Lanzado el dado & veces"; (cont - 1)
    Print "jugador1 con"; dado1; " dados de"; lado1; " lados"
    Print "jugador2 con"; dado2; " dados de"; lado2; " lados"
    Print Using "Total victorias jugador1 = &  => #.#######"; total1; (total2 / total1)
    Print "Total victorias jugador2 ="; total2
    Print (cont - 1) - (total1 + total2); !" empates\n"
    
    dado1 = 5: lado1 = 10
    dado2 = 6: lado2 = 7
    total1 = 0: total2 = 0
Next i

Sleep
Output:
Lanzado el dado 100000 veces
jugador1 con 9 dados de 4 lados
jugador2 con 6 dados de 6 lados
Total victorias jugador1 = 57274  => 0.6237211
Total victorias jugador2 = 35723
 7003 empates

Lanzado el dado 100000 veces
jugador1 con 5 dados de 10 lados
jugador2 con 6 dados de 7 lados
Total victorias jugador1 = 64093  => 0.4893826
Total victorias jugador2 = 31366
 4541 empates

Yabasic

Translation of: FreeBASIC
dado1 = 9: lado1 = 4
dado2 = 6: lado2 = 6
total1 = 0: total2 = 0

for i = 0 to 1
    for cont = 1 to 100000
        jugador1 = lanzamiento(dado1, lado1)
        jugador2 = lanzamiento(dado2, lado2)
        if jugador1 > jugador2 then
            total1 = total1 + 1
        elseif jugador1 <> jugador2 then
            total2 = total2 + 1
        endif
    next cont
    
    print "Lanzado el dado ", (cont - 1), " veces" 
    print "jugador1 con ", dado1, " dados de ", lado1, " lados"
    print "jugador2 con ", dado2, " dados de ", lado2, " lados"
    print "Total victorias jugador1 = ", total1, " => ", (total2 / total1)
    print "Total victorias jugador2 = ", total2
    print (cont - 1) - (total1 + total2), " empates\n"
    
    dado1 = 5: lado1 = 10
    dado2 = 6: lado2 = 7
    total1 = 0: total2 = 0
next i

sub lanzamiento(dado, lado)
    local lanza, total
    
    for lanza = 1 to dado
        total = total + int(ran(lado)) + 1
    next lanza
    return total
end sub
Output:
Igual que la entrada de FreeBASIC.

C 

#include <stdio.h>
#include <stdint.h>

typedef uint32_t uint;
typedef uint64_t ulong;

ulong ipow(const uint x, const uint y) {
    ulong result = 1;
    for (uint i = 1; i <= y; i++)
        result *= x;
    return result;
}

uint min(const uint x, const uint y) {
    return (x < y) ? x : y;
}

void throw_die(const uint n_sides, const uint n_dice, const uint s, uint counts[]) {
    if (n_dice == 0) {
        counts[s]++;
        return;
    }

    for (uint i = 1; i < n_sides + 1; i++)
        throw_die(n_sides, n_dice - 1, s + i, counts);
}

double beating_probability(const uint n_sides1, const uint n_dice1,
                           const uint n_sides2, const uint n_dice2) {
    const uint len1 = (n_sides1 + 1) * n_dice1;
    uint C1[len1];
    for (uint i = 0; i < len1; i++)
        C1[i] = 0;
    throw_die(n_sides1, n_dice1, 0, C1);

    const uint len2 = (n_sides2 + 1) * n_dice2;
    uint C2[len2];
    for (uint j = 0; j < len2; j++)
        C2[j] = 0;
    throw_die(n_sides2, n_dice2, 0, C2);

    const double p12 = (double)(ipow(n_sides1, n_dice1) * ipow(n_sides2, n_dice2));

    double tot = 0;
    for (uint i = 0; i < len1; i++)
        for (uint j = 0; j < min(i, len2); j++)
            tot += (double)C1[i] * C2[j] / p12;
    return tot;
}

int main() {
    printf("%1.16f\n", beating_probability(4, 9, 6, 6));
    printf("%1.16f\n", beating_probability(10, 5, 7, 6));
    return 0;
}
Output:
0.5731440767829801
0.6427886287176260

C#

Translation of: Java
using System;

class Program
{
    static uint MinOf(uint x, uint y)
    {
        return x < y ? x : y;
    }

    static void ThrowDie(uint nSides, uint nDice, uint s, uint[] counts)
    {
        if (nDice == 0)
        {
            counts[s]++;
            return;
        }
        for (uint i = 1; i <= nSides; i++)
        {
            ThrowDie(nSides, nDice - 1, s + i, counts);
        }
    }

    static double BeatingProbability(uint nSides1, uint nDice1, uint nSides2, uint nDice2)
    {
        uint len1 = (nSides1 + 1) * nDice1;
        uint[] c1 = new uint[len1]; // initialized to zero by default
        ThrowDie(nSides1, nDice1, 0, c1);

        uint len2 = (nSides2 + 1) * nDice2;
        uint[] c2 = new uint[len2];
        ThrowDie(nSides2, nDice2, 0, c2);
        double p12 = Math.Pow(nSides1, nDice1) * Math.Pow(nSides2, nDice2);

        double tot = 0.0;
        for (uint i = 0; i < len1; i++)
        {
            for (uint j = 0; j < MinOf(i, len2); j++)
            {
                tot += (double)(c1[i] * c2[j]) / p12;
            }
        }
        return tot;
    }

    static void Main(string[] args)
    {
        Console.WriteLine(BeatingProbability(4, 9, 6, 6));
        Console.WriteLine(BeatingProbability(10, 5, 7, 6));
    }
}
Output:
0.573144076782982
0.642788628717627

C++ 

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <map>

// Returns map from sum of faces to number of ways it can happen
std::map<uint32_t, uint32_t> get_totals(uint32_t dice, uint32_t faces) {
    std::map<uint32_t, uint32_t> result;
    for (uint32_t i = 1; i <= faces; ++i)
        result.emplace(i, 1);
    for (uint32_t d = 2; d <= dice; ++d) {
        std::map<uint32_t, uint32_t> tmp;
        for (const auto& p : result) {
            for (uint32_t i = 1; i <= faces; ++i)
                tmp[p.first + i] += p.second;
        }
        tmp.swap(result);
    }
    return result;
}

double probability(uint32_t dice1, uint32_t faces1, uint32_t dice2, uint32_t faces2) {
    auto totals1 = get_totals(dice1, faces1);
    auto totals2 = get_totals(dice2, faces2);
    double wins = 0;
    for (const auto& p1 : totals1) {
        for (const auto& p2 : totals2) {
            if (p2.first >= p1.first)
                break;
            wins += p1.second * p2.second;
        }
    }
    double total = std::pow(faces1, dice1) * std::pow(faces2, dice2);
    return wins/total;
}

int main() {
    std::cout << std::setprecision(10);
    std::cout << probability(9, 4, 6, 6) << '\n';
    std::cout << probability(5, 10, 6, 7) << '\n';
    return 0;
}
Output:
0.5731440768
0.6427886287

Common Lisp

DRAFT

1. Note

We use total probability rule.

2. Program

Player i gets si points as the sum of rolling ni dices, each with fi faces. 

(defun n-ways (n f s &optional (mem (make-array '(999 99) :initial-element -1)))
  (cond ((and (= n 0) (= s 0)) 1)
        ((or (= n 0) (<= s 0)) 0)
        ((>= (aref mem s n) 0) (aref mem s n))
        (t (loop for i from 1 to f
                 sum (n-ways (1- n) f (- s i) mem) into total
                 finally (return (setf (aref mem s n) total))))))

(defun winning-probability (n1 f1 n2 f2 &aux (w 0))
  (loop for i from n1 to (* n1 f1)
        do (loop for j from n2 to (* n2 f2)
                 do (if (> i j)
                        (setf w (+ w (* (n-ways n1 f1 i) (n-ways n2 f2 j))))))
        finally (return (/ w (* (expt f1 n1) (expt f2 n2))))))

3. Execution

(winning-probability 9 4 6 6)
(winning-probability 5 10 6 7)
Output:
48679795/84934656      ; ≈ 0.573
3781171969/5882450000  ; ≈ 0.643

That's all Folks !

cyril nocton (cyril.nocton@gmail.com) w/ google translate

D

version 1

import std.stdio, std.range, std.algorithm;

void throwDie(in uint nSides, in uint nDice, in uint s, uint[] counts)
pure nothrow @safe @nogc {
    if (nDice == 0) {
        counts[s]++;
        return;
    }

    foreach (immutable i; 1 .. nSides + 1)
        throwDie(nSides, nDice - 1, s + i, counts);
}

real beatingProbability(uint nSides1, uint nDice1,
                        uint nSides2, uint nDice2)()
pure nothrow @safe /*@nogc*/ {
    uint[(nSides1 + 1) * nDice1] C1;
    throwDie(nSides1, nDice1, 0, C1);

    uint[(nSides2 + 1) * nDice2] C2;
    throwDie(nSides2, nDice2, 0, C2);

    immutable p12 = real((ulong(nSides1) ^^ nDice1) *
                         (ulong(nSides2) ^^ nDice2));

    return cartesianProduct(C1[].enumerate, C2[].enumerate)
           .filter!(p => p[0][0] > p[1][0])
           .map!(p => real(p[0][1]) * p[1][1] / p12)
           .sum;
}

void main() @safe {
    writefln("%1.16f", beatingProbability!(4, 9, 6, 6));
    writefln("%1.16f", beatingProbability!(10, 5, 7, 6));
}
Output:
0.5731440767829801
0.6427886287176262

version 2 (Faster Alternative Version)

Translation of: Python
import std.stdio, std.range, std.algorithm;

ulong[] combos(R)(R sides, in uint n) pure nothrow @safe
if (isForwardRange!R) {
    if (sides.empty)
        return null;
    if (!n)
        return [1];
    auto ret = new typeof(return)(reduce!max(sides[0], sides[1 .. $]) * n + 1);
    foreach (immutable i, immutable v; enumerate(combos(sides, n - 1))) {
        if (!v)
            continue;
        foreach (immutable s; sides)
            ret[i + s] += v;
    }
    return ret;
}

real winning(R)(R sides1, in uint n1, R sides2, in uint n2)
pure nothrow @safe if (isForwardRange!R) {
    static void accumulate(T)(T[] arr) pure nothrow @safe @nogc {
        foreach (immutable i; 1 .. arr.length)
            arr[i] += arr[i - 1];
    }

    immutable p1 = combos(sides1, n1);
    auto p2 = combos(sides2, n2);
    immutable s = p1.sum * p2.sum;
    accumulate(p2);
    ulong win = 0; // 'win' is 1 beating 2.
    foreach (immutable i, immutable x1; p1.dropOne.enumerate)
        win += x1 * p2[min(i, $ - 1)];
    return win / real(s);
}

void main() @safe {
    writefln("%1.16f", winning(iota(1u, 5u),  9, iota(1u, 7u), 6));
    writefln("%1.16f", winning(iota(1u, 11u), 5, iota(1u, 8u), 6));
}
Output:
0.5731440767829801
0.6427886287176262

EasyLang

Translation of: Phix
proc throw_die n_sides n_dice s . counts[] .
   if n_dice = 0
      counts[s] += 1
      return
   .
   for i = 1 to n_sides
      throw_die n_sides (n_dice - 1) (s + i) counts[]
   .
.
func proba n_sides1 n_dice1 n_sides2 n_dice2 .
   len c1[] (n_sides1 + 1) * n_dice1
   len c2[] (n_sides2 + 1) * n_dice2
   throw_die n_sides1 n_dice1 0 c1[]
   throw_die n_sides2 n_dice2 0 c2[]
   p12 = pow n_sides1 n_dice1 * pow n_sides2 n_dice2
   for i = 1 to len c1[]
      for j = 1 to lower (i - 1) len c2[]
         tot += c1[i] * c2[j] / p12
      .
   .
   return tot
.
numfmt 5 0
print proba 4 9 6 6
print proba 10 5 7 6

Factor 

USING: dice generalizations kernel math prettyprint sequences ;
IN: rosetta-code.dice-probabilities

: winning-prob ( a b c d -- p )
    [ [ random-roll ] 2bi@ > ] 4 ncurry [ 100000 ] dip replicate
    [ [ t = ] count ] [ length ] bi /f ;

9 4 6 6 winning-prob
5 10 6 7 winning-prob [ . ] bi@
Output:
0.57199
0.64174

Forth

Works with: gforth version 0.7.3


#! /usr/bin/gforth

\ Dice game probabilities

: min? ( addr -- min )
    @
;

: max? ( addr -- max )
    cell + @
;

: max+1-min? ( addr -- max+1 min )
    dup max? 1+ swap min?
;

: addr? ( addr x -- addr' )
    over min? - 2 + cells +
;

: weight? ( addr x -- w )
    2dup swap min? < IF
        2drop 0
    ELSE
        2dup swap max? > IF
            2drop 0
        ELSE
            addr? @
        THEN
    THEN
;

: total-weight?   ( addr -- w )
    dup max? 1+   ( addr max+1 )
    over min?     ( addr max+1 min )
    0 -rot ?DO ( adrr 0 max+1 min )
        over i weight? +
    LOOP
    nip
;

: uniform-aux ( min max x -- addr )
    >r 2dup
    2dup swap - 3 + cells allocate throw ( min max min max addr )
    tuck cell + !                        ( min max min addr )
    tuck !                               ( min max addr )
    -rot swap                            ( addr max min )
    r> -rot                              ( addr x max min )
    - 3 + 2 ?DO                          ( addr x )
        2dup swap i cells + !
    LOOP
    drop
;

: convolve { addr1 addr2 -- addr }
    addr1 min? addr2 min? +
    addr1 max? addr2 max? +
    0 uniform-aux                    { addr }
    addr1 max+1-min? ?DO
        addr2 max+1-min? ?DO
            addr1 j weight?
            addr2 i weight? *
            addr i j + addr? +!
        LOOP
    LOOP
    addr
;

: even? ( n -- f )
    2 mod 0=
;

: power ( addr exp -- addr' )
    dup 1 = IF
        drop
    ELSE
        dup even? IF
            2/ recurse dup convolve
        ELSE
            over swap 2/ recurse dup convolve convolve
        THEN
    THEN
;

: .dist { addr -- }
    addr total-weight? { tw }
    addr max+1-min? ?DO
        i 10 .r
        addr i weight? dup 20 .r
        0 d>f tw 0 d>f f/ ."  " f. cr
    LOOP
;

: dist-cmp { addr1 addr2 xt -- p }
    0
    addr1 max+1-min? ?DO
        addr2 max+1-min? ?DO
            j i xt execute IF
                addr1 j weight?
                addr2 i weight?
                * +
            THEN
        LOOP
    LOOP
    0 d>f
    addr1 total-weight? addr2 total-weight? um* d>f
    f/
;

: dist> ( addr1 addr2 -- p )
    ['] > dist-cmp
;

\ creates the uniform distribution with outcomes from min to max
: uniform ( min max -- addr )
    1 uniform-aux
;

\ example

1 4 uniform 9 power
1 6 uniform 6 power
dist> f. cr

1 10 uniform 5 power
1  7 uniform 6 power
dist> f. cr

bye
Output:
time ./dice-game-probabilities.fs 
0.57314407678298 
0.642788628717626 

real	0m0,008s
user	0m0,005s
sys	    0m0,003s


FutureBasic 

local fn RollTheDice( dice as long, sides as long ) as long
  long i, total = 0
  
  for i = 1 to dice
    total += rnd(sides)
  next
end fn = total


void local fn CaclulateProbabilities( dice1 as long, sides1 as long, dice2 as long, sides2 as long )
  long  i, player1, player2
  float total1 = 0, total2 = 0, total3 = 0
  
  for i = 1 to 100000
    player1 = fn RollTheDice( dice1, sides1 )
    player2 = fn RollTheDice( dice2, sides2 )
    if ( player1  > player2 ) then total1++ : continue
    if ( player1  < player2 ) then total2++ : continue
    if ( player1 == player2 ) then total3++ : continue
  next
  
  printf @"Dice cast %ld times.", i - 1
  printf @"Player 1 with %ld dice of %ld sides.", dice1, sides1
  printf @"Player 2 with %ld dice of %ld sides.", dice2, sides2
  printf @"Total wins Player 1 = %.f. Probability of winning: %0.4f%%.", total1, (total1 / i ) * 100
  printf @"Total wins Player 2 = %.f", total2
  printf @"Number of ties: %.f.", total3
end fn

randomize
print
fn CaclulateProbabilities( 9, 4, 6, 6 )
print
fn CaclulateProbabilities( 5, 10, 6, 7 )
print

HandleEvents
Output:
Dice cast 100000 times.
Player 1 with 9 dice of 4 sides.
Player 2 with 6 dice of 6 sides.
Total wins Player 1 = 57301. Probability of winning: 57.3004%.
Total wins Player 2 = 35733
Number of ties: 6966.

Dice cast 100000 times.
Player 1 with 5 dice of 10 sides.
Player 2 with 6 dice of 7 sides.
Total wins Player 1 = 64289. Probability of winning: 64.2884%.
Total wins Player 2 = 31296
Number of ties: 4415.


Gambas

Click this link to run this code 

' Gambas module file

Public Sub Main()
Dim iSides, iPlayer1, iPlayer2, iTotal1, iTotal2, iCount, iCount0 As Integer
Dim iDice1 As Integer = 9
Dim iDice2 As Integer = 6
Dim iSides1 As Integer = 4
Dim iSides2 As Integer = 6

Randomize

For iCount0 = 0 To 1
  For iCount = 1 To 100000
    iPlayer1 = Roll(iDice1, iSides1)
    iPlayer2 = Roll(iDice2, iSides2)
    If iPlayer1 > iPlayer2 Then
      iTotal1 += 1
    Else If iPlayer1 <> iPlayer2 Then
      iTotal2 += 1
    Endif
  Next

  Print "Tested " & Str(iCount - 1) & " times"
  Print "Player1 with " & Str(iDice1) & " dice of " & Str(iSides1) & " sides"
  Print "Player2 with " & Str(iDice2) & " dice of " & Str(iSides2) & " sides"
  Print "Total wins Player1 = " & Str(iTotal1) & " - " & Str(iTotal2 / iTotal1)
  Print "Total wins Player2 = " & Str(iTotal2) 
  Print Str((iCount - 1) - (iTotal1 + iTotal2)) & " draws" & gb.NewLine

  iDice1 = 5
  iDice2 = 6
  iSides1 = 10
  iSides2 = 7
  iTotal1 = 0
  iTotal2 = 0
Next

End

Public Sub Roll(iDice As Integer, iSides As Integer) As Integer
Dim iCount, iTotal As Short

For iCount = 1 To iDice
  iTotal += Rand(1, iSides)
Next

Return iTotal

End

Output: 

Tested 100000 times
Player1 with 9 dice of 4 sides
Player2 with 6 dice of 6 sides
Total wins Player1 = 56980 - 0.62823797823798
Total wins Player2 = 35797
7223 draws

Tested 100000 times
Player1 with 5 dice of 10 sides
Player2 with 6 dice of 7 sides
Total wins Player1 = 64276 - 0.48548447320928
Total wins Player2 = 31205
4519 draws

Go

Translation of: C
package main

import(
    "math"
    "fmt"
)

func minOf(x, y uint) uint {
    if x < y {
        return x
    }
    return y
}

func throwDie(nSides, nDice, s uint, counts []uint) {
    if nDice == 0 {
        counts[s]++
        return
    }
    for i := uint(1); i <= nSides; i++ {
        throwDie(nSides, nDice - 1, s + i, counts)
    }
}

func beatingProbability(nSides1, nDice1, nSides2, nDice2 uint) float64 {
    len1 := (nSides1 + 1) * nDice1
    c1 := make([]uint, len1)  // all elements zero by default
    throwDie(nSides1, nDice1, 0, c1)

    len2 := (nSides2 + 1) * nDice2
    c2 := make([]uint, len2)
    throwDie(nSides2, nDice2, 0, c2)
    p12 := math.Pow(float64(nSides1), float64(nDice1)) *
           math.Pow(float64(nSides2), float64(nDice2))

    tot := 0.0
    for i := uint(0); i < len1; i++ {
        for j := uint(0); j < minOf(i, len2); j++ {
            tot += float64(c1[i] * c2[j]) / p12
        }
    }
    return tot
}

func main() {
    fmt.Println(beatingProbability(4, 9, 6, 6))
    fmt.Println(beatingProbability(10, 5, 7, 6))
}
Output:
0.5731440767829815
0.6427886287176273

More idiomatic go: 

package main

import (
	"fmt"
	"math/rand"
)

type set struct {
	n, s int
}

func (s set) roll() (sum int) {
	for i := 0; i < s.n; i++ {
		sum += rand.Intn(s.s) + 1
	}
	return
}

func (s set) beats(o set, n int) (p float64) {
	for i := 0; i < n; i++ {
		if s.roll() > o.roll() {
			p = p + 1.0
		}
	}
	p = p / float64(n)
	return
}

func main() {
	fmt.Println(set{9, 4}.beats(set{6, 6}, 1000))
	fmt.Println(set{5, 10}.beats(set{6, 7}, 1000))
}
Output:
0.576
0.639

Haskell

Translation of: Python
import Control.Monad (replicateM)
import Data.List (group, sort)

succeeds :: (Int, Int) -> (Int, Int) -> Double
succeeds p1 p2 =
  sum
    [ realToFrac (c1 * c2) / totalOutcomes
    | (s1, c1) <- countSums p1 
    , (s2, c2) <- countSums p2 
    , s1 > s2 ]
  where
    totalOutcomes = realToFrac $ uncurry (^) p1 * uncurry (^) p2
    countSums (nFaces, nDice) = f [1 .. nFaces]
      where
        f =
          fmap (((,) . head) <*> (pred . length)) .
          group . sort . fmap sum . replicateM nDice

main :: IO ()
main = do
  print $ (4, 9) `succeeds` (6, 6)
  print $ (10, 5) `succeeds` (7, 6)
Output:
0.5731440767829815
0.6427886287176273

J

Solution: 

gen_dict =: (({. , #)/.~@:,@:(+/&>)@:{@:(# <@:>:@:i.)~ ; ^)&x:

beating_probability =: dyad define
 'C0 P0' =. gen_dict/ x
 'C1 P1' =. gen_dict/ y
 (C0 +/@:,@:(>/&:({."1) * */&:({:"1)) C1) % (P0 * P1)
)

Example Usage: 

   10 5 (;x:inv)@:beating_probability 7 6
┌─────────────────────┬────────┐
│3781171969r5882450000│0.642789│
└─────────────────────┴────────┘
   4 9 (;x:inv)@:beating_probability 6 6
┌─────────────────┬────────┐
│48679795r84934656│0.573144│
└─────────────────┴────────┘

gen_dict explanation:
gen_dict is akin to gen_dict in the python solution and returns a table and total number of combinations, order matters. The table has 2 columns. The first column is the pip count on all dice, the second column is the number of ways this many pips can occur.

({. , #)/.~ make a vector having items head and tally of each group of like items in this case pip count and occurrences operating on
, raveled data (a vector) made of
+/&> the sum of each of the
{ Cartesian products of the
(# <@:>:@:i.) equi-probable die values
&x: but first use extended integers
; ^ links the total possibilities to the result.


The verb beating_probability is akin to the python solution function having same name.
C0 >/&:({."1) C1 is a binary table where the pips of first player exceed pips of second player. "Make a greater than table but first take the head of each item."
C0 */&:({:"1) C1 is the corresponding table of occurrences
* naturally we multiply the two tables (atom by atom, not a matrix product)
+/@:,@: sum the raveled table
% (P0 * P1) after which divide by the all possible rolls.

Java 

import java.util.Random;

public class Dice{
	private static int roll(int nDice, int nSides){
		int sum = 0;
		Random rand = new Random();
		for(int i = 0; i < nDice; i++){
			sum += rand.nextInt(nSides) + 1;
		}
		return sum;
	}
	
	private static int diceGame(int p1Dice, int p1Sides, int p2Dice, int p2Sides, int rolls){
		int p1Wins = 0;
		for(int i = 0; i < rolls; i++){
			int p1Roll = roll(p1Dice, p1Sides);
			int p2Roll = roll(p2Dice, p2Sides);
			if(p1Roll > p2Roll) p1Wins++;
		}
		return p1Wins;
	}
	
	public static void main(String[] args){
		int p1Dice = 9; int p1Sides = 4;
		int p2Dice = 6; int p2Sides = 6;
		int rolls = 10000;
		int p1Wins = diceGame(p1Dice, p1Sides, p2Dice, p2Sides, rolls);
		System.out.println(rolls + " rolls, p1 = " + p1Dice + "d" + p1Sides + ", p2 = " + p2Dice + "d" + p2Sides);
		System.out.println("p1 wins " + (100.0 * p1Wins / rolls) + "% of the time");
		
		System.out.println();
		
		p1Dice = 5; p1Sides = 10;
		p2Dice = 6; p2Sides = 7;
		rolls = 10000;
		p1Wins = diceGame(p1Dice, p1Sides, p2Dice, p2Sides, rolls);
		System.out.println(rolls + " rolls, p1 = " + p1Dice + "d" + p1Sides + ", p2 = " + p2Dice + "d" + p2Sides);
		System.out.println("p1 wins " + (100.0 * p1Wins / rolls) + "% of the time");
		
		System.out.println();
		
		p1Dice = 9; p1Sides = 4;
		p2Dice = 6; p2Sides = 6;
		rolls = 1000000;
		p1Wins = diceGame(p1Dice, p1Sides, p2Dice, p2Sides, rolls);
		System.out.println(rolls + " rolls, p1 = " + p1Dice + "d" + p1Sides + ", p2 = " + p2Dice + "d" + p2Sides);
		System.out.println("p1 wins " + (100.0 * p1Wins / rolls) + "% of the time");
		
		System.out.println();
		
		p1Dice = 5; p1Sides = 10;
		p2Dice = 6; p2Sides = 7;
		rolls = 1000000;
		p1Wins = diceGame(p1Dice, p1Sides, p2Dice, p2Sides, rolls);
		System.out.println(rolls + " rolls, p1 = " + p1Dice + "d" + p1Sides + ", p2 = " + p2Dice + "d" + p2Sides);
		System.out.println("p1 wins " + (100.0 * p1Wins / rolls) + "% of the time");
	}
}
Output:
10000 rolls, p1 = 9d4, p2 = 6d6
p1 wins 57.56% of the time

10000 rolls, p1 = 5d10, p2 = 6d7
p1 wins 64.28% of the time

1000000 rolls, p1 = 9d4, p2 = 6d6
p1 wins 57.3563% of the time

1000000 rolls, p1 = 5d10, p2 = 6d7
p1 wins 64.279% of the time

JavaScript

simulating

let Player = function(dice, faces) {
  this.dice = dice;
  this.faces = faces;
  this.roll = function() {
    let results = [];
    for (let x = 0; x < dice; x++)
      results.push(Math.floor(Math.random() * faces +1));
    return eval(results.join('+'));
  }
}

function contest(player1, player2, rounds) {
  let res = [0, 0, 0];
  for (let x = 1; x <= rounds; x++) {
    let a = player1.roll(),
        b = player2.roll();
    switch (true) {
      case (a > b): res[0]++; break;
      case (a < b): res[1]++; break;
      case (a == b): res[2]++; break;
    }
  }
  document.write(`
      <p>
        <b>Player 1</b> (${player1.dice} × d${player1.faces}): ${res[0]} wins<br>
        <b>Player 2</b> (${player2.dice} × d${player2.faces}): ${res[1]} wins<br>
        <b>Draws:</b> ${res[2]}<br>
        Chances for Player 1 to win:
        ~${Math.round(res[0] / eval(res.join('+')) * 100)} %
      </p>
  `);
}

let p1, p2;

p1 = new Player(9, 4),
p2 = new Player(6, 6);
contest(p1, p2, 1e6);

p1 = new Player(5, 10);
p2 = new Player(6, 7);
contest(p1, p2, 1e6);
Output:
Player 1 (9 × d4): 572753 wins
Player 2 (6 × d6): 356478 wins
Draws: 70769
Chances for Player 1 to win: ~57 %

Player 1 (5 × d10): 643127 wins
Player 2 (6 × d7): 312151 wins
Draws: 44722
Chances for Player 1 to win: ~64 %

jq

Translation of: Wren
Works with: jq

Works with gojq, the Go implementation of jq 

# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

# Input: an array (aka: counts)
def throwDie($nSides; $nDice; $s):
  if $nDice == 0
  then .[$s] +=  1
  else reduce range(1; $nSides + 1) as $i (.;
    throwDie($nSides; $nDice-1; $s + $i) )
  end ;

def beatingProbability(nSides1; nDice1; nSides2; nDice2):
  def a: [range(0; .) | 0];

  ((nSides1 + 1) * nDice1) as $len1
  | ($len1 | a | throwDie(nSides1; nDice1; 0)) as $c1

  | ((nSides2 + 1) * nDice2) as $len2
  | ($len2 | a | throwDie(nSides2; nDice2; 0)) as $c2

  |((nSides1|power(nDice1)) * (nSides2|power(nDice2))) as $p12

  | reduce range(0; $len1) as $i (0;
      reduce range(0; [$i, $len2] | min) as $j (.;
          . + ($c1[$i] * $c2[$j] / $p12) ) ) ;

beatingProbability(4; 9; 6; 6),
beatingProbability(10; 5; 7; 6)
Output:
0.5731440767829815
0.6427886287176273

Julia

Works with: Julia version 0.6
play(ndices::Integer, nfaces::Integer) = (nfaces, ndices) ∋ 0 ? 0 : sum(rand(1:nfaces) for i in 1:ndices)

simulate(d1::Integer, f1::Integer, d2::Integer, f2::Integer; nrep::Integer=1_000_000) =
    mean(play(d1, f1) > play(d2, f2) for _ in 1:nrep)

println("\nPlayer 1: 9 dices, 4 faces\nPlayer 2: 6 dices, 6 faces\nP(Player1 wins) = ", simulate(9, 4, 6, 6))
println("\nPlayer 1: 5 dices, 10 faces\nPlayer 2: 6 dices, 7 faces\nP(Player1 wins) = ", simulate(5, 10, 6, 7))
Output:
Player 1: 9 dices, 4 faces
Player 2: 6 dices, 6 faces
P(Player1 wins) = 0.572805

Player 1: 5 dices, 10 faces
Player 2: 6 dices, 7 faces
P(Player1 wins) = 0.642727

Kotlin

Translation of: C
// version 1.1.2

fun throwDie(nSides: Int, nDice: Int, s: Int, counts: IntArray) {
    if (nDice == 0) {
        counts[s]++
        return
    }
    for (i in 1..nSides) throwDie(nSides, nDice - 1, s + i, counts)
}

fun beatingProbability(nSides1: Int, nDice1: Int, nSides2: Int, nDice2: Int): Double {
    val len1 = (nSides1 + 1) * nDice1
    val c1 = IntArray(len1)  // all elements zero by default
    throwDie(nSides1, nDice1, 0, c1)

    val len2 = (nSides2 + 1) * nDice2
    val c2 = IntArray(len2)
    throwDie(nSides2, nDice2, 0, c2)

    val p12 = Math.pow(nSides1.toDouble(), nDice1.toDouble()) *
              Math.pow(nSides2.toDouble(), nDice2.toDouble())

    var tot = 0.0
    for (i in 0 until len1) {
        for (j in 0 until minOf(i, len2)) {
            tot += c1[i] * c2[j] / p12
        }
    }
    return tot
}

fun main(args: Array<String>) {
    println(beatingProbability(4, 9, 6, 6))
    println(beatingProbability(10, 5, 7, 6))
}
Output:
0.5731440767829815
0.6427886287176273

Lua

Simulated

local function simu(ndice1, nsides1, ndice2, nsides2)
  local function roll(ndice, nsides)
    local result = 0;
    for i = 1, ndice do
      result = result + math.random(nsides)
    end
    return result
  end
  local wins, coms = 0, 1000000
  for i = 1, coms do
    local roll1 = roll(ndice1, nsides1)
    local roll2 = roll(ndice2, nsides2)
    if (roll1 > roll2) then
      wins = wins + 1
    end
  end
  print("simulated:  p1 = "..ndice1.."d"..nsides1..", p2 = "..ndice2.."d"..nsides2..",  prob = "..wins.." / "..coms.." = "..(wins/coms))
end

simu(9, 4, 6, 6)
simu(5, 10, 6, 7)
Output:
simulated:  p1 = 9d4, p2 = 6d6,  prob = 573372 / 1000000 = 0.573372
simulated:  p1 = 5d10, p2 = 6d7,  prob = 642339 / 1000000 = 0.642339

Computed

local function comp(ndice1, nsides1, ndice2, nsides2)
  local function throws(ndice, nsides)
    local sums = {}
    for i = 1, ndice*nsides do
      sums[i] = 0
    end
    local function throw(ndice, nsides, s)
      if (ndice==0) then
        sums[s] = sums[s] + 1
      else
        for i = 1, nsides do
          throw(ndice-1, nsides, s+i)
        end
      end
      return sums
    end
    return throw(ndice, nsides, 0)
  end
  local p1 = throws(ndice1, nsides1)
  local p2 = throws(ndice2, nsides2)
  local wins, coms = 0, nsides1^ndice1 * nsides2^ndice2
  for k1,v1 in pairs(p1) do
    for k2,v2 in pairs(p2) do
      if (k1 > k2) then
        wins = wins + v1 * v2
      end
    end
  end
  print("computed:  p1 = "..ndice1.."d"..nsides1..", p2 = "..ndice2.."d"..nsides2..", prob = "..wins.." / "..coms.." = "..(wins/coms))
end

comp(9, 4, 6, 6)
comp(5, 10, 6, 7)
Output:
computed:  p1 = 9d4, p2 = 6d6, prob = 7009890480 / 12230590464 = 0.57314407678298
computed:  p1 = 5d10, p2 = 6d7, prob = 7562343938 / 11764900000 = 0.64278862871763

M2000 Interpreter

Using a buffer of random numbers for all games. 


Declare random1 lib "advapi32.SystemFunction036" {long lpbuffer, long length}
Flush
Prob(4, 9, 6, 6, 55555)
Prob(10, 5, 7, 6, 55555)
Sub Prob(face1 as long=4 ,dice1=9, face2 as long=6,dice2=6, games=1000)
	if face1<1 or face1>256 or face2<1 or face2>256 then Error "Faces out of limits"
	profiler
	local m=@PlayAll(dice1,dice2, games), tt=timecount
	print "Time to fill the buffer with random bytes ";str$(tt, "#0.00 ms")
	local s1=dice1-1, s2=dice2-1, k, l, f1=face1/256, f2=face2/256, i, n1=0&, n2=0&
	print "Buffer size ";str$(len(m)/1024,"#0.#");" Kbyte": Refresh
	profiler
	for i=0 to games-1
		long n1=dice1:for j=0to s1{n1+=eval(m, i!n1!j)*f1}
		long n2=dice2:for j=0to s2{n2+=eval(m, i!n2!j)*f2}
		if n1>n2 then k++ else if n1<n2 then l++
	next
	print "Games Total "; games
	print "Player 1 wins ";k; " probability of winning:" ;str$(k/games, "#0.00 %")
	print "Player 2 wins ";l; " probability of winning:" ;str$(l/games, "#0.00 %")
	print "Number of ties ";games-k-l
	print "Execution time ";str$(timecount/1000, "#0.0 s")
End Sub
Function PlayAll(dice1, dice2, games as long)
      if games<1 then error "games<1"
      local onegame
      structure onegame {
            n1 as byte * dice1
            n2 as byte * dice2
      }
      buffer Alfa as onegame*games
      call void random1(alfa(0), len(alfa))
      = Alfa
End Function
Output:
Time to fill the buffer with random bytes 2.00 ms
Buffer size 813.8 Kbyte
Games Total 55555
Player 1 wins 31870 probability of winning:57.37 %
Player 2 wins 19729 probability of winning:35.51 %
Number of ties 3956
Execution time 179.9 s

Time to fill the buffer with random bytes 1.65 ms
Buffer size 596.8 Kbyte
Games Total 55555
Player 1 wins 35441 probability of winning:63.79 %
Player 2 wins 17605 probability of winning:31.69
Number of ties 2509
Execution time 138.7 s

Mathematica / Wolfram Language 

ClearAll[GetProbability]
GetProbability[{dice1_List, n_Integer}, {dice2_List, m_Integer}] := 
 Module[{a, b, lena, lenb},
  a = Tuples[dice1, n];
  a = Plus @@@ a;
  lena = a // Length;
  a = Tally[a];
  a[[All, 2]] /= lena;
  
  b = Tuples[dice2, m];
  b = Plus @@@ b;
  lenb = b // Length;
  b = Tally[b];
  b[[All, 2]] /= lenb;
  
  Total[If[#[[1, 1]] > #[[2, 1]], #[[1, 2]] #[[2, 2]], 0] & /@ 
    Tuples[{a, b}]]
  ]
GetProbability[{Range[4], 9}, {Range[6], 6}]
GetProbability[{Range[10], 5}, {Range[7], 6}]
Output:
48679795/84934656
3781171969/5882450000

Nim

Library: bignum
import bignum
from math import sum

proc counts(ndices, nsides: int): seq[int] =
  result.setlen(ndices * nsides + 1)
  for i in 1..nsides:
    result[i] = 1
  for i in 1..<ndices:
    var c = newSeq[int](result.len)
    for sum in i..(i * nsides):
      for val in 1..nsides:
        inc c[sum + val], result[sum]
    result = move(c)

proc probabilities(counts: seq[int]): seq[Rat] =
  result.setLen(counts.len)
  let total = sum(counts)
  for i, n in counts:
    result[i] = newRat(n, total)

proc beatingProbability(ndices1, nsides1, ndices2, nsides2: int): Rat =
  let counts1 = counts(ndices1, nsides1)
  let counts2 = counts(ndices2, nsides2)
  var p1 = counts1.probabilities()
  var p2 = counts2.probabilities()

  result = newRat(0)
  for sum1 in ndices1..p1.high:
    var p = newRat(0)
    for sum2 in ndices2..min(sum1 - 1, p2.high):
      p += p2[sum2]
    result += p1[sum1] * p

echo beatingProbability(9, 4, 6, 6).toFloat
echo beatingProbability(5, 10, 6, 7).toFloat
Output:
0.5731440767829801
0.6427886287176261

ooRexx

Algorithm

Numeric Digits 30
Call test '9 4 6 6'
Call test '5 10 6 7'
Exit
test:
Parse Arg w1 s1 w2 s2
p1.=0
p2.=0
Call pp 1,w1,s1,p1.,p2.
Call pp 2,w2,s2,p1.,p2.
p2low.=0
Do x=w1 To w1*s1
  Do y=0 To x-1
    p2low.x+=p2.y
    End
  End
pwin1=0
Do x=w1 To w1*s1
  pwin1+=p1.x*p2low.x
  End
Say 'Player 1 has' w1 'dice with' s1 'sides each'
Say 'Player 2 has' w2 'dice with' s2 'sides each'
Say 'Probability for player 1 to win:' pwin1
Say ''
Return

pp: Procedure
/*---------------------------------------------------------------------
* Compute and assign probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
* k=1 sets p1.*, k=2 sets p2.*
*--------------------------------------------------------------------*/
Use Arg k,w,s,p1.,p2.
str=''
cnt.=0
Do wi=1 To w
  str=str||'Do v'wi'=1 To' s';'
  End
str=str||'sum='
Do wi=1 To w-1
  str=str||'v'wi'+'
  End
str=str||'v'w';'
str=str||'cnt.'sum'+=1;'
Do wi=1 To w
  str=str||'End;'
  End
Interpret str
psum=0
Do x=0 To w*s
  If k=1 Then
    p1.x=cnt.x/(s**w)
  Else
    p2.x=cnt.x/(s**w)
  psum+=p1.x
  End
Return
Output:
Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 0.642788628717626159168373721835

Algorithm using rational arithmetic

Numeric Digits 30
Call test '9 4 6 6'
Call test '5 10 6 7'
Exit

test:
Parse Arg w1 s1 w2 s2
p1.='0/1'
p2.='0/1'
Call pp 1,w1,s1,p1.,p2.
Call pp 2,w2,s2,p1.,p2.
p2low.='0/1'
Do x=w1 To w1*s1
  Do y=0 To x-1
    p2low.x=fr_add(p2low.x,p2.y)
    End
  End
pwin1='0/1'
Do x=w1 To w1*s1
  pwin1=fr_add(pwin1,fr_Mult(p1.x,p2low.x))
  End
Say 'Player 1 has' w1 'dice with' s1 'sides each'
Say 'Player 2 has' w2 'dice with' s2 'sides each'
Say 'Probability for player 1 to win:' pwin1
Parse Var pwin1 nom '/' denom
Say '                              ->' (nom/denom)
Say ''
Return

pp: Procedure
/*---------------------------------------------------------------------
* Compute and assign probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
* k=1 sets p1.*, k=2 sets p2.*
*--------------------------------------------------------------------*/
Use Arg k,w,s,p1.,p2.
str=''
cnt.=0
Do wi=1 To w
  str=str||'Do v'wi'=1 To' s';'
  End
str=str||'sum='
Do wi=1 To w-1
  str=str||'v'wi'+'
  End
str=str||'v'w';'
str=str||'cnt.sum+=1;'
Do wi=1 To w
  str=str||'End;'
  End
Interpret str

psum='0/1'
Do x=0 To w*s
  If k=1 Then
    p1.x=cnt.x'/'||(s**w)
  Else
    p2.x=cnt.x'/'||(s**w)
  psum=fr_add(psum,p1.x)
  End
Return


fr_add: Procedure
Parse Arg a,b
parse Var a an '/' az
parse Var b bn '/' bz
rn=an*bz+bn*az
rz=az*bz
res=fr_cancel(rn','rz)
Return res

fr_div: Procedure
Parse Arg a,b
parse Var a an '/' az
parse Var b bn '/' bz
rn=an*bz
rz=az*bn
res=fr_cancel(rn','rz)
Return res

fr_mult: Procedure
Parse Arg a,b
parse Var a an '/' az
parse Var b bn '/' bz
rn=an*bn
rz=az*bz
res=fr_cancel(rn','rz)
Return res

fr_cancel: Procedure
Parse Arg n ',' z
k=ggt(n,z)
Return n%k'/'z%k

ggt: Procedure
/**********************************************************************
* ggt (gcd) Greatest common Divisor
* Recursive procedure as shown in PL/I
**********************************************************************/
Parse Arg a,b
if b = 0 then return abs(a)
return ggt(b,a//b)
Output:
Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 48679795/84934656
                              -> 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 3781171969/5882450000
                              -> 0.642788628717626159168373721834

Algorithm using class fraction

Class definition adapted from Arithmetic/Raional. 

Numeric Digits 50
Call test '9 4 6 6'
Call test '5 10 6 7'
Exit

test:
Parse Arg w1 s1 w2 s2
p1.=.fraction~new(0,1)
p2.=.fraction~new(0,1)
Call pp 1,w1,s1,p1.,p2.
Call pp 2,w2,s2,p1.,p2.
p2low.=.fraction~new(0,1)
Do x=w1 To w1*s1
  Do y=0 To x-1
    p2low.x=p2low.x+p2.y
    End
  End
pwin1=.fraction~new(0,1)
Do x=w1 To w1*s1
  pwin1=pwin1+(p1.x*p2low.x)
  End
Say 'Player 1 has' w1 'dice with' s1 'sides each'
Say 'Player 2 has' w2 'dice with' s2 'sides each'
Say 'Probability for player 1 to win:' pwin1~string
Say '                              ->' pwin1~tonumber
Say ''
Return

pp: Procedure
/*---------------------------------------------------------------------
* Compute and assign probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
* k=1 sets p1.*, k=2 sets p2.*
*--------------------------------------------------------------------*/
Use Arg k,w,s,p1.,p2.
str=''
cnt.=0
Do wi=1 To w
  str=str||'Do v'wi'=1 To' s';'
  End
str=str||'sum='
Do wi=1 To w-1
  str=str||'v'wi'+'
  End
str=str||'v'w';'
str=str||'cnt.sum+=1;'
Do wi=1 To w
  str=str||'End;'
  End
Interpret str

psum=.fraction~new(0,1)
Do x=0 To w*s
  If k=1 Then
    p1.x=.fraction~new(cnt.x,s**w)
  Else
    p2.x=.fraction~new(cnt.x,s**w)
  psum=psum+p1.x
  End
Return

::class fraction inherit orderable
::options Digits 50
::method init
  expose numerator denominator
  use strict arg numerator, denominator = 1
  --Trace ?R
  --if numerator == 0 then denominator = 0
  --else if denominator == 0 then raise syntax 98.900 array("Fraction denominator cannot be zero")

  -- if the denominator is negative, make the numerator carry the sign
  if denominator < 0 then do
      numerator = -numerator
      denominator = - denominator
  end


  -- find the greatest common denominator and reduce to
  -- the simplest form
  gcd = self~gcd(numerator~abs, denominator~abs)

  numerator /= gcd
  denominator /= gcd

-- fraction instances are immutable, so these are
-- read only attributes

-- calculate the greatest common denominator of a numerator/denominator pair
::method gcd private
  use arg x, y
  --Say 'gcd:' x y digits()
  loop while y \= 0
      -- check if they divide evenly
      temp = x // y
      x = y
      y = temp
  end
  return x

-- calculate the least common multiple of a numerator/denominator pair
::method lcm private
  use arg x, y
  return x / self~gcd(x, y) * y

::method abs
  expose numerator denominator
  -- the denominator is always forced to be positive
  return self~class~new(numerator~abs, denominator)

-- convert a fraction to regular Rexx number
::method toNumber
  expose numerator denominator
  if numerator == 0 then return 0
  return numerator/denominator

::method add
  expose numerator denominator
  use strict arg other
  -- convert to a fraction if a regular number
  if \other~isa(.fraction) then other = self~class~new(other, 1)

  multiple = self~lcm(denominator, other~denominator)
  newa = numerator * multiple / denominator
  newb = other~numerator * multiple / other~denominator
  return self~class~new(newa + newb, multiple)

::method times
  expose numerator denominator
  use strict arg other
  -- convert to a fraction if a regular number
  if \other~isa(.fraction) then other = self~class~new(other, 1)
  return self~class~new(numerator * other~numerator, denominator * other~denominator)

-- some operator overrides -- these only work if the left-hand-side of the
-- subexpression is a quaternion
::method "*"
  forward message("TIMES")

::method "+"
  -- need to check if this is a prefix plus or an addition
  if arg() == 0 then
      return self  -- we can return this copy since it is imutable
  else
      forward message("ADD")

::method string
  expose numerator denominator
  if denominator == 1 then return numerator
  return numerator"/"denominator

-- override hashcode for collection class hash uses
::method hashCode
  expose numerator denominator
  return numerator~hashcode~bitxor(numerator~hashcode)

::attribute numerator GET
::attribute denominator GET

::requires rxmath library
Output:
Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 48679795/84934656
                              -> 0.57314407678298008294753086419753086419753086419753

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 3781171969/5882450000
                              -> 0.64278862871762615916837372183358974576919480828566

Test

Result from 10 million tries. 

oid='diet.xxx'; Call sysFileDelete oid
Call test '9  4 6 6'
Call test '5 10 6 7'
Exit
test:
Parse Arg n1 s1 n2 s2
Call o 'Player 1:' n1 'dice with' s1 'sides each'
Call o 'Player 2:' n2 'dice with' s2 'sides each'
cnt1.=0
cnt2.=0
win.=0
nn=10000000
Call time 'R'
Do i=1 To nn
  sum1=sum(n1 s1) ; cnt1.sum1+=1
  sum2=sum(n2 s2) ; cnt2.sum2+=1
  Select
    When sum1>sum2 Then win.1+=1
    When sum1<sum2 Then win.2+=1
    Otherwise           win.0+=1
    End
  End
Call o win.1/nn 'player 1 wins'
Call o win.2/nn 'player 2 wins'
Call o win.0/nn 'draws'
/*
Do i=min(n1,n2) To max(n1*s1,n2*s2)
  Call o right(i,2) format(cnt1.i,7) format(cnt2.i,7)
  End
*/
Call o time('E') 'seconds elapsed'
Return            

sum: Parse Arg n s
sum=0
Do k=1 To n
  sum+=rand(s)
  End
Return sum

rand: Parse Arg n
 Return random(n-1)+1

o:
Say arg(1)
Return lineout(oid,arg(1))
Output:
Player 1: 9 dice with 4 sides each
Player 2: 6 dice with 6 sides each
0.5730344 player 1 wins
0.3562464 player 2 wins
0.0707192 draws
186.794000 seconds elapsed
Player 1: 5 dice with 10 sides each
Player 2: 6 dice with 7 sides each
0.6425906 player 1 wins
0.312976 player 2 wins
0.0444334 draws
149.784000 seconds elapsed

PascalABC.NET

Translation of: Python
##
uses numlibabc;

function combos(sides: array of integer; n: integer): array of biginteger;
begin
  if n = 0 then 
  begin
    result := |1bi|; 
    exit
  end;
  result := arrfill(sides.Max * n + 1, 0bi);
  foreach var v in combos(sides, n - 1) index i do
  begin
    if v = 0 then continue;
    foreach var s in sides do 
      result[i + s] += v;
  end;
end;

procedure winning(sides1: array of integer; n1: integer; sides2: array of integer; n2: integer);
begin
  var p1 := combos(sides1, n1);
  var p2 := combos(sides2, n2);
  var (win, loss, tie) := (0bi, 0bi, 0bi); 
  foreach var x1 in p1 index i do
  begin
    win += x1 * (p2?[:i]).sum;
    tie += x1 * (p2?[i:i + 1]).sum;
    loss += x1 * (p2?[i + 1:]).sum;
  end;
  var s := p1.Sum * p2.Sum;
  var winrat := Frc(win, s);
  var tierat := Frc(tie, s);
  var lossrat := Frc(loss, s);
  println('P1 win: ', winrat.ToReal, winrat);
  println('tie   : ', tierat.ToReal, tierat);
  println('P2 win: ', lossrat.ToReal, lossrat);
  println;
end;

winning(arr(1..4), 9, arr(1..6), 6); 
winning(arr(1..10), 5, arr(1..7), 6); 
//winning(|1, 2, 3, 5, 9|, 700, |1, 2, 3, 4, 5, 6|, 800);
Output:
P1 win:  0.57314407678298 48679795/84934656
tie   :  0.070766169838454 144252007/2038431744
P2 win:  0.356089753378566 725864657/2038431744

P1 win:  0.642788628717626 3781171969/5882450000
tie   :  0.0444960303104999 523491347/11764900000
P2 win:  0.312715340971874 735812943/2352980000

Perl

Translation of: Python
use List::Util qw(sum0 max);

sub comb {
    my ($s, $n) = @_;
    $n || return (1);
    my @r = (0) x ($n - max(@$s) + 1);
    my @c = comb($s, $n - 1);
    foreach my $i (0 .. $#c) {
        $c[$i] || next;
        foreach my $k (@$s) {
            $r[$i + $k] += $c[$i];
        }
    }
    return @r;
}

sub winning {
    my ($s1, $n1, $s2, $n2) = @_;

    my @p1 = comb($s1, $n1);
    my @p2 = comb($s2, $n2);

    my ($win, $loss, $tie) = (0, 0, 0);

    foreach my $i (0 .. $#p1) {
        $win  += $p1[$i] * sum0(@p2[0    .. $i - 1]);
        $tie  += $p1[$i] * sum0(@p2[$i   .. $i    ]);
        $loss += $p1[$i] * sum0(@p2[$i+1 .. $#p2  ]);
    }
    my $total = sum0(@p1) * sum0(@p2);
    map { $_ / $total } ($win, $tie, $loss);
}

print '(', join(', ', winning([1 ..  4], 9, [1 .. 6], 6)), ")\n";
print '(', join(', ', winning([1 .. 10], 5, [1 .. 7], 6)), ")\n";
Output:
(0.57314407678298, 0.070766169838454, 0.356089753378566)
(0.642788628717626, 0.0444960303104999, 0.312715340971874)

Phix

Translation of: Go
with javascript_semantics
function throwDie(integer nSides, nDice, s, sequence counts)
    if nDice == 0 then
        counts[s] += 1
    else
        for i=1 to nSides do
            counts = throwDie(nSides, nDice-1, s+i, counts)
        end for
    end if
    return counts
end function
 
function beatingProbability(integer nSides1, nDice1, nSides2, nDice2)
    integer len1 := (nSides1 + 1) * nDice1,
            len2 := (nSides2 + 1) * nDice2
    sequence c1 = throwDie(nSides1, nDice1, 0, repeat(0,len1)),
             c2 = throwDie(nSides2, nDice2, 0, repeat(0,len2))
    atom p12 := power(nSides1, nDice1) * power(nSides2, nDice2),
         tot := 0.0
    for i=1 to len1 do
        for j=1 to min(i-1,len2) do
            tot += (c1[i] * c2[j]) / p12
        end for
    end for
    return tot
end function
 
printf(1,"%0.16f\n",beatingProbability(4, 9, 6, 6))
printf(1,"%0.16f\n",beatingProbability(10, 5, 7, 6))
Output:

(aside: the following tiny discrepancies are to be expected when using IEEE-754 64/80-bit floats; if you want to read anything into them, it should just be that we are all using the same hardware, and probably showing a couple too many digits of (in)accuracy on 32-bit.)
32 bit, same as Go, Kotlin 

0.5731440767829815
0.6427886287176273

64 bit, same as D, Python[last], Ruby, Tcl 

0.5731440767829801
0.6427886287176262

PL/I

version 1

*process source attributes xref;
 dicegame: Proc Options(main);
 Call test(9, 4,6,6);
 Call test(5,10,6,7);
 test: Proc(w1,s1,w2,s2);
 Dcl (w1,s1,w2,s2,x,y) Bin Fixed(31);
 Dcl p1(100)    Dec Float(18) Init((100)0);
 Dcl p2(100)    Dec Float(18) Init((100)0);
 Dcl p2low(100) Dec Float(18) Init((100)0);
 Call pp(w1,s1,p1);
 Call pp(w2,s2,p2);
 Do x=w1 To w1*s1;
   Do y=0 To x-1;
     p2low(x)+=p2(y);
     End;
   End;
 pwin1=0;
 Do x=w1 To w1*s1;
  pwin1+=p1(x)*p2low(x);
  End;
 Put Edit('Player 1 has ',w1,' dice with ',s1,' sides each')
         (Skip,3(a,f(2)));
 Put Edit('Player 2 has ',w2,' dice with ',s2,' sides each')
         (Skip,3(a,f(2)));
 Put Edit('Probability for player 1 to win: ',pwin1)(Skip,a,f(20,18));
 Put Edit('')(Skip,a);
 End;
 pp: Proc(w,s,p);
 /*--------------------------------------------------------------------
 * Compute and assign probabilities to get a sum x
 * when throwing w dice each having s sides (marked from 1 to s)
 *-------------------------------------------------------------------*/
 Dcl (w,s)    Bin Fixed(31);
 Dcl p(100)   Dec Float(18);
 Dcl cnt(100) Bin Fixed(31);
 Dcl (a(12),e(12),v(12),sum,i,n) Bin Fixed(31);
 a=0;
 e=0;
 Do i=1 To w;
   a(i)=1;
   e(i)=s;
   End;
 n=0;
 cnt=0;
 Do v(1)=a(1) To e(1);
   Do v(2)=a(2) To e(2);
     Do v(3)=a(3) To e(3);
       Do v(4)=a(4) To e(4);
         Do v(5)=a(5) To e(5);
           Do v(6)=a(6) To e(6);
             Do v(7)=a(7) To e(7);
               Do v(8)=a(8) To e(8);
                 Do v(9)=a(9) To e(9);
                   Do v(10)=a(10) To e(10);
                     sum=0;
                     Do i=1 To 10;
                       sum=sum+v(i);
                       End;
                     cnt(sum)+=1;
                     n+=1;
                     End;
                   End;
                 End;
               End;
             End;
           End;
         End;
       End;
     End;
   End;
 Do k=w To w*s;
   p(k)=divide(cnt(k),n,18,16);
   End;
 End;
 End;
Output:
Player 1 has  9 dice with  4 sides each
Player 2 has  6 dice with  6 sides each
Probability for player 1 to win: 0.573013663291931152

Player 1 has  5 dice with 10 sides each
Player 2 has  6 dice with  7 sides each
Probability for player 1 to win: 0.642703175544738770

version 2 using rational arithmetic

*process source attributes xref;
 dgf: Proc Options(main);
 Call test(9, 4,6,6);
 Call test(5,10,6,7);
 test: Proc(w1,s1,w2,s2);
 Dcl (w1,s1,w2,s2,x,y) Dec Float(18);
 Dcl 1 p1(100),
      2 nom      Dec Float(18) Init((100)0),
      2 denom    Dec Float(18) Init((100)1);
 Dcl 1 p2(100),
      2 nom      Dec Float(18) Init((100)0),
      2 denom    Dec Float(18) Init((100)1);
 Dcl 1 p2low(100),
      2 nom      Dec Float(18) Init((100)0),
      2 denom    Dec Float(18) Init((100)1);
 Dcl 1 pwin1,
      2 nom      Dec Float(18) Init(0),
      2 denom    Dec Float(18) Init(1);
 Dcl 1 prod Like pwin1;
 Call pp(w1,s1,p1);
 Call pp(w2,s2,p2);
 Do x=w1 To w1*s1;
   Do y=0 To x-1;
     Call fr_add(p2low(x),p2(y),p2low(x));
     End;
   End;
 Do x=w1 To w1*s1;
  Call fr_mult(p1(x),p2low(x),prod);
  Call fr_add(pwin1,prod,pwin1);
  End;
 Put Edit('Player 1 has ',w1,' dice with ',s1,' sides each')
         (Skip,3(a,f(2)));
 Put Edit('Player 2 has ',w2,' dice with ',s2,' sides each')
         (Skip,3(a,f(2)));
 Put Edit('Probability for player 1 to win: ',
          str(pwin1.nom),'/',str(pwin1.denom))(Skip,4(a));
 Put Edit('                              -> ',
          pwin1.nom/pwin1.denom)(Skip,a,f(20,18));
 Put Edit('')(Skip,a);
 End;

 pp: Proc(w,s,p);
 /*--------------------------------------------------------------------
 * Compute and assign probabilities to get a sum x
 * when throwing w dice each having s sides (marked from 1 to s)
 *-------------------------------------------------------------------*/
 Dcl (w,s)    Dec Float(18);
 Dcl 1 p(100),
      2 nom   Dec Float(18),
      2 denom Dec Float(18);
 Dcl cnt(100) Dec Float(18);
 Dcl (a(12),e(12),v(12),sum,i,n) Dec Float(18);
 a=0;
 e=0;
 Do i=1 To w;
   a(i)=1;
   e(i)=s;
   End;
 n=0;
 cnt=0;
 Do v(1)=a(1) To e(1);
   Do v(2)=a(2) To e(2);
     Do v(3)=a(3) To e(3);
       Do v(4)=a(4) To e(4);
         Do v(5)=a(5) To e(5);
           Do v(6)=a(6) To e(6);
             Do v(7)=a(7) To e(7);
               Do v(8)=a(8) To e(8);
                 Do v(9)=a(9) To e(9);
                   Do v(10)=a(10) To e(10);
                     sum=0;
                     Do i=1 To 10;
                       sum=sum+v(i);
                       End;
                     cnt(sum)+=1;
                     n+=1;
                     End;
                   End;
                 End;
               End;
             End;
           End;
         End;
       End;
     End;
   End;
 Do k=w To w*s;
   p(k).nom=cnt(k);
   p(k).denom=n;
   End;
 End;

 fr_add: Proc(a,b,res);
 Dcl 1 a,
      2 nom   Dec Float(18),
      2 denom Dec Float(18);
 Dcl 1 b Like a;
 Dcl res like a;
 /* Put Edit('fr_add',a.nom,a.denom,b.nom,b.denom)(Skip,a,4(f(15))); */
 res.nom=a.nom*b.denom+b.nom*a.denom;
 res.denom=a.denom*b.denom;
 Call fr_cancel(res,res);
 End;

 fr_mult: Proc(a,b,res);
 Dcl 1 a,
      2 nom   Dec Float(18),
      2 denom Dec Float(18);
 Dcl 1 b Like a;
 Dcl res like a;
 /* Put Edit('fr_mult',a.nom,a.denom,b.nom,b.denom)(Skip,a,4(f(15)));*/
 res.nom=a.nom*b.nom;
 res.denom=a.denom*b.denom;
 Call fr_cancel(res,res);
 End;

 fr_cancel: Proc(a,res);
 Dcl 1 a,
      2 nom   Dec Float(18),
      2 denom Dec Float(18);
 Dcl k Dec Float(18);
 Dcl 1 res like a;
 /* Put Edit('fr_cancel',a.nom,a.denom)(Skip,a,4(f(15)));            */
 k=ggt(a.nom,a.denom);
 res=a/k;
 End;

 ggt: Proc(a,b) Recursive Returns(Dec Float(18));
 /**********************************************************************
 * ggt (gcd) Greatest common Divisor
 * Recursive Proc(a,b)) as shown in PL/I
 **********************************************************************/
 Dcl (a,b) Dec Float(18);
 if b = 0 then return (abs(a));
 return (ggt(b,mod(a,b)));
 End;

 str: Proc(x) Returns(Char(20) Var);
 Dcl x Dec Float(18);
 Dcl res Char(20) Var;
 Put String(res) Edit(x)(f(20));
 Return (trim(res));
 End;

 End;
Output:
Player 1 has  9 dice with  4 sides each
Player 2 has  6 dice with  6 sides each
Probability for player 1 to win: 48679795/84934656
                              -> 0.573144076782980083

Player 1 has  5 dice with 10 sides each
Player 2 has  6 dice with  7 sides each
Probability for player 1 to win: 3781171969/5882450000
                              -> 0.642788628717626159

Python 

from itertools import product

def gen_dict(n_faces, n_dice):
    counts = [0] * ((n_faces + 1) * n_dice)
    for t in product(range(1, n_faces + 1), repeat=n_dice):
        counts[sum(t)] += 1
    return counts, n_faces ** n_dice

def beating_probability(n_sides1, n_dice1, n_sides2, n_dice2):
    c1, p1 = gen_dict(n_sides1, n_dice1)
    c2, p2 = gen_dict(n_sides2, n_dice2)
    p12 = float(p1 * p2)

    return sum(p[1] * q[1] / p12
               for p, q in product(enumerate(c1), enumerate(c2))
               if p[0] > q[0])

print beating_probability(4, 9, 6, 6)
print beating_probability(10, 5, 7, 6)
Output:
0.573144076783
0.642788628718

To handle larger number of dice (and faster in general): 

from __future__ import print_function, division

def combos(sides, n):
    if not n: return [1]
    ret = [0] * (max(sides)*n + 1)
    for i,v in enumerate(combos(sides, n - 1)):
        if not v: continue
        for s in sides: ret[i + s] += v
    return ret

def winning(sides1, n1, sides2, n2):
    p1, p2 = combos(sides1, n1), combos(sides2, n2)
    win,loss,tie = 0,0,0 # 'win' is 1 beating 2
    for i,x1 in enumerate(p1):
        # using accumulated sum on p2 could save some time
        win += x1*sum(p2[:i])
        tie += x1*sum(p2[i:i+1])
        loss+= x1*sum(p2[i+1:])
    s = sum(p1)*sum(p2)
    return win/s, tie/s, loss/s

print(winning(range(1,5), 9, range(1,7), 6))
print(winning(range(1,11), 5, range(1,8), 6)) # this seem hardly fair

# mountains of dice test case
# print(winning((1, 2, 3, 5, 9), 700, (1, 2, 3, 4, 5, 6), 800))
Output:
(0.5731440767829801, 0.070766169838454, 0.3560897533785659)
(0.6427886287176262, 0.044496030310499875, 0.312715340971874)

If we further restrict die faces to be 1 to n instead of arbitrary values, the combo generation can be made much faster: 

from __future__ import division, print_function
from itertools import accumulate # Python3 only

def combos(sides, n):
    ret = [1] + [0]*(n + 1)*sides # extra length for negative indices
    for p in range(1, n + 1):
        rolling_sum = 0
        for i in range(p*sides, p - 1, -1):
            rolling_sum += ret[i - sides] - ret[i]
            ret[i] = rolling_sum
        ret[p - 1] = 0
    return ret

def winning(d1, n1, d2, n2):
    c1, c2 = combos(d1, n1), combos(d2, n2)
    ac = list(accumulate(c2 + [0]*(len(c1) - len(c2))))

    return sum(v*a for  v,a in zip(c1[1:], ac)) / (ac[-1]*sum(c1))


print(winning(4, 9, 6, 6))
print(winning(5, 10, 6, 7))

#print(winning(6, 700, 8, 540))
Output:
0.5731440767829801
0.6427886287176262

R

Solving these sorts of problems by simulation is trivial in R. 

probability <- function(facesCount1, diceCount1, facesCount2, diceCount2)
{
  mean(replicate(10^6, sum(sample(facesCount1, diceCount1, replace = TRUE)) > sum(sample(facesCount2, diceCount2, replace = TRUE))))
}
cat("Player 1's probability of victory is", probability(4, 9, 6, 6),
    "in the first game and", probability(10, 5, 7, 6), "in the second.")
Output:
Player 1's probability of victory is 0.572652 in the first game and 0.642817 in the second.

Racket 

#lang racket

(define probs# (make-hash))

(define (NdD n d)
  (hash-ref!
   probs# (cons n d)
   (λ ()
     (cond
       [(= n 0) ; every chance of nothing!
        (hash 0 1)]
       [else
        (for*/fold ((hsh (hash))) (((i p) (in-hash (NdD (sub1 n) d))) (r (in-range 1 (+ d 1))))
          (hash-update hsh (+ r i) (λ (p+) (+ p+ (/ p d))) 0))]))))

(define (game-probs N1 D1 N2 D2)
  (define P1 (NdD N1 D1))
  (define P2 (NdD N2 D2))  
  (define-values (W D L)
    (for*/fold ((win 0) (draw 0) (lose 0)) (((r1 p1) (in-hash P1)) ((r2 p2) (in-hash P2)))
      (define p (* p1 p2))
      (cond
        [(< r1 r2) (values win draw (+ lose p))]
        [(= r1 r2) (values win (+ draw p) lose)]
        [(> r1 r2) (values (+ win p) draw lose)])))
  
  (printf "P(P1 win): ~a~%" (real->decimal-string W 6))
  (printf "P(draw):   ~a~%" (real->decimal-string D 6))
  (printf "P(P2 win): ~a~%" (real->decimal-string L 6))  
  (list W D L))

(printf "GAME 1 (9D4 vs 6D6)~%")
(game-probs 9 4 6 6)
(newline)

(printf "GAME 2 (5D10 vs 6D7) [what is a D7?]~%")
(game-probs 5 10 6 7)
Output:
GAME 1 (9D4 vs 6D6)
P(P1 win): 0.573144
P(draw):   0.070766
P(P2 win): 0.356090
(48679795/84934656 144252007/2038431744 725864657/2038431744)

GAME 2 (5D10 vs 6D7) [what is a D7?]
P(P1 win): 0.642789
P(draw):   0.044496
P(P2 win): 0.312715
(3781171969/5882450000 523491347/11764900000 735812943/2352980000)

Raku

(formerly Perl 6)

Works with: Rakudo version 2020.08.1
sub likelihoods ($roll) {
    my ($dice, $faces) = $roll.comb(/\d+/);
    my @counts;
    @counts[$_]++ for [X+] |(1..$faces,) xx $dice;
    return [@counts[]:p], $faces ** $dice;
}
 
sub beating-probability ([$roll1, $roll2]) {
    my (@c1, $p1) := likelihoods $roll1;
    my (@c2, $p2) := likelihoods $roll2;
    my $p12 = $p1 * $p2;
 
    [+] gather for flat @c1 X @c2 -> $p, $q {
	take $p.value * $q.value / $p12 if $p.key > $q.key;
    }
}

# We're using standard DnD notation for dice rolls here.
say .gist, "\t", .raku given beating-probability < 9d4 6d6 >;
say .gist, "\t", .raku given beating-probability < 5d10 6d7 >;
Output:
0.573144077	<48679795/84934656>
0.64278862872	<3781171969/5882450000>

Note that all calculations are in integer and rational arithmetic, so the results in fractional notation are exact.

REXX

version 1

Translation of: ooRexx

(adapted for Classic Rexx) 

/* REXX */
Numeric Digits 30
Call test '9 4 6 6'
Call test '5 10 6 7'
Exit
test:
Parse Arg w1 s1 w2 s2
plist1=pp(w1,s1)
p1.=0
Do x=w1 To w1*s1
  Parse Var plist1 p1.x plist1
  End
plist2=pp(w2,s2)
p2.=0
Do x=w2 To w2*s2
  Parse Var plist2 p2.x plist2
  End
p2low.=0
Do x=w1 To w1*s1
  Do y=0 To x-1
    p2low.x=p2low.x+p2.y
    End
  End
pwin1=0
Do x=w1 To w1*s1
  pwin1=pwin1+p1.x*p2low.x
  End
Say 'Player 1 has' w1 'dice with' s1 'sides each'
Say 'Player 2 has' w2 'dice with' s2 'sides each'
Say 'Probability for player 1 to win:' pwin1
Say ''
Return

pp: Procedure
/*---------------------------------------------------------------------
* Compute and return the probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
*--------------------------------------------------------------------*/
Parse Arg w,s
str=''
cnt.=0
Do wi=1 To w
  str=str||'Do v'wi'=1 To' s';'
  End
str=str||'sum='
Do wi=1 To w-1
  str=str||'v'wi'+'
  End
str=str||'v'w';'
str=str||'cnt.'sum'=cnt.'sum'+1;'
Do wi=1 To w
  str=str||'End;'
  End
Interpret str
psum=0
Do x=0 To w*s
  p.x=cnt.x/(s**w)
  psum=psum+p.x
  End
res=''
Do x=w To s*w
  res=res p.x
  End
Return res
Output:
Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 0.642788628717626159168373721835

version 2

/* REXX */
oid='diet.xxx'; 'erase' oid
Call test '9  4 6 6'
Call test '5 10 6 7'
Exit
test:
Parse Arg n1 s1 n2 s2
Call o 'Player 1:' n1 'dice with' s1 'sides each'
Call o 'Player 2:' n2 'dice with' s2 'sides each'
cnt1.=0
cnt2.=0
win.=0
nn=10000
Call time 'R'
Do i=1 To nn
  sum1=sum(n1 s1) ; cnt1.sum1=cnt1.sum1+1
  sum2=sum(n2 s2) ; cnt2.sum2=cnt2.sum2+1
  Select
    When sum1>sum2 Then win.1=win.1+1
    When sum1<sum2 Then win.2=win.2+1
    Otherwise           win.0=win.0+1
    End
  End
Call o win.1/nn 'player 1 wins'
Call o win.2/nn 'player 2 wins'
Call o win.0/nn 'draws'
/*
Do i=min(n1,n2) To max(n1*s1,n2*s2)
  Call o right(i,2) format(cnt1.i,7) format(cnt2.i,7)
  End
*/
Call o time('E') 'seconds elapsed'
Return

sum: Parse Arg n s
sum=0
Do k=1 To n
  sum=sum+rand(s)
  End
Return sum

rand: Parse Arg n
 Return random(n-1)+1

o:
Say arg(1)
Return lineout(oid,arg(1))
Output:
Player 1: 9 dice with 4 sides each
Player 2: 6 dice with 6 sides each
0.574 player 1 wins
0.3506 player 2 wins
0.0754 draws
0.109000 seconds elapsed
Player 1: 5 dice with 10 sides each
Player 2: 6 dice with 7 sides each
0.6411 player 1 wins
0.3115 player 2 wins
0.0474 draws
0.078000 seconds elapsed

optimized

This REXX version is an optimized and reduced version of the first part of the 1st REXX example. 

/*REXX pgm computes and displays the probabilities of a two─player S─sided, N─dice game.*/
numeric digits 100                               /*increase/decrease to heart's desire. */
call game  9  4, 6  6   /*1st player:  9 dice,  4 sides;   2nd player:  6 dice,  6 sides*/
call game  5 10, 6  7   /* "     "     5   "   10   "       "     "     6   "    7   "  */
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
game: parse arg  w.1  s.1,   w.2  s.2            /*1st player(dice sides), 2nd player···*/
      p.= 0
                 do   j=1  for 2;         @@.j= prob(w.j, s.j)
                   do k=w.j  to w.j*s.j;  parse var  @@.j   p.j.k  @@.j;  end  /*k*/
                 end   /*j*/
      low.= 0
                 do   j=w.1  to w.1*s.1
                   do k=0  for j;         low.j= low.j + p.2.k;           end  /*k*/
                 end   /*j*/
      say '   Player  1  has '       w.1       " dice with "       s.1      ' sides each.'
      say '   Player  2  has '       w.2       " dice with "       s.2      ' sides each.'
      winP= 0
                 do   j=w.1  to w.1*s.1;  winP= winP   + p.1.j * low.j
                 end   /*j*/
      say 'The probability for first player to win is '  format(winP*100,,digits()%2) "%."
      say                                        /*                               ↑     */
      return                                     /*show 1/2 of 100 dec. digits────┘     */
/*──────────────────────────────────────────────────────────────────────────────────────*/
prob: procedure; parse arg n,s,,@ $;     #.= 0;              pow= s**n
                                 do j=1  for n;     @= @'DO _'j"=1 FOR" s';';   end  /*j*/
      @= @'_=';                  do k=1  for n-1;   @= @"_"k'+'             ;   end  /*k*/
      interpret  @'_'n";#."_'=#.'_"+1"copies(';END', k)
      ns= n*s;                   do j=0  to ns;     p.j= #.j / pow;             end  /*j*/
                                 do k=n  to ns;     $= $ p.k;                   end  /*k*/
      return $   /* ◄──────────────── probability of 1st player to win, S─sided, N dice.*/
output   when using the default inputs:
   Player  1  has  9  dice with  4  sides each.
   Player  2  has  6  dice with  6  sides each.
The probability for first player to win is  57.31440767829800829475308641975308641975308641975309 %.

   Player  1  has  5  dice with  10  sides each.
   Player  2  has  6  dice with  7  sides each.
The probability for first player to win is  64.27886287176261591683737218335897457691948082856633 %.

Ruby 

def roll_dice(n_dice, n_faces)
  return [[0,1]] if n_dice.zero?
  one  = [1] * n_faces
  zero = [0] * (n_faces-1)
  (1...n_dice).inject(one){|ary,_|
    (zero + ary + zero).each_cons(n_faces).map{|a| a.inject(:+)}
  }.map.with_index(n_dice){|n,sum| [sum,n]}  # sum: total of the faces
end

def game(dice1, faces1, dice2, faces2)
  p1 = roll_dice(dice1, faces1)
  p2 = roll_dice(dice2, faces2)
  p1.product(p2).each_with_object([0,0,0]) do |((sum1, n1), (sum2, n2)), win|
    win[sum1 <=> sum2] += n1 * n2        # [0]:draw, [1]:win, [-1]:lose
  end
end

[[9, 4, 6, 6], [5, 10, 6, 7]].each do |d1, f1, d2, f2|
  puts "player 1 has #{d1} dice with #{f1} faces each"
  puts "player 2 has #{d2} dice with #{f2} faces each"
  win = game(d1, f1, d2, f2)
  sum = win.inject(:+)
  puts "Probability for player 1 to win: #{win[1]} / #{sum}",
       "                              -> #{win[1].fdiv(sum)}", ""
end
Output:
player 1 has 9 dice with 4 faces each
player 2 has 6 dice with 6 faces each
Probability for player 1 to win: 7009890480 / 12230590464
                              -> 0.5731440767829801

player 1 has 5 dice with 10 faces each
player 2 has 6 dice with 7 faces each
Probability for player 1 to win: 7562343938 / 11764900000
                              -> 0.6427886287176262

Rust 

// Returns a vector containing the number of ways each possible sum of face
// values can occur.
fn get_totals(dice: usize, faces: usize) -> Vec<f64> {
    let mut result = vec![1.0; faces + 1];
    for d in 2..=dice {
        let mut tmp = vec![0.0; d * faces + 1];
        for i in d - 1..result.len() {
            for j in 1..=faces {
                tmp[i + j] += result[i];
            }
        }
        result = tmp;
    }
    result
}

fn probability(dice1: usize, faces1: usize, dice2: usize, faces2: usize) -> f64 {
    let totals1 = get_totals(dice1, faces1);
    let totals2 = get_totals(dice2, faces2);
    let mut wins = 0.0;
    let mut total = 0.0;
    for i in dice1..totals1.len() {
        for j in dice2..totals2.len() {
            let n = totals1[i] * totals2[j];
            total += n;
            if j < i {
                wins += n;
            }
        }
    }
    wins / total
}

fn main() {
    println!("{}", probability(9, 4, 6, 6));
    println!("{}", probability(5, 10, 6, 7));
}
Output:
0.5731440767829801
0.6427886287176262

Sidef

Translation of: Python
func combos(sides, n) {
    n || return [1]
    var ret = ([0] * (n*sides.max + 1))
    combos(sides, n-1).each_kv { |i,v|
        v && for s in sides { ret[i + s] += v }
    }
    return ret
}

func winning(sides1, n1, sides2, n2) {
    var (p1, p2) = (combos(sides1, n1), combos(sides2, n2))
    var (win,loss,tie) = (0,0,0)
    p1.each_kv { |i, x|
        win  += x*p2.first(i).sum
        tie  += x*p2.slice(i).first(1).sum
        loss += x*p2.slice(i+1).sum
    }
    [win, tie, loss] »/» p1.sum*p2.sum
}

func display_results(String title, Array res) {
    say "=> #{title}"
    for name, prob in (%w(p₁\ win tie p₂\ win) ~Z res) {
        say "P(#{'%6s' % name}) =~ #{prob.round(-11)} (#{prob.as_frac})"
    }
    print "\n"
}

display_results('9D4 vs 6D6',  winning(range(1, 4), 9, range(1,6), 6))
display_results('5D10 vs 6D7', winning(range(1,10), 5, range(1,7), 6))
Output:
=> 9D4 vs 6D6
P(p₁ win) =~ 0.57314407678 (48679795/84934656)
P(   tie) =~ 0.07076616984 (144252007/2038431744)
P(p₂ win) =~ 0.35608975338 (725864657/2038431744)

=> 5D10 vs 6D7
P(p₁ win) =~ 0.64278862872 (3781171969/5882450000)
P(   tie) =~ 0.04449603031 (523491347/11764900000)
P(p₂ win) =~ 0.31271534097 (735812943/2352980000)

Tcl

To handle the nested loop in NdK, Tcl's metaprogramming abilities are exercised. The goal is to produce a script that looks like: 

foreach d0 {1 2 3 4 5 6} {
  foreach d1 {1 2 3 4 5 6} {
    ...
      foreach dN {1 2 3 4 5 6} {
        dict incr sum [::tcl::mathop::+ $n $d0 $d1 ... $DN]
      }
    ...
  }
}

See the comments attached to that procedure for a more thorough understanding of how that is achieved (with the caveat that $d0..$dN are reversed).

Such metaprogramming is a very powerful technique in Tcl for building scripts where other approaches (in this case, recursion) might not be appealing, and should be in every programmer's toolbox! 

proc range {b} {    ;# a common standard proc:  [range 5] -> {0 1 2 3 4}
    set a 0
    set res {}
    while {$a < $b} {
        lappend res $a
        incr a
    }
    return $res
}

# This proc builds up a nested foreach call, then evaluates it.
#
# The script is built up in $script, starting with the body using "%%" as 
# a placeholder.
#
# For each die, a level is wrapped around it as follows:
#   set script {foreach d0 {1 2 3 4 5 6} $script}
#   set script {foreach d1 {1 2 3 4 5 6} $script}
#
# .. and {$d0 $d1 ..} are collected in the variable $vars, which is used
# to replace "%%" at the end.
#
# The script is evaluated with [try] - earlier Tcl's could use [catch] or [eval]
proc NdK {n {k 6}} {    ;# calculate a score histogram for $n dice of $k faces
    set sum {}
    set script {
        dict incr sum [::tcl::mathop::+ $n %%]   ;# add $n because ranges are 0-based
    }   ;# %% is a placeholder
    set vars ""
    for {set i 0} {$i < $n} {incr i} {
        set script [list foreach d$i [range $k] $script]
        append vars " \$d$i"
    }
    set script [string map [list %% $vars] $script]
    try $script
    return $sum
}

proc win_pr {p1 p2} {    ;# calculate the winning probability of player 1 given two score histograms
    set P 0
    set N 0
    dict for {d1 k1} $p1 {
        dict for {d2 k2} $p2 {
            set k [expr {$k1 * $k2}]
            incr N $k
            incr P [expr {$k * ($d1 > $d2)}]
        }
    }
    expr {$P * 1.0 / $N}
}

foreach {p1 p2} {
    {9 4}   {6 6}
    {5 10}  {6 7}
} {
    puts [format "p1 has %dd%d; p2 has %dd%d" {*}$p1 {*}$p2]
    puts [format " p1 wins with Pr(%s)" [win_pr [NdK {*}$p1] [NdK {*}$p2]]]
}
Output:
p1 has 9d4; p2 has 6d6
 p1 wins with Pr(0.5731440767829801)
p1 has 5d10; p2 has 6d7
 p1 wins with Pr(0.6427886287176262)

Factoring out foreach*

The nested-loop generation illustrated above is useful to factor out as a routine by itself. Here it is abstracted as foreach*, with NdK modified to suit.

I include this to emphasise the importance and power of metaprogramming in a Tcler's toolbox, as well as sharing a useful proc. 

package require Tcl 8.6    ;# for [tailcall] - otherwise use [uplevel 1 $script]
# This proc builds up a nested foreach call, then evaluates it.
#
# this:
#   foreach* a {1 2 3} b {4 5 6} {
#     puts "$a + $b"
#   }
#
# becomes:
#   foreach a {1 2 3} {
#     foreach b {4 5 6} {
#       puts "$a + $b"
#     }
#   }
proc foreach* {args} {
    set script [lindex $args end]
    set args [lrange $args 0 end-1]
    foreach {b a} [lreverse $args] {
        set script [list foreach $a $b $script]
    }
    tailcall {*}$script
}

proc NdK {n {k 6}} {    ;# calculate a score histogram for $n dice of $k faces

    set args {}     ;# arguments to [foreach*]
    set vars {}     ;# variables used in [foreach*] arguments that need to be added to sum
    set sum {}      ;# this will be the result dictionary

    for {set i 0} {$i < $n} {incr i} {
        lappend args d$i [range $k]
        lappend vars "\$d$i"
    }

    set vars [join $vars +]

    # [string map] to avoid "Quoting Hell"
    set script [string map [list %% $vars] {
        dict incr sum [expr {$n + %%}]  ;# $n because [range] is 0-based
    }]

    foreach* {*}$args $script
    return $sum
}

V (Vlang)

Translation of: Go
import math

fn min_of(x int, y int) int {
    if x < y {
        return x
    }
    return y
}
 
fn throw_die(n_sides int, n_dice int, s int, mut counts []int) {
    if n_dice == 0 {
        counts[s]++
        return
    }
    for i := int(1); i <= n_sides; i++ {
        throw_die(n_sides, n_dice - 1, s + i, mut counts)
    }
}
 
fn beating_probability(n_sides1 int, n_dice1 int, n_sides2 int, n_dice2 int) f64 {
    len1 := (n_sides1 + 1) * n_dice1
    mut c1 := []int{len: len1}  // all elements zero by default
    throw_die(n_sides1, n_dice1, 0, mut c1)
 
    len2 := (n_sides2 + 1) * n_dice2
    mut c2 := []int{len: len2}
    throw_die(n_sides2, n_dice2, 0, mut c2)
    p12 := math.pow(f64(n_sides1), f64(n_dice1)) *
           math.pow(f64(n_sides2), f64(n_dice2))
 
    mut tot := 0.0
    for i := int(0); i < len1; i++ {
        for j := int(0); j < min_of(i, len2); j++ {
            tot += f64(c1[i] * c2[j]) / p12
        }
    }
    return tot
}
 
fn main() {
    println(beating_probability(4, 9, 6, 6))
    println(beating_probability(10, 5, 7, 6))
}
Output:
0.57314407678298
0.64278862871763

Wren

Translation of: Kotlin
var throwDie // recursive
throwDie = Fn.new { |nSides, nDice, s, counts|
    if (nDice == 0) {
        counts[s] = counts[s] + 1
        return
    }
    for (i in 1..nSides) throwDie.call(nSides, nDice-1, s + i, counts)
}

var beatingProbability = Fn.new { |nSides1, nDice1, nSides2, nDice2|
    var len1 = (nSides1 + 1) * nDice1
    var c1 = List.filled(len1, 0)
    throwDie.call(nSides1, nDice1, 0, c1)

    var len2 = (nSides2 + 1) * nDice2
    var c2 = List.filled(len2, 0)
    throwDie.call(nSides2, nDice2, 0, c2)

    var p12 = nSides1.pow(nDice1) * nSides2.pow(nDice2)
    var tot = 0
    for (i in 0...len1) {
        for (j in 0...i.min(len2)) {
            tot = tot + c1[i] * c2[j] / p12
        }
    }
    return tot
}

System.print(beatingProbability.call(4, 9, 6, 6))
System.print(beatingProbability.call(10, 5, 7, 6))
Output:
0.57314407678298
0.64278862871763

zkl

Translation of: Python
fcn combos(sides, n){
   if(not n) return(T(1));
   ret:=((0).max(sides)*n + 1).pump(List(),0);
   foreach i,v in (combos(sides, n - 1).enumerate()){
      if(not v) continue;
      foreach s in (sides){ ret[i + s] += v }
   }
   ret
}
 
fcn winning(sides1,n1, sides2,n2){
   p1, p2 := combos(sides1, n1), combos(sides2, n2);
   win,loss,tie := 0,0,0; # 'win' is 1 beating 2
   foreach i,x1 in (p1.enumerate()){
      # using accumulated sum on p2 could save some time
      win += x1*p2[0,i].sum(0);
      tie += x1*p2[i,1].sum(0);  // i>p2.len() but p2[bigi,?]-->[]
      loss+= x1*p2[i+1,*].sum(0);
   }
   s := p1.sum(0)*p2.sum(0);
   return(win.toFloat()/s, tie.toFloat()/s, loss.toFloat()/s);
}
println(winning([1..4].walk(), 9, [1..6].walk(),6));
println(winning([1..10].walk(),5, [1..7].walk(),6)); # this seem hardly fair
Output:
L(0.573144,0.0707662,0.35609)
L(0.642789,0.044496,0.312715)
Retrieved from "https://rosettacode.org/wiki/Dice_game_probabilities?oldid=371096"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information