Never been to DZone Snippets before?

Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

« Newer Snippets
Older Snippets »
Showing 1-10 of 63 total  RSS 

Fast factorial approximations

Choose and factorial algorithms
   1  
   2    # choose and factorial were stolen from: 
   3    # http://bluebones.net/2007/09/combinatorics-in-ruby/
   4    def choose(m)
   5      self.factorial / (m.factorial * (self - m).factorial)
   6    end
   7    def factorial
   8      (2..self).inject(1) { |f, n| f * n }
   9    end
  10  
  11    def fastchoose(m)
  12      self.gosper / (m.to_f.gosper * (self - m).gosper)
  13    end
  14  
  15    # quick, accurate approximation to factorial
  16    # gleaned from http://mathworld.wolfram.com/StirlingsApproximation.html
  17    def gosper
  18      Math.sqrt( ((2*self) + (1.0/3))*Math::PI ) * ((self/Math::E)**self).round
  19    end
  20  
  21    def stirling
  22      Math::E ** (self*Math.log(self) - self)
  23    end

Bit field class

// 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.

   1  
   2  <bitfield.h>
   3  // This class is a code snippet, it can be freely used and distributed.
   4  // Author: irfan[dot]hamid[at]gmail[dot]com
   5  //
   6  // Bit fields are a widely used mechanism in computing applications.
   7  // Typical applications include memory allocators and Bloom filters. This class
   8  // provides a generic bit field implementation for an arbitrary number of bits.
   9  //
  10  // Implementation notes:
  11  // The use of C++ allows the clean separation of the data structure from the
  12  // interface. The bit field is implemented as a vector of unsigned int that is
  13  // allocated at object creation.
  14  //
  15  // field: The vector that represents the bit field itself
  16  // bit_count: The total number of bits in the bit field
  17  
  18  #ifndef __BITFIELD_H__
  19  #define __BITFIELD_H__
  20  #include <math.h>
  21  #define SZ_UINT sizeof(unsigned int)
  22  
  23  class Bitfield
  24  {
  25  protected:
  26  	unsigned int *field;
  27  	unsigned int bit_count;
  28  
  29  public:
  30  	Bitfield (int bc);
  31  	~Bitfield ();
  32  
  33  	int set (int bit);
  34  	int get (int bit);
  35  	int reset (int bit);
  36  };
  37  
  38  #endif
  39  </bitfield.h>
  40  
  41  <bitfield.cpp>
  42  #include <stdlib.h>
  43  #include "bitfield.h"
  44  
  45  // The constructor takes an int param which gives the number of bits in this
  46  // bit field.
  47  Bitfield::Bitfield (int bc)
  48  {
  49  	bit_count = bc;
  50  
  51  	// E.g, ceil(257/32*8) = 65, 64 ints fully used, last one partially used
  52  	field = (unsigned int*) malloc(((int) ceil(bc/SZ_UINT*8)));
  53  }
  54  
  55  Bitfield::~Bitfield ()
  56  {
  57  	if (field)
  58  		free(field);
  59  }
  60  
  61  // This function sets the corresponding bit in the field equal to 1.
  62  //
  63  // Returns: 0 on success, -1 on error.
  64  int Bitfield::set (int bit)
  65  {
  66  	// Sanity check
  67  	if (bit >= bit_count || !field)
  68  		return -1;
  69  
  70  	// The correct index into the vector will be given by bit/SZ_UNIT. The
  71  	// index into the correct vector element is given by bit%SZ_UNIT, to
  72  	// achieve this, simply left shift 0x01 the appropriate number of times.
  73  	field[bit/SZ_UINT] |= (0x00000001 << (bit%(SZ_UINT*8)-1));
  74  	return 0;
  75  }
  76  
  77  // This function sets the corresponding bit in the field equal to 0.
  78  //
  79  // Returns: 0 on success, -1 on error.
  80  int Bitfield::reset (int bit)
  81  {
  82  	if (bit >= bit_count || !field)
  83  		return -1;
  84  	field[bit/SZ_UINT] &= ~(0x00000001 << (bit%(SZ_UINT*8)-1));
  85  	return 0;
  86  }
  87  
  88  // This function returns the value of the corresponding bit.
  89  //
  90  // Returns: The value of the bit, if the field is initialized and bit index is
  91  // within bounds, -1 otherwise.
  92  int Bitfield::get (int bit)
  93  {
  94  	if (bit >= bit_count || !field)
  95  		return -1;
  96  	return (field[bit/SZ_UINT] & (0x00000001 << (bit%(SZ_UINT*8)-1)) ? 1 : 0);
  97  }
  98  </bitfield.cpp>

