/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/math_prime_numbers.h

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_MATH_PRIME_NUMBERS_H
#define OPENTISSUE_CORE_MATH_MATH_PRIME_NUMBERS_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/math/math_constants.h>
#include <OpenTissue/core/math/math_random.h>

#include <cmath>   //--- for floor and sqrt


namespace OpenTissue
{

  namespace math
  {
    
        
    namespace detail
    {
      
      inline bool withness(int a,int n)
      {
        int ndec = n-1;int d  = 1; int k = ndec;
        if(k!=0)
        {
          int c=0;
          while((k&0x80000000)==0)
            k<<=1;c++;
          while(c<32)
          {
            int x =d;
            d=(d*d)%n;
            if((d==1)&&(x!=1)&&(x!=ndec))
              return true;//--- Notrival square root of 1 was discovered.
            if((k&0x80000000)!=0)
              d=(d*a)%n;
            k<<=1;c++;
          }
        }
        if(d!=1)
          return true;
        return false;
      }

    }//namespace detail


    inline bool trial_division(int n)
    {
      using std::sqrt;
      using std::floor;

      double sqrN = sqrt(static_cast<double>(n));
      int j = static_cast<int>( floor(sqrN) );
      for(int i=2;i<=j;++i)
      {
        if((n%i)==0)
          return false;
      }
      return true;
    }

    inline int modular_exponentiation(int a,int b,int n)
    {
      int d = 1;
      while(b!=0)
      {
        if((b&1)!=0)
          d=(d*a)%n;
        a = (a*a)%n;
        b>>=1;
      }
      return d;
    }

    inline bool pseudo_prime(int n)
    {
      if((n==1)||(n==2))//--- trivial cases
        return true;
      if(modular_exponentiation(2,n-1,n)!=1)
        return false;//--- definitely not a prime, i.e. composite
      return true;//--- possible prime or a base-2 pseudoprime
    }


    inline bool miller_rabin(int n,int s)
    {
      using std::floor;

      static Random<double> random;

      int ndec = n-2;
      for(int j=s;j>0;--j)
      {
        int a = static_cast<int>(  floor( random()*ndec )  ) + 1;//--- random(1,n-1)
        if(detail::withness(a,n))
          return false;//--- n was composite.
      }
      return true;//--- n is almost surely prime.
    }

    inline int prime_search(int n)
    {
      int l = n;
      if((l%2)==0)
        l--;
      for(;l>1;l=l-2)
      {
        //      if(trial_division(l))
        //      if(pseudo_prime(l))
        if(miller_rabin(l,4))
          break;
      }
      int o = n;
      if((o%2)==0)
        o++;
      int max_val = detail::highest<int>();
      for(;o<max_val;o=o+2)
      {
        //      if(trial_division(o))
        //      if(pseudo_prime(o))
        if(miller_rabin(o,4))
          break;
      }
      //--- now l and o must contain two closest
      //--- primes to n such that l <= n <= o
      //--- figure out which one that is closest
      //--- to n and return it.
      if((n-l)<(o-n))
        return l;
      return o;
    }

    inline int gcd_euclid_algorithm(int a,int b)
    {
      if(b==0)
        return a;
      else
        return gcd_euclid_algorithm(b,a%b);
    }
    
    inline bool is_relative_prime(int a,int b)
    {
      if(gcd_euclid_algorithm(a,b)==1)
        return true;
      return false;
    }

  } // namespace math

} // namespace OpenTissue

//OPENTISSUE_CORE_MATH_MATH_PRIME_NUMBERS_H
#endif