DZone Snippets: math code

Bit field class

phaana (Irfan Hamid) — Thu, 15 May 2008 11:48:59 GMT

// A bit field is an important data structure in computer science, major uses include memory allocators and Bloom filters.
// The code is in the form of a header file and an implementation file.





// This class is a code snippet, it can be freely used and distributed.

// Author: irfan[dot]hamid[at]gmail[dot]com

//

// Bit fields are a widely used mechanism in computing applications.

// Typical applications include memory allocators and Bloom filters. This class

// provides a generic bit field implementation for an arbitrary number of bits.

//

// Implementation notes:

// The use of C++ allows the clean separation of the data structure from the

// interface. The bit field is implemented as a vector of unsigned int that is

// allocated at object creation.

//

// field: The vector that represents the bit field itself

// bit_count: The total number of bits in the bit field



#ifndef __BITFIELD_H__

#define __BITFIELD_H__

#include 

#define SZ_UINT sizeof(unsigned int)



class Bitfield

{

protected:

	unsigned int *field;

	unsigned int bit_count;



public:

	Bitfield (int bc);

	~Bitfield ();



	int set (int bit);

	int get (int bit);

	int reset (int bit);

};



#endif







#include 

#include "bitfield.h"



// The constructor takes an int param which gives the number of bits in this

// bit field.

Bitfield::Bitfield (int bc)

{

	bit_count = bc;



	// E.g, ceil(257/32*8) = 65, 64 ints fully used, last one partially used

	field = (unsigned int*) malloc(((int) ceil(bc/SZ_UINT*8)));

}



Bitfield::~Bitfield ()

{

	if (field)

		free(field);

}



// This function sets the corresponding bit in the field equal to 1.

//

// Returns: 0 on success, -1 on error.

int Bitfield::set (int bit)

{

	// Sanity check

	if (bit >= bit_count || !field)

		return -1;



	// The correct index into the vector will be given by bit/SZ_UNIT. The

	// index into the correct vector element is given by bit%SZ_UNIT, to

	// achieve this, simply left shift 0x01 the appropriate number of times.

	field[bit/SZ_UINT] |= (0x00000001 << (bit%(SZ_UINT*8)-1));

	return 0;

}



// This function sets the corresponding bit in the field equal to 0.

//

// Returns: 0 on success, -1 on error.

int Bitfield::reset (int bit)

{

	if (bit >= bit_count || !field)

		return -1;

	field[bit/SZ_UINT] &= ~(0x00000001 << (bit%(SZ_UINT*8)-1));

	return 0;

}



// This function returns the value of the corresponding bit.

//

// Returns: The value of the bit, if the field is initialized and bit index is

// within bounds, -1 otherwise.

int Bitfield::get (int bit)

{

	if (bit >= bit_count || !field)

		return -1;

	return (field[bit/SZ_UINT] & (0x00000001 << (bit%(SZ_UINT*8)-1)) ? 1 : 0);

}

Calculate Pi

jakob () — Sat, 26 Apr 2008 14:57:34 GMT

This script presents a way to compute Pi using a statistic method.
See http://en.wikipedia.org/wiki/Computing_π for further algorithms.



import random, math



class calcPi:

    def __init__(self):

        self.times = pow(10,6)

        self.i = 0

        self.isnot = 0



    def IsOnCircle(self,x,y):

        if math.sqrt(x**2+y**2) < 1:

            return True

        else:

            return False

    

    def run(self):

        for x in range(self.times):

            x,y = random.random(),random.random()

            if self.IsOnCircle(x,y):

                self.i+=1

            else:

                self.isnot+=1

    

    def getResults(self):

        return (float(self.i), float(self.isnot))



    def getPi(self):

        self.run()

        r = self.getResults()

        return r[0]/(r[0]+r[1])*4



if __name__ == '__main__':

    print calcPi().getPi()

First n primes in Ruby

ntk () — Thu, 17 Apr 2008 20:32:45 GMT





#!/usr/local/bin/ruby -w



require 'benchmark' 



class Integer    

   def prime?         # cf. http://snippets.dzone.com/posts/show/4636

     n = self.abs()

     return true if n == 2

     return false if n == 1 || n & 1 == 0

     d = n-1

     d >>= 1 while d & 1 == 0

     20.times do                               # 20 = k from above

       a = rand(n-2) + 1

       t = d

       y = ModMath.pow(a,t,n)                  # implemented below

       while t != n-1 && y != 1 && y != n-1

         y = (y * y) % n

         t <<= 1

       end

       return false if y != n-1 && t & 1 == 0

     end

     return true

   end

