factoring.hpp

Go to the documentation of this file.

#ifndef LIBJMMCG_CORE_FACTORING_HPP
#define LIBJMMCG_CORE_FACTORING_HPP
 
/******************************************************************************
** Copyright © 1993 by J.M.McGuiness, coder@hussar.me.uk
**
** This library is free software; you can redistribute it and/or
** modify it under the terms of the GNU Lesser General Public
** License as published by the Free Software Foundation; either
** version 2.1 of the License, or (at your option) any later version.
**
** This library is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
** Lesser General Public License for more details.
**
** You should have received a copy of the GNU Lesser General Public
** License along with this library; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
 
#include "gcd.hpp"
#include "blatant_old_msvc_compiler_hacks.hpp"
 
#include <set>
#include <vector>
 
namespace jmmcg { namespace LIBJMMCG_VER_NAMESPACE { namespace factoring {
 
using collection_type=std::vector<unsigned long>;
 
/// Factoring by division.
/**
 * \param   x  The number to be factored. Must be greater than zero.
 * \return  The array that has the factors placed in it.
 * 
 * Thanks to Knuth, et al for this routine. See his book "Seminumerical Algorithms. Vol 2".
 * 
 * \see factoring::monte_carlo()
 */
inline collection_type __fastcall
division(const collection_type::value_type x) noexcept(false) {
   assert(x>collection_type::value_type{});
   std::set<collection_type::value_type> factors;
   collection_type::value_type count=1lu, tst_divs=2lu, n=x;
   while (n!=1) {
      if (n%tst_divs) {
         if (n>1) {
            if (count==1) {
               ++tst_divs;
               ++count;
            } else if (count<3) {
               tst_divs+=2;
               ++count;
            } else {
               count++&1 ? (tst_divs+=2) : (tst_divs+=4);
               if ((count>8) && (!(tst_divs%5) || !(tst_divs%7))) {
                  count++&1 ? (tst_divs+=2) : (tst_divs+=4);
               }
            }
         } else {
            factors.insert(n);
            return collection_type{factors.begin(), factors.end()};
         }
      } else {
         factors.insert(tst_divs);
         n/=tst_divs;
      }
   }
   return collection_type{factors.begin(), factors.end()};
}
 
/// Monte Carlo factorisation.
/**
 * \todo FIXME
 * 
 * \param   y  The number to be factored.
 * \return  The array that has the factors placed in it.
 * 
 * Thanks to Knuth, et al for this routine. See his book "Seminumerical Algorithms. Vol 2".
 * 
 * \see factoring::division()
 */
inline collection_type __fastcall
monte_carlo(const collection_type::value_type y) noexcept(false) {
   collection_type factors;
   collection_type::value_type x=5, x_=2, k=1, l=1, n=y, g;
b2:
   if (n&1) {
      factors.emplace_back(n);
      return factors;
   }
b3:
   if ((g=euclid::gcd((x_-x),n))==1) {
      goto b4;
   } else if ((g=n)!=0) {
      return factors;
   } else {
      n/=g;
      x%=n;
      x_%=n;
      goto b2;
   }
b4:
   if (!--k) {
      x_=x;
      l<<=1;
      k=l;
   }
   x=(x*x+1)%n;
   goto b3;
}
 
} } }
 
#endif