stc/exp/raphus/bbs.py~

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