end

 

module ModMath

   def ModMath.pow(base, power, mod)

     result = 1

     while power > 0

       result = (result * base) % mod if power & 1 == 1

       base = (base * base) % mod

       power >>= 1;

     end

     result

   end

end







class Integer



   def primes   # cf. http://snippets.dzone.com/posts/show/3734



      sieve = []

      3.step(self, 2) { |i| sieve[i] = i }

      sieve[1] = nil

      sieve[2] = 2



      3.step(Math.sqrt(self).floor, 2) do |i| 

         next unless sieve[i]

         (i*i).step(self, i) do |j|

            sieve[j] = nil

         end

      end



      sieve.compact!



   end 





   def primes2       # cf. http://snippets.dzone.com/posts/show/3734



      primes = [nil, nil].concat((2..self).to_a)



      (2 .. Math.sqrt(self)).each do |i|

         next unless primes[i]

            (i*i).step(self, i) do |j|

               primes[j] = nil

            end

      end

	

      primes.compact!



   end



end







class Integer



   @primes_cache ||= 100.primes

   #@primes_cache ||= [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]

   class << self; attr_accessor :primes_cache; end



   @nthprime_limit ||= 5_000_000

   class << self; attr_reader :nthprime_limit; end

   #class << self; attr_accessor :nthprime_limit; end





   def sieve_size

      num = self.to_i.abs

      logn = Math.log(num)

      if num < 15985     # cf. http://primes.utm.edu/howmany.shtml

         ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor

      else

         ( num * (logn + Math.log(logn) - 0.9427) ).floor

      end

   end





   def nthprime



      num = self.to_i.abs



      # set a limit; cf. http://primes.utm.edu/nthprime/: The 5,000,000th prime is 86,028,121.

      raise "#{num}: number too big for Integer#nthprime" if num > Integer.nthprime_limit      



      primes_cache_size = Integer.primes_cache.size



      if num > primes_cache_size



         # optional; cf. Integer#nthprime_add and Integer#nthprime_add_mr below

         if num >= 50_000 && num < 100_000 && (num - primes_cache_size) < 5000 then return num.nthprime_add end

         if num >= 100_000 && num < 200_000 && (num - primes_cache_size) < 15000 then return num.nthprime_add end

         if num >= 200_000 && (num - primes_cache_size) < 35000 then return num.nthprime_add end

         

         logn = Math.log(num)



         if num < 15985       # cf. http://primes.utm.edu/howmany.shtml

            limit = ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor

            Integer.primes_cache = limit.primes

         else

            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor

            Integer.primes_cache = limit.primes

         end



=begin

         elsif num >= 15985 && num <= 1_000_000

            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor

            Integer.primes_cache = limit.primes

         elsif num > 1_000_000

            Integer.primes_cache = 20_000_000.primes

         end

=end



      else

         return Integer.primes_cache.at(num-1)

      end





      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      if Integer.primes_cache.size > 2_500_000

         num.nthprime_add_mr 

      else

         num.nthprime_add   # faster for a (relatively) small prime cache

      end



   end





   def nthprime_add       #  add primes to Integer.primes_cache up to the nth prime



      num = self.to_i.abs



      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop      # avoid a duplicate prime

      else

         i = last_prime_divmod.first + 1

      end





      while Integer.primes_cache.size < num



         n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and

         n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            next_num_sqrt = Math.sqrt(p).floor



            Integer.primes_cache.each do |d| 

               if d > next_num_sqrt 

                  Integer.primes_cache << p

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

                  break

               elsif p % d == 0

                  break

               end  

            end

         end   # each  



      end  # while



      return Integer.primes_cache.at(num-1)



   end





   def nthprime_add_mr       #  add Miller-Rabin primes to Integer.primes_cache up to the nth prime



      num = self.to_i.abs



      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop

      else

         i = last_prime_divmod.first + 1

      end





      while Integer.primes_cache.size < num



         search_next_prime = true



         while search_next_prime



            n1 = 6*i+1       

            n2 = n1+4        

            i += 1



            [n1, n2].each do |p| 



               next if p % 5 == 0 || p % 7 == 0



               if p.prime?

                  Integer.primes_cache << p

                  search_next_prime = false

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

               end

            end

         end



      end  #  first while loop



      return Integer.primes_cache.at(num-1)



   end