Calculate Pi

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

   1  
   2  import random, math
   3  
   4  class calcPi:
   5      def __init__(self):
   6          self.times = pow(10,6)
   7          self.i = 0
   8          self.isnot = 0
   9  
  10      def IsOnCircle(self,x,y):
  11          if math.sqrt(x**2+y**2) < 1:
  12              return True
  13          else:
  14              return False
  15      
  16      def run(self):
  17          for x in range(self.times):
  18              x,y = random.random(),random.random()
  19              if self.IsOnCircle(x,y):
  20                  self.i+=1
  21              else:
  22                  self.isnot+=1
  23      
  24      def getResults(self):
  25          return (float(self.i), float(self.isnot))
  26  
  27      def getPi(self):
  28          self.run()
  29          r = self.getResults()
  30          return r[0]/(r[0]+r[1])*4
  31  
  32  if __name__ == '__main__':
  33      print calcPi().getPi()

First n primes in Ruby

   1  
   2  
   3  #!/usr/local/bin/ruby -w
   4  
   5  require 'benchmark' 
   6  
   7  class Integer    
   8     def prime?         # cf. http://snippets.dzone.com/posts/show/4636
   9       n = self.abs()
  10       return true if n == 2
  11       return false if n == 1 || n & 1 == 0
  12       d = n-1
  13       d >>= 1 while d & 1 == 0
  14       20.times do                               # 20 = k from above
  15         a = rand(n-2) + 1
  16         t = d
  17         y = ModMath.pow(a,t,n)                  # implemented below
  18         while t != n-1 && y != 1 && y != n-1
  19           y = (y * y) % n
  20           t <<= 1
  21         end
  22         return false if y != n-1 && t & 1 == 0
  23       end
  24       return true
  25     end
  26  end
  27   
  28  module ModMath
  29     def ModMath.pow(base, power, mod)
  30       result = 1
  31       while power > 0
  32         result = (result * base) % mod if power & 1 == 1
  33         base = (base * base) % mod
  34         power >>= 1;
  35       end
  36       result
  37     end
  38  end
  39  
  40  
  41  
  42  class Integer
  43  
  44     def primes   # cf. http://snippets.dzone.com/posts/show/3734
  45  
  46        sieve = []
  47        3.step(self, 2) { |i| sieve[i] = i }
  48        sieve[1] = nil
  49        sieve[2] = 2
  50  
  51        3.step(Math.sqrt(self).floor, 2) do |i| 
  52           next unless sieve[i]
  53           (i*i).step(self, i) do |j|
  54              sieve[j] = nil
  55           end
  56        end
  57  
  58        sieve.compact!
  59  
  60     end 
  61  
  62  
  63     def primes2       # cf. http://snippets.dzone.com/posts/show/3734
  64  
  65        primes = [nil, nil].concat((2..self).to_a)
  66  
  67        (2 .. Math.sqrt(self)).each do |i|
  68           next unless primes[i]
  69              (i*i).step(self, i) do |j|
  70                 primes[j] = nil
  71              end
  72        end
  73  	
  74        primes.compact!
  75  
  76     end
  77  
  78  end
  79  
  80  
  81  
  82  class Integer
  83  
  84     @primes_cache ||= 100.primes
  85     #@primes_cache ||= [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
  86     class << self; attr_accessor :primes_cache; end
  87  
  88     @nthprime_limit ||= 5_000_000
  89     class << self; attr_reader :nthprime_limit; end
  90     #class << self; attr_accessor :nthprime_limit; end
  91  
  92  
  93     def sieve_size
  94        num = self.to_i.abs
  95        logn = Math.log(num)
  96        if num < 15985     # cf. http://primes.utm.edu/howmany.shtml
  97           ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
  98        else
  99           ( num * (logn + Math.log(logn) - 0.9427) ).floor
 100        end
 101     end
 102  
 103  
 104     def nthprime
 105  
 106        num = self.to_i.abs
 107  
 108        # set a limit; cf. http://primes.utm.edu/nthprime/: The 5,000,000th prime is 86,028,121.
 109        raise "#{num}: number too big for Integer#nthprime" if num > Integer.nthprime_limit      
 110  
 111        primes_cache_size = Integer.primes_cache.size
 112  
 113        if num > primes_cache_size
 114  
 115           # optional; cf. Integer#nthprime_add and Integer#nthprime_add_mr below
 116           if num >= 50_000 && num < 100_000 && (num - primes_cache_size) < 5000 then return num.nthprime_add end
 117           if num >= 100_000 && num < 200_000 && (num - primes_cache_size) < 15000 then return num.nthprime_add end
 118           if num >= 200_000 && (num - primes_cache_size) < 35000 then return num.nthprime_add end
 119           
 120           logn = Math.log(num)
 121  
 122           if num < 15985       # cf. http://primes.utm.edu/howmany.shtml
 123              limit = ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor
 124              Integer.primes_cache = limit.primes
 125           else
 126              limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
 127              Integer.primes_cache = limit.primes
 128           end
 129  
 130  =begin
 131           elsif num >= 15985 && num <= 1_000_000
 132              limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor
 133              Integer.primes_cache = limit.primes
 134           elsif num > 1_000_000
 135              Integer.primes_cache = 20_000_000.primes
 136           end
 137  =end
 138  
 139        else
 140           return Integer.primes_cache.at(num-1)
 141        end
 142  
 143  
 144        if num <= Integer.primes_cache.size
 145           return Integer.primes_cache.at(num-1)
 146        end
 147  
 148        if Integer.primes_cache.size > 2_500_000
 149           num.nthprime_add_mr 
 150        else
 151           num.nthprime_add   # faster for a (relatively) small prime cache
 152        end
 153  
 154     end
 155  
 156  
 157     def nthprime_add       #  add primes to Integer.primes_cache up to the nth prime
 158  
 159        num = self.to_i.abs
 160  
 161        if num <= Integer.primes_cache.size
 162           return Integer.primes_cache.at(num-1)
 163        end
 164  
 165        last_prime = Integer.primes_cache.last
 166        last_prime_divmod = last_prime.divmod(6)
 167  
 168        if last_prime_divmod.last == 1
 169           i = last_prime_divmod.first
 170           Integer.primes_cache.pop      # avoid a duplicate prime
 171        else
 172           i = last_prime_divmod.first + 1
 173        end
 174  
 175  
 176        while Integer.primes_cache.size < num
 177  
 178           n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
 179           n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369
 180           i += 1
 181  
 182           [n1, n2].each do |p| 
 183  
 184              next if p % 5 == 0 || p % 7 == 0
 185  
 186              next_num_sqrt = Math.sqrt(p).floor
 187  
 188              Integer.primes_cache.each do |d| 
 189                 if d > next_num_sqrt 
 190                    Integer.primes_cache << p
 191                    #print "\r\e[0K#{Integer.primes_cache.size}"
 192                    #$stdout.flush
 193                    break
 194                 elsif p % d == 0
 195                    break
 196                 end  
 197              end
 198           end   # each  
 199  
 200        end  # while
 201  
 202        return Integer.primes_cache.at(num-1)
 203  
 204     end
 205  
 206  
 207     def nthprime_add_mr       #  add Miller-Rabin primes to Integer.primes_cache up to the nth prime
 208  
 209        num = self.to_i.abs
 210  
 211        if num <= Integer.primes_cache.size
 212           return Integer.primes_cache.at(num-1)
 213        end
 214  
 215        last_prime = Integer.primes_cache.last
 216        last_prime_divmod = last_prime.divmod(6)
 217  
 218        if last_prime_divmod.last == 1
 219           i = last_prime_divmod.first
 220           Integer.primes_cache.pop
 221        else
 222           i = last_prime_divmod.first + 1
 223        end
 224  
 225  
 226        while Integer.primes_cache.size < num
 227  
 228           search_next_prime = true
 229  
 230           while search_next_prime
 231  
 232              n1 = 6*i+1       
 233              n2 = n1+4        
 234              i += 1
 235  
 236              [n1, n2].each do |p| 
 237  
 238                 next if p % 5 == 0 || p % 7 == 0
 239  
 240                 if p.prime?
 241                    Integer.primes_cache << p
 242                    search_next_prime = false
 243                    #print "\r\e[0K#{Integer.primes_cache.size}"
 244                    #$stdout.flush
 245                 end
 246              end
 247           end
 248  
 249        end  #  first while loop
 250  
 251        return Integer.primes_cache.at(num-1)
 252  
 253     end
 254  
 255  
 256  #------------------------------------------- some additional prime methods
 257  
 258  
 259     def next_primes_in_cache    # next primes in Integer.primes_cache
 260  
 261        search_next_primes = self.to_i.abs
 262  
 263        last_prime = Integer.primes_cache.last
 264        last_prime_divmod = last_prime.divmod(6)
 265  
 266        if last_prime_divmod.last == 1
 267           i = last_prime_divmod.first
 268           Integer.primes_cache.pop
 269        else
 270           i = last_prime_divmod.first + 1
 271        end
 272  
 273        while search_next_primes > 0
 274  
 275           n1 = 6*i+1       
 276           n2 = n1+4        
 277           i += 1
 278  
 279           [n1, n2].each do |p| 
 280  
 281              next if p % 5 == 0 || p % 7 == 0
 282              next_num_sqrt = Math.sqrt(p).floor
 283  
 284              Integer.primes_cache.each do |d| 
 285                 if d > next_num_sqrt 
 286                    Integer.primes_cache << p
 287                    search_next_primes -= 1
 288                    #print "\r\e[0K#{Integer.primes_cache.size}"
 289                    #$stdout.flush
 290                    break
 291                 elsif p % d == 0
 292                    break
 293                 end  
 294              end
 295           end 
 296        end   # while
 297  
 298        Integer.primes_cache.last(self.to_i.abs)
 299  
 300     end
 301  
 302  
 303     def next_mr_prime     # next Miller-Rabin prime
 304       
 305        last_num = self.to_i.abs
 306        last_num_divmod = last_num.divmod(6)
 307  
 308        if last_num_divmod.last == 1
 309           i = last_num_divmod.first
 310        else
 311           i = last_num_divmod.first + 1
 312        end
 313  
 314        next_prime = nil
 315        search_next_prime = true
 316  
 317        while search_next_prime
 318  
 319           n1 = 6*i+1       
 320           n2 = n1+4        
 321           i += 1
 322  
 323           [n1, n2].each do |p| 
 324  
 325              next if p % 5 == 0 || p % 7 == 0
 326  
 327              if p > last_num && p.prime?
 328                 next_prime = p
 329                 search_next_prime = false 
 330                 break
 331              end
 332           end
 333        end
 334  
 335           next_prime
 336  
 337     end
 338  
 339  
 340  
 341     def next_mr_primes(n)     # next Miller-Rabin primes
 342  
 343        last_num = self.to_i.abs
 344        last_num_divmod = last_num.divmod(6)
 345  
 346        if last_num_divmod.last == 1
 347           i = last_num_divmod.first
 348        else
 349           i = last_num_divmod.first + 1
 350        end
 351  
 352        next_primes = []
 353        search_next_primes = n.to_i.abs
 354  
 355        while search_next_primes > 0
 356  
 357           n1 = 6*i+1      
 358           n2 = n1+4       
 359           i += 1
 360  
 361           [n1, n2].each do |p| 
 362  
 363              next if p % 5 == 0 || p % 7 == 0
 364  
 365              if p > last_num && p.prime?
 366                 next_primes << p
 367                 search_next_primes -= 1 
 368              end
 369           end
 370        end
 371  
 372        next_primes
 373  
 374     end
 375  
 376  end
 377  
 378  
 379  
 380  num1 = 10_000
 381  num1 = 1_000
 382  num1 = 5_000
 383  
 384  num2 = 210_349
 385  num2 = 100_125
 386  num2 = 55_000
 387  
 388  ret1 = nil
 389  ret2 = nil
 390  ret3 = nil
 391  ret4 = nil
 392  
 393  
 394  Benchmark.bm(16) do |t| 
 395  
 396     t.report("first #{num1} primes: ") do
 397        ret1 = num1.nthprime_add
 398     end 
 399  
 400     Integer.primes_cache.clear
 401     Integer.primes_cache.concat(100.primes)
 402  
 403     t.report("first #{num1} primes: ") do
 404        ret2 = num1.nthprime_add_mr
 405     end 
 406  
 407     Integer.primes_cache.clear
 408     Integer.primes_cache.concat(100.primes)
 409  
 410     t.report("first #{num1} primes: ") do
 411        ret3 = num1.nthprime
 412     end 
 413  
 414     Integer.primes_cache.clear
 415     Integer.primes_cache.concat(100.primes)
 416  
 417     t.report("first #{num2} primes: ") do
 418        ret4 = num2.nthprime
 419     end 
 420  
 421  end
 422  
 423  
 424  puts
 425  puts "the #{num1}th prime number: #{ret1}"
 426  puts "the #{num1}th prime number: #{ret2}"
 427  puts "the #{num1}th prime number: #{ret3}"
 428  puts "the #{num2}th prime number: #{ret4}"
 429  puts
 430  puts "the #{num1}th prime: #{Integer.primes_cache.at(num1-1)}"
 431  puts "the #{num2}th prime: #{Integer.primes_cache.at(num2-1)}"
 432  puts
 433  p Integer.primes_cache.first(10)
 434  p Integer.primes_cache.last(10)
 435  p Integer.primes_cache.size
 436  puts
 437  
 438  
 439  p 15.next_primes_in_cache
 440  
 441  puts 594_213.next_mr_prime
 442  
 443  p 149_137.next_mr_primes(5)
 444  
 445  puts
 446  p 1_000.sieve_size
 447  p 1_000_000.sieve_size
 448  p 2_500_000.sieve_size
 449  p 5_000_000.sieve_size
 450  p 10_000_000.sieve_size
 451  p 100_000_000.sieve_size
 452  
 453  
 454  =begin
 455  
 456  # check with primegen.c, http://cr.yp.to/primegen.html
 457  primes 2 48611 | nl | tail -n 1
 458  primes 2 104729 | nl | tail -n 1
 459  primes 2 679277 | nl | tail -n 1
 460  primes 2 1301423 | nl | tail -n 1
 461  primes 2 2903603 | nl | tail -n 1
 462  primes 1303129 1303283 | nl
 463  primes 594213 $((594213 + 500)) | head -n 1
 464  primes 149137 $((149137 + 1000)) | head -n 5
 465  
 466  # further prime tools from http://libtom.org
 467  - LibTomMath
 468     bn_mp_prime_is_prime.c
 469     bn_mp_prime_miller_rabin.c
 470  - TomsFastMath
 471     fp_isprime.c
 472     fp_prime_miller_rabin.c
 473  
 474  =end

