不循环求GCD - R

Question

所以我正在尝试学习 R 并使用许多资源，包括一本名为“使用 R 发现统计数据”的书和一堆其他很酷的电子书。

我知道编程中一个很好的方法是欧几里得算法。

循环实现可以这样实现：

 gcd(x,y) //assuming x is the largest value
 //do
 r = x%y;
 x = y;
 y = r;
 //while r != 0;
 return x;

在 Google、SO 和 Youtube 上多次搜索刷新了我对 gcd 算法的记忆后，我找不到一个不使用循环的算法。甚至递归方法似乎也使用循环。

在 R 中如何在不使用循环或 if 语句的情况下实现这一点？

提前致谢。

Answer 1

从字面上使用语句“没有循环或

if

语句”，这是一个使用

ifelse

的递归版本：

gcd <- function(x,y) {
  r <- x%%y;
  return(ifelse(r, gcd(y, r), y))
}

人们可能没有想到，但这实际上是矢量化的：

gcd(c(1000, 10), c(15, 10))
[1]  5 10

使用

if

的解决方案无法处理长度大于 1 的向量。

Answer 2

减少两个整数的 GCD 使您能够计算任何整数序列（排序或未排序）的 GCD：

gcd2 <- function(a, b) {
  if (b == 0) a else Recall(b, a %% b)
}

gcd <- function(...) Reduce(gcd2, c(...))

Answer 3

你可以递归地解决它。

euclids <- function(x,y){
        theMax = max(x,y)
        theMin = min(x,y)

        if (theMax == theMin) return (theMax)
        else return (euclids(theMin, theMax-theMin))
}

Answer 4

我不确定是否有“无循环”版本，但这里有一个使用 cpp11 包的“无R循环”版本。

C++ 函数

library(cpp11)

cpp_source(code = '
#include <cpp11.hpp>
#include <Rinternals.h>

// #include <numeric>

#define R_NO_REMAP

[[cpp11::register]]
double cpp_gcd2(double x, double y, double tol = 0, bool na_rm = true){
  double zero = 0.0;
  if (!na_rm && ( !(x == x) || !(y == y) )){
      return NA_REAL;
  }
  // GCD(0,0)=0
  if (x == zero && y == zero){
      return zero;
  }
  // GCD(a,0)=a
  if (x == zero){
      return y;
  }
  // GCD(a,0)=a
  if (y == zero){
      return x;
  }
  double r;
  // Taken from number theory lecture notes
  while(std::fabs(y) > tol){
      r = std::fmod(x, y);
      x = y;
      y = r;
  }
  return x;
}

[[cpp11::register]]
SEXP cpp_gcd2_vectorised(SEXP x, SEXP y, double tol, bool na_rm){
  if (tol < 0 || tol >= 1){
      Rf_error("tol must be >= 0 and < 1");
  }
  int xn = Rf_length(x);
  int yn = Rf_length(y);
  int n = std::max(xn, yn);
  if (xn == 0 || yn == 0){
      n = 0;
  }
  Rf_protect(x = Rf_coerceVector(x, REALSXP));
  Rf_protect(y = Rf_coerceVector(y, REALSXP));
  SEXP out = Rf_protect(Rf_allocVector(REALSXP, n));
  double *p_out = REAL(out);
  double *p_x = REAL(x);
  double *p_y = REAL(y);
  int xi;
  int yi;
  for (int i = 0; i < n; ++i) {
      xi = i % xn;
      yi = i % yn;
      p_out[i] = cpp_gcd2(p_x[xi], p_y[yi], 0, true);
  }
  Rf_unprotect(3);
  return out;
}

[[cpp11::register]]
SEXP cpp_gcd(SEXP x, double tol, bool na_rm, int start, bool round){
  if (tol < 0 || tol >= 1){
      Rf_error("tol must be >= 0 and < 1");
  }
  int n = Rf_length(x);
  double *p_x = REAL(x);
  SEXP out = Rf_protect(Rf_allocVector(REALSXP, std::min(n, 1)));
  double *p_out = REAL(out);
  double gcd = p_x[start - 1];
  for (int i = start; i < n; ++i) {
      gcd = cpp_gcd2(gcd, p_x[i], tol, na_rm);
  }
  if (round && tol > 0){
      double factor = std::pow(10, std::ceil(std::fabs(std::log10(tol))) + 1);
      gcd = std::round(gcd * factor) / factor;
  }
  p_out[0] = gcd;
  Rf_unprotect(1);
  return out;
}')

基本用法

# Let's make helpers..

# GCD of 2 single numbers
gcd2 <- function(x, y, tol = sqrt(.Machine$double.eps), na_rm = TRUE){
  cpp_gcd2(x, y, tol = tol, na_rm = na_rm)
}

# GCD of a vector of numbers
gcd <- function(x, tol = sqrt(.Machine$double.eps), na_rm = TRUE, round = FALSE){
  cpp_gcd(x, as.double(tol), na_rm, start = 1L, round = round)
}

# GCD vectorised across 2 vectors of numbers
gcd2_vectorised <- function(x, y, tol = sqrt(.Machine$double.eps), na_rm = TRUE){
  cpp_gcd2_vectorised(x, y, tol = tol, na_rm = na_rm)
}

library(bench)

# Binary version
gcd2(25, 125)
#> [1] 25

# GCD of all values in a vector
gcd(c(5, 25, 125))
#> [1] 5

x <- seq(0, len = 10^6, by = 5)

gcd(x)
#> [1] 5
mark(gcd(x))
#> # A tibble: 1 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 gcd(x)       17.2ms     19ms      49.8        0B        0

# Vectorised across two vectors

all(gcd2_vectorised(x, 5) == 5)
#> [1] TRUE

^{创建于 2023-11-17，使用 reprex v2.0.2}

Answer 5

通过几次模运算很容易完成。可悲的是，我将我的个人

gcd

代码留在了另一台机器上（在遥远的星系中） - 但您可以在

numbers

或

pracma

包中找到源代码。

顺便说一句，这是查找现有代码的好方法：

library(sos); ???gcd