#------------------------------------------- some additional prime methods





   def next_primes_in_cache    # next primes in Integer.primes_cache



      search_next_primes = self.to_i.abs



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop

      else

         i = last_prime_divmod.first + 1

      end



      while search_next_primes > 0



         n1 = 6*i+1       

         n2 = n1+4        

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0

            next_num_sqrt = Math.sqrt(p).floor



            Integer.primes_cache.each do |d| 

               if d > next_num_sqrt 

                  Integer.primes_cache << p

                  search_next_primes -= 1

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

                  break

               elsif p % d == 0

                  break

               end  

            end

         end 

      end   # while



      Integer.primes_cache.last(self.to_i.abs)



   end





   def next_mr_prime     # next Miller-Rabin prime

     

      last_num = self.to_i.abs

      last_num_divmod = last_num.divmod(6)



      if last_num_divmod.last == 1

         i = last_num_divmod.first

      else

         i = last_num_divmod.first + 1

      end



      next_prime = nil

      search_next_prime = true



      while search_next_prime



         n1 = 6*i+1       

         n2 = n1+4        

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            if p > last_num && p.prime?

               next_prime = p

               search_next_prime = false 

               break

            end

         end

      end



         next_prime



   end







   def next_mr_primes(n)     # next Miller-Rabin primes



      last_num = self.to_i.abs

      last_num_divmod = last_num.divmod(6)



      if last_num_divmod.last == 1

         i = last_num_divmod.first

      else

         i = last_num_divmod.first + 1

      end



      next_primes = []

      search_next_primes = n.to_i.abs



      while search_next_primes > 0



         n1 = 6*i+1      

         n2 = n1+4       

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            if p > last_num && p.prime?

               next_primes << p

               search_next_primes -= 1 

            end

         end

      end



      next_primes



   end



end







num1 = 10_000

num1 = 1_000

num1 = 5_000



num2 = 210_349

num2 = 100_125

num2 = 55_000



ret1 = nil

ret2 = nil

ret3 = nil

ret4 = nil





Benchmark.bm(16) do |t| 



   t.report("first #{num1} primes: ") do

      ret1 = num1.nthprime_add

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num1} primes: ") do

      ret2 = num1.nthprime_add_mr

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num1} primes: ") do

      ret3 = num1.nthprime

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num2} primes: ") do

      ret4 = num2.nthprime

   end 



end





puts

puts "the #{num1}th prime number: #{ret1}"

puts "the #{num1}th prime number: #{ret2}"

puts "the #{num1}th prime number: #{ret3}"

puts "the #{num2}th prime number: #{ret4}"

puts

puts "the #{num1}th prime: #{Integer.primes_cache.at(num1-1)}"

puts "the #{num2}th prime: #{Integer.primes_cache.at(num2-1)}"

puts

p Integer.primes_cache.first(10)

p Integer.primes_cache.last(10)

p Integer.primes_cache.size

puts





p 15.next_primes_in_cache



puts 594_213.next_mr_prime



p 149_137.next_mr_primes(5)



puts

p 1_000.sieve_size

p 1_000_000.sieve_size

p 2_500_000.sieve_size

p 5_000_000.sieve_size

p 10_000_000.sieve_size

p 100_000_000.sieve_size





=begin



# check with primegen.c, http://cr.yp.to/primegen.html

primes 2 48611 | nl | tail -n 1

primes 2 104729 | nl | tail -n 1

primes 2 679277 | nl | tail -n 1

primes 2 1301423 | nl | tail -n 1

primes 2 2903603 | nl | tail -n 1

primes 1303129 1303283 | nl

primes 594213 $((594213 + 500)) | head -n 1

primes 149137 $((149137 + 1000)) | head -n 5



# further prime tools from http://libtom.org

- LibTomMath

   bn_mp_prime_is_prime.c

   bn_mp_prime_miller_rabin.c

- TomsFastMath

   fp_isprime.c

   fp_prime_miller_rabin.c



=end

ruby is_numeric? extension to String

sikelianos (Zeke Sikelianos) — Fri, 18 Jan 2008 07:46:24 GMT



class String

  def is_numeric?

    Float self rescue false

  end

end

Pascal's triangle in Ruby

ntk () — Fri, 26 Oct 2007 17:02:28 GMT





# cf. http://www.ruby-forum.com/topic/97105



