Braun's player's hard 19 vs. dealer's 10; getting it right

This regards Julian Braun's calculations of how many times a player will be dealt a hand vs. a dealer's upcard. In my main simulation, I run all 340 possible combinations, and for the most part, it gets it right. But hard 19 vs. dealer's 10, I get weird outcomes. So I wrote this code to simulate just this one possibility, and it matches the number Braun published. So here's the code, but it doesn't have a GUI (my main 340 possibility code does), this just runs from the command line. And also, if you're not using a C++11 compatible GNU compiler, just ignore the rest of this. So here's the code, it will deal 15,000,000,000 hands, but it will either match or be within one of Braun's published number of 1448:

#include <vector>
#include <algorithm>
#include <iostream>
#include <random>
#include <chrono>
#include <iterator>

/* #define SHOWSHUFFLE */

// functions
void simulate();
void shuffle(std::vector<unsigned short>&, std::vector<unsigned short>&);
void shuffle(std::vector<unsigned short>&);

int main()
{
    simulate();
    return 0;
}

void simulate()
{
    using std::vector;
    using std::cout;
    using std::endl;

    unsigned long  int numsims = 150000;
    unsigned long  int hard19 = 0; // number of times player's hand is 19
    unsigned short int p1;         // player's first card
    unsigned short int p2;         // player's second card
    unsigned short int d1;         // dealter's first card (is ignored)
    unsigned short int du;         // dealer's up card
    unsigned short int k;          // loop counter for the shuffle
    int s;                         // variable for the randomizer seed
    unsigned long int i;
    unsigned long int j;
    // suits are ignored in blackjack except under special casino promotions
    vector<unsigned short int> shoe =
    {  //Ace Deuce                          Jack Queen King
         1,  2,    3, 4, 5, 6, 7, 8, 9, 10, 10,  10,   10, // spades?
         1,  2,    3, 4, 5, 6, 7, 8, 9, 10, 10,  10,   10, // hearts?
         1,  2,    3, 4, 5, 6, 7, 8, 9, 10, 10,  10,   10, // clubs?
         1,  2,    3, 4, 5, 6, 7, 8, 9, 10, 10,  10,   10  // diamonds?
    };

    shuffle(shoe);

    vector <unsigned short int> discards;

    for(i = 0; i < numsims; i++)
    {
        for(j = 0; j < 100000; j++) // Braun's sim was 100,000 hands
        {
            p1 = shoe.back(); // player's first card
            shoe.pop_back();
            discards.push_back(p1);

            d1 = shoe.back();  //dealer's first card (is ignored)
            shoe.pop_back();
            discards.push_back(d1);
            
            p2 = shoe.back(); // player's second card
            shoe.pop_back();
            discards.push_back(p2);
            
            du = shoe.back();  //dealer's up card
            shoe.pop_back();
            discards.push_back(du);
            
            if(du == 10 && (p1 == 10 && p2 == 9)) {
                hard19++;
            }
            else if(du == 10 && (p1 == 9 && p2 ==10)) {
                hard19++;
            }

            if( shoe.empty() )
            {
#ifdef SHOWSHUFFLE
                cout << "\nSize of shoe is " << shoe.size() << endl;
                cout << "Size of discards is " << discards.size() << endl;
                cout << "\nShuffling, sim number " << j << "\n" << endl;
#endif
                for (k = 0; k < 52; k++)
                {
                    shoe.push_back(discards[k]);
                }
                discards.clear();
#if defined(__GNUC__)
                // MIPSPro (Silicon Graphics compiler) doesn't have this stuff
                s = std::chrono::system_clock::now().time_since_epoch().count();
                shuffle(shoe.begin(),shoe.end(), std::default_random_engine(s));
#endif

#if defined __SGI

#endif

#ifdef SHOWSHUFFLE
                cout << "Size of shoe " << shoe.size() << endl;
                cout << "Size of discards " << discards.size() << "\n" << endl;
#endif
            }
        }
        cout << i << endl;
    }
    cout << (hard19/numsims) << endl;
}

void shuffle(std::vector<unsigned short>& shoe) 
{
    static int s; //seed for the random number generator

    s = std::chrono::system_clock::now().time_since_epoch().count();
    shuffle(shoe.begin(), shoe.end(), std::default_random_engine(s));
}