What is the most efficient way to simulate a 7-sided die with a 6-sided die? I’ve put some thought into it but I’m not sure I get somewhere specifically.
To create a 7-sided die we can use a rejection technique. 3-bits give uniform 1-8 and we need uniform 1-7 which means that we have to reject 1/8 i.e. 12.5% rejection probability.
To create $n * 7$-sided die rolls we need $\lceil log_2( 7^n ) \rceil$ bits. This means that our rejection probability is $p_r(n)=1-\frac{7^n}{2^{\lceil log_2( 7^n ) \rceil}}$.
It turns out that the rejection probability varies wildly but for $n=26$ we get $p_r(26) = 1 – \frac{7^{26}}{2^{\lceil log_2(7^{26}) \rceil}} = 1-\frac{7^{26}}{2^{73}} \approx 0.6\%$ rejection probability which is quite good. This means that we can generate with good odds 26 7-die rolls out of 73 bits.
Similarly, if we throw a fair die $n$ times we get number from $0…(6^n-1)$ which gives us $\lfloor log_2(6^{n}) \rfloor$ bits by rejecting everything which is above $2^{\lfloor log_2(6^{n}) \rfloor}$. Consequently the rejection probability is $p_r(n)=1-\frac{2^{\lfloor log_2( 6^{n} ) \rfloor}}{6^n}$.
Again this varies wildly but for $n = 53$, we get $p_r(53) = 1-\frac{2^{137}}{6^{53}} \approx 0.2\%$ which is excellent. As a result, we can roll the 6-face die 53 times and get ~137 bits.
This means that we get about $\frac{137}{53} * \frac{26}{73} = 0.9207$ 7-face die rolls out of 6-face die rolls which is close to the optimum $\frac{log 7}{log6} = 0.9208$.
Is there a way to get the optimum? Is there an way to find those $n$ numbers as above that minimize errors? Is there relevant theory I could have a look at?
P.S. Relevant python expressions:
min([ (i, round(1000*(1-( 7**i ) / (2**ceil(log(7**i,2)))) )/10) for i in xrange(1,100)], key=lambda x: x[1])
min([ (i, round(1000*(1- ((2**floor(log(6**i,2))) / ( 6**i )) ) )/10) for i in xrange(1,100)], key=lambda x: x[1])
P.S.2 Thanks to @Erick Wong for helping me get the question right with his great comments.
Related question: Is there a way to simulate any $n$-sided die using a fixed set of die types for all $n$?
Okay, here’s a pretty good edit: really @#$%ing awesome algorithm! Like andre’s, it uses only integer arithmetic. Like Erick Wong’s, it achieves (nearly) the ideal entropy rate. There are three modifications to andre’s algorithm:
Here’s some terse C-ish code:
unsigned Roll() {
static unsigned rMod = 1, rValue = 0;
while (true) {
while (rMod < (unsigned)-1 / 6) {
rValue = rValue * 6 + (unsigned)rand() % 6;
rMod *= 6; }
unsigned unused = rMod % 7;
if (rValue < unused) rMod = unused;
else {
rValue -= unused; rMod -= unused;
unsigned answer = rValue % 7;
rValue /= 7; rMod /= 7;
return answer; } } }
Here’s the output:
0: 1428604207
1: 1428626732
2: 1428565440
3: 1428630701
4: 1428475676
5: 1428538428
6: 1428558816
die.InRollCount: 10860331523
die.OutRollCount: 10000000000
Die roll rate (lower is better):
Actual: 1.0860331523
Ideal: 1.086033132502
And the full C++ test program that generates the above output:
#include <iostream>
#include <cmath>
#include <random>
std::default_random_engine generator;
template<unsigned InMod, unsigned OutMod>
struct ConvertedDie
{
ConvertedDie() :
InRollCount(0),
OutRollCount(0),
_inDie(0, InMod - 1),
_reserveMod(1),
_reserveValue(0) { }
unsigned Next()
{
++OutRollCount;
while (true)
{
// Roll the input die as many times as we can
while (_reserveMod < std::numeric_limits<typeof(_reserveMod)>::max() / InMod)
{
++InRollCount;
_reserveValue = _reserveValue * InMod + _inDie(generator);
_reserveMod *= InMod;
}
const unsigned unusableValues = _reserveMod % OutMod;
// This comparison loses entropy,
// which is minimized by making it fail almost always.
if (_reserveValue < unusableValues)
{
_reserveMod = unusableValues;
}
else
{
_reserveValue -= unusableValues;
_reserveMod -= unusableValues;
const unsigned answer = _reserveValue % OutMod;
_reserveValue /= OutMod;
_reserveMod /= OutMod;
return answer;
}
}
}
unsigned long InRollCount;
unsigned long OutRollCount;
private:
std::uniform_int_distribution<unsigned char> _inDie;
unsigned _reserveMod;
unsigned _reserveValue;
};
using namespace std;
int main()
{
constexpr unsigned inMod = 6;
constexpr unsigned outMod = 7;
ConvertedDie<inMod, outMod> die;
unsigned bins[outMod] = {};
for (long i = 0; i < 10000000000; ++i)
{
++bins[die.Next()];
}
for (int i = 0; i < outMod; ++i)
{
cout << i << ": " << bins[i] << endl;
}
cout << endl;
cout << "die.InRollCount: " << die.InRollCount << endl;
cout << "die.OutRollCount: " << die.OutRollCount << endl;
cout.precision(log10(die.OutRollCount) + 1);
cout << endl << "Die roll rate (lower is better):" << endl;
cout << "Actual: " << (double)die.InRollCount / (double)die.OutRollCount << endl;
cout.precision(log10(die.OutRollCount) + 3);
cout << "Ideal: " << log((double)outMod) / log((double)inMod) << endl;
return 0;
}
Roll the D6 twice. Order the pairs $(1,1), (1,2), \ldots, (6,5)$ and associate them with the set $\{1,2,\ldots,35\}$. If one of these pairs is rolled, take the associated single value and reduce it mod $7$. So far you have a uniform distribution on $1$ through $7$.
If the pair $(6,6)$ is rolled, you are required to start over. This procedure will probabilistically end at some point. Since it has a $\frac{35}{36}$ chance of not requiring a repeat, the expected number of repeats is only $\frac{36}{35}$. More specifically, the number of iterations of this 2-dice-toss process has a geometric distribution with $p=\frac{35}{36}$. So $P(\mbox{requires $n$ iterations})=\frac{35}{36}\left(\frac{1}{36}\right)^{n-1}$
Counting by die rolls instead of iterations, with this method, $$\{P(\mbox{exactly $n$ die rolls are needed})\}_{n=1,2,3,\ldots}=\left\{0,\frac{35}{36},0,\frac{35}{36}\left(\frac{1}{36}\right),0,\frac{35}{36}\left(\frac{1}{36}\right)^{2},0,\frac{35}{36}\left(\frac{1}{36}\right)^{3},\ldots\right\}$$
The method that uses base-6 representations of real numbers (@Erick Wong’s answer) has
$$\{P(\mbox{exactly $n$ die rolls are needed})\}_{n=1,2,3,\ldots}=\left\{0,\frac{5}{6},\frac{5}{6}\left(\frac{1}{6}\right),\frac{5}{6}\left(\frac{1}{6}\right)^{2},\frac{5}{6}\left(\frac{1}{6}\right)^{3},\frac{5}{6}\left(\frac{1}{6}\right)^{4},\ldots\right\}$$
Put another way, let $Q(n)=P(\mbox{at most $n$ die rolls are needed using this method})$ and $R(n)=P(\mbox{at most $n$ die rolls are needed using the base-6 method})$. Then we sum the above term by term, and for ease of comparison, I’ll use common denominators:
$$\begin{align}
\{Q(n)\}_{n=1,2,\ldots} &= \left\{0,\frac{45360}{36^3},\frac{45360}{36^3},\frac{46620}{36^3},\frac{46620}{36^3},\frac{46655}{36^3},\frac{46655}{36^3},\ldots\right\}\\
\{R(n)\}_{n=1,2,\ldots} &= \left\{0,\frac{38880}{36^3},\frac{45360}{36^3},\frac{46440}{36^3},\frac{46620}{36^3},\frac{46650}{36^3},\frac{46655}{36^3},\ldots\right\}\\
\end{align}$$
So on every other $n$, the base-6 method ties with this method, and otherwise this method is “winning”.
EDIT Ah, I first understood the question to be about simulating one D7 roll with $n$ D6 rolls, and minimizing $n$. Now I understand that the problem is about simulating $m$ D7 rolls with $n$ D6 rolls, and minimizing $n$.
So alter this by keeping track of “wasted” random information. Here is the recursion in words. I am sure that it could be coded quite compactly in perl:
Going into the first roll, we will have $6$ uniformly distributed outcomes (and so not enough to choose from $\{1,…,7\}$.) This is the initial condition for the recursion.
Generally, going into the $k$th roll, we have some number $t$ of uniformly distributed outcomes to consider. Partition $\{1,\ldots,t\}$ into $$\{1,\ldots,7; 8,\ldots,14;\ldots,7a;7a+1,\ldots,t\}$$ where $a=t\operatorname{div}7$. Agree that if the outcome from the $k$th roll puts us in $\{1,\ldots,7a\}$, we will consider the result mod $7$ to be a D7 roll result. If the $k$th roll puts us in the last portion of the partition, we must re-roll.
Now, either
we succeeded with finding a D7 roll. Which portion of the partition were we in? We could have uniformly been in the 1st, 2nd, …, or $a$th. The next roll therefore will give us $6a$ options to consider, and the recursion repeats.
we did not find another D7 roll. So our value was among $7a+1, \ldots, t$, which is at most $6$ options. However many options it is, call it $b$ for now ($b=t\operatorname{mod}7$). Going into the next roll, we will have $6b$ options to consider, and the recursion repeats.
In the long run, just skip the binary conversion altogether and go with some form of arithmetic coding: use the $6$-dice rolls to generate a uniform base-$6$ real number in $[0,1]$ and then extract base-$7$ digits from that as they resolve. For instance:
int rand7()
{
static double a=0, width=7; // persistent state
while ((int)(a+width) != (int)a)
{
width /= 6;
a += (rand6()-1)*width;
}
int n = (int)a;
a -= n;
a *= 7; width *= 7;
return (n+1);
}
A test run of $10000$ outputs usually requires exactly $10861$ dice rolls, and occasionally needs one or two more. Note that the uniformity of this implementation is not exact (even if rand6
is perfect) due to floating-point truncation, but should be pretty good overall.
A variation on @steveverrill and @alex.jordan answers. Roll two D6s:
Die 1
1 2 3 4 5 6
===========
1|1 2 3 4 5 6
D 2|1 2 3 4 5 6
i 3|1 2 3 4 5 6
e 4|1 2 3 4 5 6
2 5|1 2 3 4 5 6
6|7 7 7 7 7 X
Apply the following rules:
As mentioned in the comments, the probability of hitting 1,2,3,4,5,6 or 7 appears to be $\frac{5}{35} = \frac{1}{7}$. Taking into account of the reroll, the probability is more than $\frac{5}{36}$. It’s actually an infinite series that converges to $\frac{1}{7}$:
$$\sum_{i=1}^n \frac{5}{36^i} = \frac{5}{36} + \frac{5}{36^2} + \frac{5}{36^3} + … = \frac{1}{7}$$
Cut the field where you will throw the die into seven symmetric regions. Throw the die, ignore what number arises, and look to where it landed. Rethrow for border landings. Most of the time you’ll only need to throw once ðŸ™‚
I would use the following Las Vegas algorithm.
The “%7” below means “the rest of the division with 7”.
int roll_dice7()
{
a0 = roll_dice6();
a1 = roll_dice6();
a = 6*(a1-1) + (a0-1) + 1;
if (a > 35) {
return roll_dice7();
} else {
return (a-1)%7+1;
}
}
If you want to keep as much entropy as you can, you can look at powers of 6 in base 7. For instance $6^{25} = 101620613015632362263436_7$ so you can extract 23 7-sided die rolls from 25 6-sided die rolls although you will have to re-roll about 4% of the time.
Another possibility is to use the die to generate a base 6 decimal (heximal?) If you roll a 1, 2, 3, 4, 5, 6 respectively, write down 1, 2, 3, 4, 5, 0 respectively as the next digit to the right of heximal point. For instance if your first two rolls were 6,2, then the first two digits of your heximal number would be $.02$. After a while you would be able to identify your heximal as being in one of the intervals $(0,\frac17)$ or $(\frac17, \frac27)$,…, or $(\frac67,1)$. At this point you could assign one of your desired 7 outcomes according to the numerator of the right endpoint of the interval you landed in.
I give this software answer here – which represents the process described in the question. It performs “poorly” in the sense of giving 0.9131 7-dice rolls per 6-dice roll BUT has the important property that its mathematics more or less make sense as described in the question. Other, way more efficient solutions like @Erick Wong’s are excellent and generic which is brilliant – but I think they leave aspects of randomness to C’s or machine’s arbitrary characteristics like float truncation. Don’t get me wrong – I’m ready to accept it since it is superior, but if possible I would like to see a bit of mathematics behind them. e.g. how it’s guaranteed (?) that after a few million repeats it won’t get stack-at 0 or it won’t suddenly get severely degraded performance?
import random
from collections import defaultdict
class MetaDice:
def __init__(self):
self.bits = 0
self.r = 0
self.in_rejection_count = 0
self.in_count = 0
self.in_rounds = 0
self.out_rejection_count = 0
self.out_count = 0
def _get_more_6_dice(self):
while True:
self.in_rounds += 1
t = 0
for i in xrange(0,53):
t *= 6
t += random.randint(0,5)
self.in_count += 1
if t < 2**137:
break
self.in_rejection_count += 1
self.r *= 2**137
self.r += t
self.bits += 137
def _give_me_73_bits(self):
if self.bits < 73:
self._get_more_6_dice()
self.bits -= 73
t = self.r % (2**73)
self.r /= (2**73)
return t
def give_me_26_7_dice(self):
self.out_count += 1
while True:
inp = self._give_me_73_bits()
if inp < 7**26:
return inp
self.out_rejection_count += 1
a = MetaDice()
hist = defaultdict(int)
for i in xrange(0, 100000):
v = a.give_me_26_7_dice()
for j in xrange(0,26):
r = v % 7
hist[r] +=1
v /= 7
print hist
print "p_r in %3.2f%%" % (float(100 * a.in_rejection_count) / float(a.in_rounds))
print "p_r out %3.2f%%" % (float(100 * a.out_rejection_count) / float(a.out_count))
print "7-face die rolls out of 6-face die rolls: %3.4f" % (float(26 * a.out_count) / float(a.in_count))
output:
defaultdict(<type 'int'>, {0L: 370900, 1L: 372017, 2L: 371918, 3L: 371621, 4L: 371577, 5L: 370468, 6L: 371499})
p_r in 0.21%
p_r out 0.61%
7-face die rolls out of 6-face die rolls: 0.9131
A general algorithm
Rolls on a D6 with numbers from 0 through 5 is equivalent to the digits in a uniformly random real number in the range $[0;1]$ represented in base 6.
Rolls on a D7 with numbers from 0 through 6 is equivalent to the digits in a uniformly random real number in the range $[0;1]$ represented in base 7.
You can roll the D6 enough times to know what the first digit in base 7 to simulate the first D7. In order to simulate another D7 you roll the D6 enough times to know what the second digit in base 7 is.
This algorithm is optimal in terms of number of times you need to roll the D6. But the calculations you need to perform quickly blow up.
You can always throw away the entire state of the algorithm and start over. You lose some entropy, but the outputs will not get skewed.
Range coding/arithmetic coding
Instead of throwing away the entire state, it can be reduced through some clever rounding, which does not skew the output. And that is essentially what the arithmetic coding/range coding algorithms others have mentioned are doing.
Enumerating a small set of states
One possible variation is to enumerate a set of intermediate states and clearly define, what state is thrown away at each step. And algorithm with 21 intermediate states could work like this:
Rather than simply encoding the 21 possible intermediate states as an integer 1 to 21, it simplifies the algorithm to encode the intermediate state in 2 parts, a “current group” value from 1 to 6, and a “state within the group” value that ranges from 1 up to, at most, the current group value.
The intermediate state is in one of six groups of states, the probability distribution among the groups is not significant, all states within a group must have equal probability.
First group has 1 state, this is the starting state. Second group has 2 states, etc.
At each step you roll one D6, combined with the number of possible states in the current group, this produces from 6 to 36 possible outcomes. Each outcome is assigned a new state and whenever possible, a simulated D7.
In group 1 there is only 1 state, so there are 6 outcomes. A D7 can not be simulated. The 6 outcomes take us to the 6 possible states in group 6 with equal probability.
From group 2 there are 12 outcomes. 7 outcomes produce a simulated D7 and take us back to the starting state. 5 outcomes take us to group 5.
From group 3 there are 18 outcomes. 14 outcomes produce a simulated D7 and take us to a state in group 2. 4 outcomes take us to a state in group 4.
From group 4 there are 24 outcomes. 21 outcomes produce a simulated D7 and take us to a state in group 3. 3 outcomes take us to a state in group 3.
From group 5 there are 30 outcomes. 28 outcomes produce a simulated D7 and take us to a state in group 4. 2 outcomes take us to a state in group 2.
from group 6 there are 36 outcomes. 35 outcomes produce a simulated D7 and take us to a state in group 5. 1 outcome takes us back to the starting state.
I thought the algorithm described by kasperd was pretty clever, so I went ahead and implemented it.
(Still learning Python, so I’d be happy for people to point out how this can be made easier to read and more Pythonic).
I’m pretty sure this is effectively the same algorithm as the one described by alex.jordan .
(The test_transitions() function compares the outputs of one transition_a() function that closely follows kasperd’s description, vs. the outputs of another transition() function that closely follows alex.jordan’s description, and verifies that both give identical output for every possible input).
#! /usr/bin/env python
# 2014-08-18: David Cary started.
# licensed under cc by-sa 3.0 with attribution required
def get_d6():
import random
return random.randint( 0, 5 )
def transition_a( numstates, substate, fresh_6_dice_value ):
assert 1 <= numstates <= 6
assert 0 <= substate < numstates
outcome = substate*6 + fresh_6_dice_value
if 1 == numstates :
assert ( 0 <= outcome < 6 )
fresh_7_dice_value = None
numstates=6
substate = outcome
elif 2 == numstates :
assert ( 0 <= outcome < 12 )
if outcome < 7:
fresh_7_dice_value = outcome
numstates = 1
substate = 0
else:
fresh_7_dice_value = None
numstates = 12 - 7
substate = outcome-7
elif 3 == numstates :
assert ( 0 <= outcome < 18 )
if outcome < 14:
fresh_7_dice_value = outcome % 7 # floor division
numstates = 2
substate = outcome // 7
else:
fresh_7_dice_value = None
numstates = 18 - 14
substate = outcome-14
elif 4 == numstates :
assert ( 0 <= outcome < 24 )
if outcome < 21:
fresh_7_dice_value = outcome % 7 # floor division
numstates = 3
substate = outcome // 7
else:
fresh_7_dice_value = None
numstates = 24 - 21
substate = outcome - 21
elif 5 == numstates :
assert ( 0 <= outcome < 30 )
if outcome < 28:
fresh_7_dice_value = outcome % 7 # floor division
numstates = 4
substate = outcome // 7
else:
fresh_7_dice_value = None
numstates = 30 - 28
substate = outcome - 28
elif 6 == numstates :
assert ( 0 <= outcome < 36 )
if outcome < 35:
fresh_7_dice_value = outcome % 7 # floor division
numstates = 5
substate = outcome // 7
else:
fresh_7_dice_value = None
numstates = 36 - 35
substate = outcome - 35
else:
print( "this should never happen." )
assert 1 <= numstates <= 6
assert 0 <= substate < numstates
return numstates, substate, fresh_7_dice_value
def transition( numstates, substate, fresh_6_dice_value ):
assert 1 <= numstates <= 6
assert 0 <= substate < numstates
outcome = substate*6 + fresh_6_dice_value
num_outcomes = numstates*6
assert ( 0 <= outcome < num_outcomes )
max_sevens = num_outcomes // 7 # floor division
if outcome < (7*max_sevens):
# actual_newvalue has equal probability of any value
# from 0 to 7*x-1,
# which is the same as all possible pairs of
# 1d7 (fresh_7_dice_value) and 1d(x) (new_group).
# Store the 1d(x) value for later use, and emit the 1d7 value.
fresh_7_dice_value = outcome % 7
numstates = max_sevens
substate = outcome // 7 # floor division
"""
# (alternate method of partitioning)
fresh_7_dice_value = outcome // max_sevens # floor division
numstates = max_sevens
substate = outcome % max_sevens
"""
else:
# actual_newvalue has equal probability of any value
# from n=7*x to max_newvalue,
# this always happens when 0 == n:
fresh_7_dice_value = None
numstates = num_outcomes - (7*max_sevens)
assert( num_outcomes % 7 == numstates )
substate = outcome - (7*max_sevens)
assert( outcome % 7 == substate )
# print( "-->", numstates, substate, fresh_7_dice_value )
assert 1 <= numstates <= 6
assert 0 <= substate < numstates
return numstates, substate, fresh_7_dice_value
def test_dice_converter():
input_6_dice = 0
output_7_dice = 0
output_passes = 0
numstates = 1
substate = 0
d6values = [0, 0, 0, 0, 0, 0]
d7values = [0, 0, 0, 0, 0, 0, 0]
substate_values_pass = [0, 0, 0, 0, 0, 0, 0]
substate_values_out = [0, 0, 0, 0, 0, 0, 0]
consecutive_passes = 0
max_consecutive_passes = 0
# at any one time, 0 <= substate <= numstates.
# Also, at any one time, all values of substate have equal probability.
while( output_7_dice < 10000 ):
fresh_6_dice_value = get_d6()
input_6_dice += 1
d6values[fresh_6_dice_value] += 1
if(consecutive_passes > 30 ):
print( "consecutive_passes: ", consecutive_passes, "numstates: ", numstates, "substates: ", substate, fresh_6_dice_value )
numstates, substate, fresh_7_dice_value = transition(
numstates, substate, fresh_6_dice_value )
if fresh_7_dice_value is None:
output_passes += 1
substate_values_pass[substate] += 1
consecutive_passes += 1
else:
output_7_dice += 1
d7values[fresh_7_dice_value] += 1
substate_values_out[substate] += 1
if consecutive_passes > max_consecutive_passes:
max_consecutive_passes = consecutive_passes
consecutive_passes = 0
# print( fresh_6_dice_value, " ", fresh_7_dice_value )
print( "input_6_dice ", input_6_dice )
print( "output_7_dice ", output_7_dice )
print( "d6 distribution:", d6values )
print( "d7 distribution:", d7values )
print( "passes:", output_passes)
print( "max_consecutive_passes:", max_consecutive_passes)
print( "substate_values_pass", substate_values_pass)
print( "substate_values_out" , substate_values_out)
assert( sum( d6values ) == input_6_dice )
assert( sum( d7values ) == output_7_dice )
assert( output_7_dice + output_passes == input_6_dice )
def main():
test_transitions()
test_dice_converter()
if __name__ == '__main__':
main()
My tests indicate that this algorithm, given a fair 6-sided die, does give evenly-distributed simulated 7-sided die.
...
('input_6_dice ', 13713)
('output_7_dice ', 10000)
('d6 distribution:', [2277, 2312, 2295, 2299, 2291, 2239])
('d7 distribution:', [1425, 1474, 1437, 1444, 1308, 1446, 1466])
('passes:', 3713)
('max_consecutive_passes:', 6)
...
Alas, this requires about 13749 rolls of a 6-sided die in order to get 10000 rolls of the simulated 7-sided die — better than the over 20000 rolls used by some other algorithms, but not as good as the about 10861 rolls that the algorithm from Chris Culter.
Both this algorithm and Culter’s algorithm attempt to preserve some leftover entropy from previous D6 values.
As supercat pointed out, there’s a small bit of entropy loss in this algorithm when we inspect the current state and decide whether or not to emit a D7 value — equivalent to flipping a biased coin, so less than 1 bit of entropy.
This algorithm simply throws away that entropy.
Culter’s algorithm does the same thing, but it adds a little code complexity to arrange things such that the “coin” is much, much more biased, so it loses much less entropy each time.
As others have said, it’s a good idea to roll again on a double six, to give you 35 possibilities instead of 36. The question then becomes how to assign the 35 possibilities to the 7 scores you require. The following formula is reasonably easy to calculate:
(7+Die1-Die2) mod 7
Die 1
1 2 3 4 5 6
===========
Die 1|0 1 2 3 4 5
2 2|6 0 1 2 3 4
3|5 6 0 1 2 3
4|4 5 6 0 1 2
5|3 4 5 6 0 1
6|2 3 4 5 6 X X=roll again
As can be seen, each score from 0 to 6 appears in the table exactly 5 times.
In plain english, this becomes: “Roll the die twice. Subtract the second roll from the first (if this is not possible because it would give a negative number, add 7 to the value of the first roll before subtracting.) If you roll a double six, discard the score and roll again.”
Let the first roll be $r1$ and the second be $r2$. Re-roll on (6,6). Otherwise, use the following formula to get the 7-sided result ($d7$):
$$d7 = ((r1*6 + r2) mod 7) + 1$$
When the number of rolls $n$ is even, this gives you a $1/36^{n/2}$ chance that you need to keep rolling.
Don’t; use a d8 instead, rerolling 8’s (or rolling 1d8-1 and rejecting 0).
No, this doesn’t answer the precise question; yes it gives the same end result as a d7, so it depends on if the question was asked for practical or academic reasons.