6.times { k=0; p $*.map!{|i|k+k=i} << 1 }



# ... or ...



ar=[]; 6.times { k=0; p ar.map!{|i|k+k=i} << 1 }





# a more general approach (for polynomials)



class Polynomial



   def triangle(nterms, row, pos=nil)



      return nil if nterms < 2 || row < 1

      nterms = nterms - 2

      num_of_rows = row



      var1 = 0 + nterms   

      var2 = 1 + nterms

      var3 = 3 + nterms 

  

      ar1 = [0, 1, 0]   # first row

      var1.times { ar1.push(0) }

      var1.times { ar1.unshift(0) }



      ar2 = []

      ar3 = []

      ar4 = [[1]]



      for num in 0..(num_of_rows - 1)  



         nextnum = ar1.size - var2



         for nextn in 1..nextnum

            sum = 0

            count = 0

            ar1.each do |n|  

               count += 1 

               if count < var3 then t = sum += n; ar2 << t else break end 

            end



            ar3 << ar2.last

            ar2.clear

            ar1.shift



         end   # second for-loop



         ar1.clear

         ar1 << ar3

         ar1.flatten!



         var2.times { ar1.push(0) }

         var2.times { ar1.unshift(0) }



         ar4 << ar3

         ar3 = []



      end  # first for-loop



      if !pos.nil?

         ret = ar4.at(row).at(pos)

         return "No such position: #{pos} in row: #{row}" if ret.nil?

         ret

      else

         ar4.map! { |r| r.join('-') }

         ar4

      end

   end 

end 





puts Polynomial.new.triangle(2, 5)

puts Polynomial.new.triangle(3, 5)

puts Polynomial.new.triangle(4, 5)

puts Polynomial.new.triangle(5, 5)

puts Polynomial.new.triangle(5, 4, 8)

puts Polynomial.new.triangle(4, 9)

puts Polynomial.new.triangle(4, 9, 10)





#------------------------





class Integer

   def fak

      f=1

      (2..self).each { |i| f *= i   }

      f

   end

end



module Enumerable

   def sum

      inject { |n, m| n + m  }

   end

end



# cf. http://blade.nagaokaut.ac.jp/~sinara/ruby/math/combinatorics/array-rep_perm.rb

class Array

  def rep_perm(n)

    if n < 0

    elsif n == 0

      yield([])

    else

      rep_perm(n - 1) do |x|

	each do |y|

	  yield(x + [y])

	end

      end

    end

  end

end





nterms = 2

exponent = 80

exponent = 8



# create the same number of variable names as there are terms 

# example: ['a', 'b'] for (a+b)**3



var_names = ('a'..'z').to_a.slice(0, nterms)

#var_names = (('a'..'z').to_a << ('A'..'Z').to_a << ('aa'..'zz').to_a).flatten!.slice(0, nterms)



ar1 = []

(0..exponent).to_a.rep_perm(nterms) { |x| p x; ar1 << x if x.sum == exponent }   # example: ... if [2,6].sum == 8

ar1.reverse!



#p ar1

#puts ar1



ar2 = []



ar1.each do |term|



   #puts "term: #{term.inspect}"  # example: term: [5, 0, 0]

   count = 0

   var1 = 1

   term.each { |i| var1 *= i.fak }

   var2 = exponent.fak / var1 

   var3 = "#{var2}   (  #{ term.join('-') << '-' } )"  # prepare term for parsing with gsub below

   ar2 << var3



end



#p ar2



result = ar2.collect do |term|

   p term

   count = -1

   term.gsub!(/(\d+)-/)  { count += 1;  "#{var_names.at(count)}" << '**' << $1 << ' ' }

   term.gsub!(/^(\d+)( +)/, '\1\2*\2')

end



puts result

Random numbers in Ruby

ntk () — Wed, 24 Oct 2007 19:01:38 GMT





# Pseudo random number generator

# From: http://gnuvince.wordpress.com/2005/11/24/pseudo-random-number-generator

# Author: Vincent Foley-Bourgon



def random_number

  t = Time.now.to_f / (Time.now.to_f % Time.now.to_i)

  random_seed = t * 1103515245 + 12345;

  (random_seed / 65536) % 32768;

end



10.times { puts random_number }



cnt = Array.new(10,0)

20000.times { cnt[(random_number % 10).to_i] += 1 }

p cnt

puts





# Ruby Gaussian Random Number Generator

# Author: Glenn

