stc/exp/raphus/bbs.py

   1
   2
   3 """
   4 Adapted from http://javarng.googlecode.com/svn/trunk/com/modp/random/BlumBlumShub.java
   5
   6 Original license follows:-
   7
   8 Copyright 2005, Nick Galbreath -- nickg [at] modp [dot] com
   9 All rights reserved.
  10
  11 Redistribution and use in source and binary forms, with or without
  12 modification, are permitted provided that the following conditions are
  13 met:
  14
  15    Redistributions of source code must retain the above copyright
  16    notice, this list of conditions and the following disclaimer.
  17
  18    Redistributions in binary form must reproduce the above copyright
  19    notice, this list of conditions and the following disclaimer in the
  20    documentation and/or other materials provided with the distribution.
  21
  22    Neither the name of the modp.com nor the names of its
  23    contributors may be used to endorse or promote products derived from
  24    this software without specific prior written permission.
  25
  26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  27 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  28 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  29 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  30 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  31 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  32 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  33 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  34 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  35 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  36 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  37
  38 This is the standard "new" BSD license:
  39 http://www.opensource.org/licenses/bsd-license.php
  40
  41 """
  42
  43 import primes;
  44 import random;
  45 import sys;
  46 from math import ceil as _ceil, log as _log
  47
  48
  49 class BlumBlumShub(object):
  50
  51     def getPrime(self, bits):
  52         """
  53         Generate appropriate prime number for use in Blum-Blum-Shub.
  54
  55         This generates the appropriate primes (p = 3 mod 4) needed to compute the
  56         "n-value" for Blum-Blum-Shub.
  57         bits - Number of bits in prime
  58
  59         """
  60         while True:
  61             p = primes.bigppr(bits);
  62             if p%4 == 3:
  63                 break
  64         return p
  65
  66     def generateN(self, bits):
  67         """
  68         This generates the "n value" for use in the Blum-Blum-Shub algorithm.
  69         bits - The number of bits of security
  70         """
  71
  72         p = self.getPrime(bits/2)
  73         q = self.getPrime(bits/2)
  74
  75         # make sure p != q (almost always true, but just in case, check)
  76         while p == q:
  77             q = getPrime(self, bits/2);
  78
  79         # print "GenerateN - p = " + repr(p) + ", q = " + repr(q)
  80         return p * q;
  81
  82
  83     def __init__(self, bits):
  84         """
  85         Constructor, specifing bits for n.
  86         bits - number of bits
  87         """
  88         self.n = self.generateN(bits)
  89         # print "n set to " + repr(self.n)
  90         length = self.bitLen(self.n)
  91         seed = random.getrandbits(length)
  92         self.setSeed(seed)
  93
  94     def setN(self, nn):
  95         """
  96         Sets or resets the self.n value
  97         """
  98         self.n = nn
  99
 100     def setSeed(self, seed):
 101         """
 102         Sets or resets the seed value and internal state.
 103         seed -The new seed
 104         """
 105         self.state = seed % self.n
 106
 107
 108
 109
 110     def bitLen(self, x):
 111         " Get the bit lengh of a positive number"
 112         assert(x>0)
 113         q = 0
 114         while x:
 115             q = q+1
 116             x = x>>1
 117         return q
 118
 119     def next(self, numBits):
 120         self.state = (self.state**2) % self.n
 121         return float(self.state)/self.n
 122
 123     def random(self):
 124         self.state = (self.state**2) % self.n
 125         return float(self.state)/self.n
 126
 127
 128     def sample(self, population, k):
 129         """Chooses k unique random elements from a population sequence.
 130
 131         Returns a new list containing elements from the population while
 132         leaving the original population unchanged.  The resulting list is
 133         in selection order so that all sub-slices will also be valid random
 134         samples.  This allows raffle winners (the sample) to be partitioned
 135         into grand prize and second place winners (the subslices).
 136
 137         Members of the population need not be hashable or unique.  If the
 138         population contains repeats, then each occurrence is a possible
 139         selection in the sample.
 140
 141         To choose a sample in a range of integers, use xrange as an argument.
 142         This is especially fast and space efficient for sampling from a
 143         large population:   sample(xrange(10000000), 60)
 144         """
 145
 146         # Sampling without replacement entails tracking either potential
 147         # selections (the pool) in a list or previous selections in a set.
 148
 149         # When the number of selections is small compared to the
 150         # population, then tracking selections is efficient, requiring
 151         # only a small set and an occasional reselection.  For
 152         # a larger number of selections, the pool tracking method is
 153         # preferred since the list takes less space than the
 154         # set and it doesn't suffer from frequent reselections.
 155
 156         n = len(population)
 157         if not 0 <= k <= n:
 158             raise ValueError("sample larger than population")
 159         random = self.random
 160         _int = int
 161         result = [None] * k
 162         setsize = 21        # size of a small set minus size of an empty list
 163         if k > 5:
 164             setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
 165         if n <= setsize or hasattr(population, "keys"):
 166             # An n-length list is smaller than a k-length set, or this is a
 167             # mapping type so the other algorithm wouldn't work.
 168             pool = list(population)
 169             for i in xrange(k):         # invariant:  non-selected at [0,n-i)
 170                 j = _int(random() * (n-i))
 171                 result[i] = pool[j]
 172                 pool[j] = pool[n-i-1]   # move non-selected item into vacancy
 173         else:
 174             try:
 175                 selected = set()
 176                 selected_add = selected.add
 177                 for i in xrange(k):
 178                     j = _int(random() * n)
 179                     while j in selected:
 180                         j = _int(random() * n)
 181                     selected_add(j)
 182                     result[i] = population[j]
 183             except (TypeError, KeyError):   # handle (at least) sets
 184                 if isinstance(population, list):
 185                     raise
 186                 return self.sample(tuple(population), k)
 187         return result
 188
 189
 190 def main():
 191
 192     bbs = BlumBlumShub(159);
 193     bbs.setN(18532395500947174450709383384936679868383424444311405679463280782405796233163977*39688644836832882526173831577536117815818454437810437210221644553381995813014959)
 194     bbs.setSeed(18532395500947174450709383384936679868383424444311)
 195
 196     #print "Generating 10 numbers"
 197     for i in range (0,10):
 198         bbs.next(159)
 199
 200     print bbs.sample(xrange(100),55)
 201
 202     #print bbs.sample(range(10),5)
 203 if __name__ == "__main__":
 204     main()