root/core/gcd.hpp

/* [<][>][^][v][top][bottom][index][help] */

INCLUDED FROM


DEFINITIONS

This source file includes following definitions.
  1. gcd
  2. gcd

   1 #ifndef LIBJMMCG_CORE_GCD_HPP
   2 #define LIBJMMCG_CORE_GCD_HPP
   3 
   4 /******************************************************************************
   5 ** $Header: svn+ssh://jmmcg@svn.code.sf.net/p/libjmmcg/code/trunk/libjmmcg/core/hp_timer.hpp 1970 2016-11-13 21:48:21Z jmmcg $
   6 **
   7 ** Copyright (c) 1993 by J.M.McGuiness, coder@hussar.me.uk
   8 **
   9 ** This library is free software; you can redistribute it and/or
  10 ** modify it under the terms of the GNU Lesser General Public
  11 ** License as published by the Free Software Foundation; either
  12 ** version 2.1 of the License, or (at your option) any later version.
  13 **
  14 ** This library is distributed in the hope that it will be useful,
  15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 ** Lesser General Public License for more details.
  18 **
  19 ** You should have received a copy of the GNU Lesser General Public
  20 ** License along with this library; if not, write to the Free Software
  21 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  22 */
  23 
  24 namespace jmmcg {
  25 
  26 namespace euclid {
  27 
  28         /// Euclid's Greatest Common Divisor Algorithm.
  29         /**
  30          * This is best used for the general case of finding the gcd.
  31          * The Binary method may be faster than Euclid's Algorithm.
  32          * 
  33          * A one-liner for fun! (But it's slower, as it's got an extra jump...)
  34          *              for (u=x,v=y;v;r=u%v,u=v,v=r);
  35          *
  36          * Thanks to Knuth, et al for Euclid's Algorithm and the Binary GCD algorithm.
  37          * See his book "Seminumerical Algorithms. Vol 2", pg 320.
  38          * 
  39          * \param x A value greater that zero.
  40          * \param y A value greater that zero.
  41          * 
  42          * \see binary::gcd()
  43          */
  44         template<class V> constexpr inline V
  45         gcd(const V x, const V y) noexcept(true) __attribute__((const));
  46         template<class V> constexpr inline V
  47         gcd(const V x, const V y) noexcept(true) {
  48                 assert(x>V{} && y>V{});
  49                 V u=x, v=y;
  50                 do {
  51                         const V r=u%v;
  52                         u=v;
  53                         v=r;
  54                 } while (v);
  55                 return u;
  56         }
  57 
  58 }
  59 
  60 namespace binary {
  61 
  62         /// Binary Greatest Common Divisor Algorithm.
  63         /**
  64          * This method is best suited to integers.
  65          * 
  66          * A one-liner for fun! (But it's slower, as it's got an extra jump...)
  67          *              for (u=x,v=y;v;r=u%v,u=v,v=r);
  68          *
  69          * Thanks to Knuth, et al for Euclid's Algorithm and the Binary GCD algorithm.
  70          * See his book "Seminumerical Algorithms. Vol 2", pg 320.
  71          * 
  72          * \param x A value greater that zero.
  73          * \param y A value greater that zero.
  74          * 
  75          * \see euclid::gcd()
  76          */
  77         template<class V> constexpr inline V
  78         gcd(const V x, const V y) noexcept(true) __attribute__((const));
  79         template<class V> constexpr inline V
  80         gcd(const V x, const V y) noexcept(true) {
  81                 assert(x>V{} && y>V{});
  82                 V k=0, u=x, v=y;
  83                 for (; !(u&1) && !(v&1); ++k, (u>>=1), (v>>=1));
  84                 V t=u&1 ? -v : u;
  85                 do {
  86                         do {
  87                                 t>>=1;
  88                         } while (!(t&1));
  89                         t>0 ? (u=t) : (v=-t);
  90                 } while ((t=u-v));
  91                 return u<<k;
  92         }
  93 
  94 }
  95 
  96 }
  97 
  98 #endif

/* [<][>][^][v][top][bottom][index][help] */