# http://webhost101.net/rails/typo/articles/2007/07/31/ruby-gaussian-random-number-generator



def gaussian_rand 

   u1 = u2 = w = g1 = g2 = 0  # declare

   begin

     u1 = 2 * rand - 1

     u2 = 2 * rand - 1

     w = u1 * u1 + u2 * u2

   end while w >= 1

   

   w = Math::sqrt( ( -2 * Math::log(w)) / w )

   g2 = u1 * w;

   g1 = u2 * w;

   # g1 is returned  

end



sum = 0

sumsq = 0

n = 1000

0.upto(n) do 

  #r = gaussian_rand

  r = gaussian_rand * 100 + 50   # new_random_number = gaussian_rand * standard_deviation + average

  #puts r

  sum += r

  sumsq += (r*r)

end



ave = sum/n

stddev = Math::sqrt(sumsq/n - ave * ave)

puts "Average = #{ave}"

puts "StdDev = #{stddev}"

puts





# Random Number Generator

# http://scutigena.sourceforge.net/test-random.html

# http://scutigena.sourceforge.net/sources/ruby-1.7.2/random.html

# $Id: random.ruby,v 1.2 2003/12/30 01:25:05 davidw Exp $



IM = 139968

IA = 3877

IC = 29573



$last = 42.0

def gen_random (max) (max * ($last = ($last * IA + IC) % IM)) / IM end



10.times do

    puts gen_random(100000.0)

end



printf "%.5f\n", gen_random(100.0)

puts





# From: http://matt.blogs.it/entries/00002641.html

# Author: Matt Mower

# cf. http://www.taygeta.com/random/gaussian.html

# cf. http://www.bearcave.com/misl/misl_tech/wavelets/hurst/random.html



def box_mueller( mean = 0.0, stddev = 1.0 )

  x1 = 0.0, x2 = 0.0, w = 0.0



  until w > 0.0 && w < 1.0

    x1 = 2.0 * rand - 1.0

    x2 = 2.0 * rand - 1.0

    w = ( x1 * x2 ) + ( x2 * x2 )

  end



  w = Math.sqrt( -2.0 * Math.log( w ) / w )

  r = x1 * w



  mean + r * stddev

end



10.times { puts box_mueller(5.0, 1.0) }

puts





# Generating random numbers with a specified distribution

# From: http://www.cs.nyu.edu/~michaels/blog/?p=24

# Author: Michael Schidlowsky



def gen

  (x=rand)>0.5 ? x : (x < rand/2.0) ? 1.0 - x : x

end



def gen2

   (x = rand) && rand ? 1.0 - x : x

end



def compute_distribution(numSamples, &block)

   samples = []

   values = 10

   numSamples.times{samples << yield}

   dist = Array.new(values, 0)

   samples.each{ |s| dist[(s*values).floor] += 1 }

   dist

end



p compute_distribution(1000){rand}

p compute_distribution(1000){gen}

p compute_distribution(1000) { gen2 }

puts





# Random number generation for a triangular distribution

# From: http://www.brighton-webs.co.uk/distributions/triangular.asp

# cf. http://en.wikipedia.org/wiki/Triangular_distribution



def rng(m, low, high)



   # cf. the parameter info at http://www.brighton-webs.co.uk/distributions/triangular.asp

   return nil unless high > low && m > low && m < high    



   u = rand



   if u <= (m-low)/(high-low)

      r = low+ Math.sqrt(u*(high-low)*(m-low))

   else

      r = high - Math.sqrt((1.0-u)*(high-low)*(high-m))

   end



end



10.times do

  #puts rng(0.0, -25.0, 25.0)

  puts rng(50.0, 25.0, 75.0)

end

puts





# From: http://www.ruby-forum.com/topic/95931

# Author: Christoffer Lernö



class Range

  def rand

    return first if exclude_end? && last == first

    Kernel::rand(last - first + (exclude_end? ? 0 : 1)) + first

  end

end



r1 = -9..9

r2 = -90...100

p r1.rand, r2.rand



r3 = 0..9



num = []

10.times do

   num << r3.rand

   if num.first.zero? then num.shift; redo end

end



p num

puts num.to_s.to_i

puts





# See:

# http://redcorundum.blogspot.com/2008/01/randomizing-array-revisited.html

# http://redcorundum.blogspot.com/2006/12/randomizing-array-and-other-faqs.html