ruby is_numeric? extension to String

   1  
   2  class String
   3    def is_numeric?
   4      Float self rescue false
   5    end
   6  end

Pascal's triangle in Ruby

   1  
   2  
   3  # cf. http://www.ruby-forum.com/topic/97105
   4  
   5  6.times { k=0; p $*.map!{|i|k+k=i} << 1 }
   6  
   7  # ... or ...
   8  
   9  ar=[]; 6.times { k=0; p ar.map!{|i|k+k=i} << 1 }
  10  
  11  
  12  # a more general approach (for polynomials)
  13  
  14  class Polynomial
  15  
  16     def triangle(nterms, row, pos=nil)
  17  
  18        return nil if nterms < 2 || row < 1
  19        nterms = nterms - 2
  20        num_of_rows = row
  21  
  22        var1 = 0 + nterms   
  23        var2 = 1 + nterms
  24        var3 = 3 + nterms 
  25    
  26        ar1 = [0, 1, 0]   # first row
  27        var1.times { ar1.push(0) }
  28        var1.times { ar1.unshift(0) }
  29  
  30        ar2 = []
  31        ar3 = []
  32        ar4 = [[1]]
  33  
  34        for num in 0..(num_of_rows - 1)  
  35  
  36           nextnum = ar1.size - var2
  37  
  38           for nextn in 1..nextnum
  39              sum = 0
  40              count = 0
  41              ar1.each do |n|  
  42                 count += 1 
  43                 if count < var3 then t = sum += n; ar2 << t else break end 
  44              end
  45  
  46              ar3 << ar2.last
  47              ar2.clear
  48              ar1.shift
  49  
  50           end   # second for-loop
  51  
  52           ar1.clear
  53           ar1 << ar3
  54           ar1.flatten!
  55  
  56           var2.times { ar1.push(0) }
  57           var2.times { ar1.unshift(0) }
  58  
  59           ar4 << ar3
  60           ar3 = []
  61  
  62        end  # first for-loop
  63  
  64        if !pos.nil?
  65           ret = ar4.at(row).at(pos)
  66           return "No such position: #{pos} in row: #{row}" if ret.nil?
  67           ret
  68        else
  69           ar4.map! { |r| r.join('-') }
  70           ar4
  71        end
  72     end 
  73  end 
  74  
  75  
  76  puts Polynomial.new.triangle(2, 5)
  77  puts Polynomial.new.triangle(3, 5)
  78  puts Polynomial.new.triangle(4, 5)
  79  puts Polynomial.new.triangle(5, 5)
  80  puts Polynomial.new.triangle(5, 4, 8)
  81  puts Polynomial.new.triangle(4, 9)
  82  puts Polynomial.new.triangle(4, 9, 10)
  83  
  84  
  85  #------------------------
  86  
  87  
  88  class Integer
  89     def fak
  90        f=1
  91        (2..self).each { |i| f *= i   }
  92        f
  93     end
  94  end
  95  
  96  module Enumerable
  97     def sum
  98        inject { |n, m| n + m  }
  99     end
 100  end
 101  
 102  # cf. http:/