# http://szeryf.wordpress.com/2007/06/19/a-simple-shuffle-that-proved-not-so-simple-after-all/



class Array

  def shuffle

    array = dup

    size.downto 2 do |j|

      r = rand j

      array[j-1], array[r] = array[r], array[j-1]

    end

    array

  end

end



array = (1..50).to_a



10.times { p array.shuffle.first(10) }





#-------------------





class RNG



   def num(min=8,max=min+5,iter=1)



      return nil unless min.is_a?(Integer) && max.is_a?(Integer) && max > min && iter.is_a?(Integer)



      ret = []

      stats = Hash.new(0)     # optional; cf. stats[random_num] += 1

      digits = Array(0..9).sort_by {rand}

      #digits = (Array(0..9) * (rand(4)+1)).sort_by {rand}

      digits_size = digits.size



      iter.times do

         count = 0

         len = min + rand(max-min+1)

         ar = []

         while count < len

            i = rand(digits_size)    # get a random array index

            random_num = digits.at(i)

            stats[random_num] += 1

            ar << random_num

            digits = digits.sort_by {rand}

            count += 1

            if count == 1 && ar.first.zero? then ar.shift; stats[0] -= 1; count = 0 end

         end

         ret << ar.to_s.to_i

      end  # iter

      #ret

      ret << stats   # optional

   end



end



puts



min = 8

max = 13

iter = 10000

result =  RNG.new.num(min, max, iter)

stats = result.pop



t = 0

stats.each_value { |v| t += v }  # get the overall frequency



n = 0

t = t * 1.0

m = t / 10.0 



stats.sort.each do |k,v| 

   i = (v / t) * 100.0

   x = v > m ? v - m : m - v

   y = (x / m) * 100.0

   n += i

   puts "#{k}  ::  #{"%.2f" % m}  ::  #{v}  ::  #{"%.2f" % i} %  ::  #{ "%.2f" % ((m-v) * -1) }  ::  #{"%.2f" % y} %"

end



puts n, m, t



puts

puts RNG.new.num(20, 25, 10)[0..-2]



puts

rn = RNG.new.num(80, 100)

puts rn[0..-2]

p rn.last

Miller-Rabin prime test in Ruby

ntk () — Thu, 11 Oct 2007 18:20:02 GMT





# From: http://en.wikipedia.org/wiki/Miller-Rabin_primality_test



class Integer

   def prime?

     n = self.abs()

     return true if n == 2

     return false if n == 1 || n & 1 == 0



     # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and

     # http://everything2.com/index.pl?node_id=1176369



     return false if n > 3 && n % 6 != 1 && n % 6 != 5     # added



     d = n-1

     d >>= 1 while d & 1 == 0

     20.times do                               # 20 = k from above

       a = rand(n-2) + 1

       t = d

       y = ModMath.pow(a,t,n)                  # implemented below

       while t != n-1 && y != 1 && y != n-1

         y = (y * y) % n

         t <<= 1

       end

       return false if y != n-1 && t & 1 == 0

     end

     return true

   end

end

 

module ModMath

   def ModMath.pow(base, power, mod)

     result = 1

     while power > 0

       result = (result * base) % mod if power & 1 == 1

       base = (base * base) % mod

       power >>= 1;

     end

     result

   end

end





#------------------------------------------------------------------





# From: http://www.h2np.net/tips/ruby-cipher-math.txt

# Author: Hironobu SUZUKI



# $Id: mathh.rb,v 1.5 2002/11/24 16:15:56 hironobu Exp $

# Simple Secure Stream Cryptography

# Mathematical Library

# Hironobu SUZUKI 

# LGPL, (C) 2002, Hironobu SUZUKI



class Random  

  def initialize()

    @random_device="/dev/urandom"

    @verbose=false

  end

  def random_device=(str)

    @random_device=str

  end

  def verbose

    @verbose=true

  end

  def get(size)

    if (@verbose) ;STDERR.print "+";end

    value=0

    r=File.open(@random_device)

    r.read(size).each_byte {|c|

      value = (value << 8) + c.to_i

    }

    r.close

    return value

    

  end

  def get_prime(lower,higher)

    bytesize = higher.bytesize

    while true

      r = self.get(bytesize)

      if (r % 2 == 0) ; r -= 1 ;   end

      if ( lower <=  r && r <= higher && r.prime? == true )

	return r

      end

    end

  end

end





class Integer

  def cdiv(d)			# ceil divide

    w=self.divmod(d)

    if w[1] > 0 ;  w[0] +=1 ; end

    return w[0]

  end

  def powm(i,n)	# square-and-multiply for exp modulo n

    if ( i == 0 )

      return 1

    end

    b=1

    a=self



    if ( (i & 1) == 1 )

      b=a

    end

    i = i>>1



    while ( i > 0 )

      a = (a * a) % n

      if ( (i & 1) == 1 )

	b = (a * b) % n 

      end

      i = i>>1

    end

    return b

  end

  def gcd(b)			# GCD

    a = self

    if ( a < b ); t=b; b=a; a=t; end

    while b > 0 

      r = a % b 

      a = b

      b = r

    end

    return a

  end

  def invert(n)  # Extension Euclid Algorithm

    a = n        # See Knuth's The Art of Computer Programming 

    b = self % n # Vol2 pp.342 -- 343  

    p = 0 ; q = 1;  v = 0;  u = 1;

    while q > 0 

      p = a / b

      q =  a % b

      w = v - (p*u)

      ## DEBUG    printf "%d = %d(%d) + %d  (%d)\n", a,p,b,q,v ###

      a = b

      b = q

      v = u

      u = w

    end

    return (v+n) % n

  end



# Miller-Rabin Test  (Prime Test)

# See, http://www.cs.albany.edu/~berg/csi445/Assignments/pa4.html

  def bitarray(n) 

    b=Array::new  

    i=0       

    v=n

    while v > 0

      b[i] = (0x1 & v)

      v = v/2

      i = i + 1

    end

    return b

  end

  def miller_rabin(n,s)

    b=bitarray(n-1)

    i=b.size 

    j =1

    while j <= s

      a = 1 + (rand(n-2).to_i)

      if witness(a,n,i,b) == true

	return false

      end

      j+=1

    end

    return true

  end

  def witness(a,n,i,b)

    d=1

    x=0

    while i > 0 

      x = d 

      d = (d**2) % n

      if ( (d == 1) && (x != 1) && (x != (n-1)) )

	return true

      end

      if ( b[i-1] == 1 )

	d = (d * a ) % n

      end

      i -= 1

    end

    if ( d != 1) 

      return true

    end

    return false

  end



  #def prime?

  def is_prime?

    a = self

    return miller_rabin(a,30)    

  end



  def bytesize

    n = self

    i=0

    while n > 0

      n = n >> 8

      i += 1

    end

    return i

  end

  def bitsize

    n = self

    i=0

    while n > 0

      n = n >> 1

      i += 1

    end

    return i

  end

end





p 107565456790871.prime?      # true

p 107565456790871.is_prime?   # true



a = 107565456790871



begin

  a += 2

end while !a.prime? && !a.is_prime?



puts a   #=> 107565456790991 (next prime)

nsieve for Ruby

ntk () — Thu, 11 Oct 2007 16:19:39 GMT

From: The Great Computer Language Shootout: nsieve-bits Ruby
Author: Glenn Parker





CharExponent = 3

BitsPerChar = 1 << CharExponent

LowMask = BitsPerChar - 1



def sieve(m)

  ret = []

  items = "\xFF" * ((m / BitsPerChar) + 1)

  masks = ""

  BitsPerChar.times do |b|

    masks << (1 << b).chr

  end



  count = 0

  pmax = m - 1

  2.step(pmax, 1) do |p|

    if items[p >> CharExponent][p & LowMask] == 1

      ret << p

      count += 1

      p.step(pmax, p) do |mult|

	a = mult >> CharExponent

	b = mult & LowMask

	items[a] -= masks[b] if items[a][b] != 0

        #p items

      end

    end

  end

  #count

  ret

end



n = (ARGV[0] || 2).to_i

n.step(n - 2, -1) do |exponent|

  break if exponent < 0

  m = 2 ** exponent * 10_000

  primes = sieve(m)

  count = primes.size

  puts

  printf "Primes up to %8d %8d\n", m, count

  puts "last ten primes: #{primes[-10..-1].join(' ')}"

  puts

end

Custom base conversion (to and from)

devijvers (Steven Devijver) — Sun, 30 Sep 2007 23:44:48 GMT

Handy to store large numbers in a shorter notation. Doesn't work with negative numbers (for negative numbers you need to store an additional byte to indicate a positive or negative value, pass the absolute value).



	private static final short BASE_40_RADIX = 40;

	private static final short BASE_40_BASE = 48;

	private static final String BASE_40_ZERO = "0";



	public static long fromBase40(char[] cs) {

		return fromBase(cs, BASE_40_BASE, BASE_40_RADIX);

	}



	public static char[] toBase40(final long i) {

		return toBase(i, BASE_40_BASE, BASE_40_RADIX, BASE_40_ZERO);

	}



	private static final short BASE_75_RADIX = 75;

	private static final short BASE_75_BASE = 48;

	private static final String BASE_75_ZERO = "0";



	public static long fromBase75(char[] cs) {

		return fromBase(cs, BASE_75_BASE, BASE_75_RADIX);

	}



	public static char[] toBase75(final long i) {

		return toBase(i, BASE_75_BASE, BASE_75_RADIX, BASE_75_ZERO);

	}



	private static final short BASE_90_RADIX = 90;

	private static final short BASE_90_BASE = 33;

	private static final String BASE_90_ZERO = "!";



	public static long fromBase90(char[] cs) {

		return fromBase(cs, BASE_90_BASE, BASE_90_RADIX);

	}



	public static char[] toBase90(final long i) {

		return toBase(i, BASE_90_BASE, BASE_90_RADIX, BASE_90_ZERO);

	}



	private static char[] toBase(final long i, final short base, final short radix, String nullCharacter) {

		long value = i;



		if (value == 0) {

			return nullCharacter.toCharArray();

		}



		String result = "";

		while (value > 0) {

			long mod = value % radix;

			value -= mod;

			if (value > 0) value /= radix;

			result = (char)(mod + base) + result;

		}



		return result.toCharArray();

	}



	private static long fromBase(char[] cs, final short base, final short radix) {

		long value = 0;

		for (int i = cs.length - 1; i >= 0; i--) {

			int digit = ((int)cs[i]) - base;

			if (digit < 0 || digit >= radix) {

				throw new IllegalArgumentException("Invalid Base" + radix + " character: " + cs[i]);

			}

			long digitBase = (long)Math.pow(radix, (cs.length - 1) - i);

			value = (digit * digitBase) + value;

		}



		return value;

	}

And here are some tests:



	public void testFromBase40() {

		assertEquals(2559999, NumberUtils.fromBase40("WWWW".toCharArray()));

		assertEquals(1117580, NumberUtils.fromBase40("ABCD".toCharArray()));



		assertEquals(0, NumberUtils.fromBase40("0".toCharArray()));

		assertEquals(1, NumberUtils.fromBase40("1".toCharArray()));



		assertEquals(40, NumberUtils.fromBase40("10".toCharArray()));

		assertEquals(80, NumberUtils.fromBase40("20".toCharArray()));



		assertEquals(41, NumberUtils.fromBase40("11".toCharArray()));

		assertEquals(81, NumberUtils.fromBase40("21".toCharArray()));



		assertEquals(81, NumberUtils.fromBase40("0021".toCharArray()));



		assertEquals(3061560805L, NumberUtils.fromBase40("MSTSD5".toCharArray()));

	}



	public void testToBase40() {

		assertEquals("0", new String(NumberUtils.toBase40(0)));

		assertEquals("1", new String(NumberUtils.toBase40(1)));



		assertEquals("10", new String(NumberUtils.toBase40(40)));

		assertEquals("20", new String(NumberUtils.toBase40(80)));



		assertEquals("11", new String(NumberUtils.toBase40(41)));

		assertEquals("21", new String(NumberUtils.toBase40(81)));



		assertEquals("ABCD", new String(NumberUtils.toBase40(1117580)));

		assertEquals("WWWW", new String(NumberUtils.toBase40(2559999)));



		assertEquals("MSTSD5", new String(NumberUtils.toBase40(3061560805L)));



	}



	public void testToBase90() throws Exception {

		assertEquals("#Q&Y.`QY#", new String(NumberUtils.toBase90(10908158098650842L)));

	}



	public void testFromBase90() throws Exception {

		assertEquals(10908158098650842L, NumberUtils.fromBase90("#Q&Y.`QY#".toCharArray()));

	}

Calculate 2 to the power of N

dubby (Dave) — Thu, 06 Sep 2007 10:21:15 GMT

// Calculate 2 to the power of N
// Returns "double" value







        private double TwoPowerN(int n)

        {

            int k;

            int i = 1;



            for (k=1; k <= n; k++)

            {

                    i = i * 2;   

            }

            return i

        }