数学問題-二重根号の簡約（R 版）

数学問題-3 で，「二重根号の簡約は，直接的に解を求める事が出来ない（simplify では解が得られない）」と書いたが，解を求めるプログラムを R で書いてみた。エラーチェックや，数式を引数にするという書式にしたので，ちょっと長くなった。

使用例

> double.root("sqrt(6-sqrt(32))")
$equation
[1] "2 - sqrt(2)"

$expression
expression(2 - sqrt(2))

$eval
[1] 0.5857864

> double.root("sqrt(6-2*sqrt(8))")
$equation
[1] "2 - sqrt(2)"

$expression
expression(2 - sqrt(2))

$eval
[1] 0.5857864

> double.root("sqrt(106+2*sqrt(1288))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+sqrt(5152))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+4*sqrt(322))")
$equation
[1] "sqrt(92) + sqrt(14)"

$expression
expression(sqrt(92) + sqrt(14))

$eval
[1] 13.33332

> double.root("sqrt(106+sqrt(368))")
[1] "can't simplify"
> double.root("sqrt(-16+sqrt(368))")
[1] "can't simplify"
> double.root("sqrt(16-sqrt(-368))")
[1] "square root of negative value is not allowed (inner sqrt)"
> double.root("sqrt(16+sqrt(-368))")
[1] "square root of negative value is not allowed (inner sqrt)"
> double.root("sqrt(16-sqrt(368))")
[1] "square root of negative value is not allowed (outer sqrt)"
> double.root("sqrt(18+2*sqrt(81))")
$equation
[1] "6"

$expression
expression(6)

$eval
[1] 6

> double.root("sqrt(13+sqrt(144))")
$equation
[1] "5"

$expression
expression(5)

$eval
[1] 5

プログラム

double.root = function(equation) {
  solve = function(a.plus.b, a.mult.b) {
    for (b in 1:sqrt(a.mult.b)) {
      if (a.mult.b %% b == 0 && a.plus.b == b + a.mult.b %/% b) {
        return(list(a = a.mult.b %/% b, b = b))
      }
    }
    return(NA)
  }
  simplify = function(a) {
    int.root = trunc(sqrt(a))
    if (int.root ^ 2 == a) {
      return(sprintf("%d", int.root))
    } else {
      return(sprintf("sqrt(%d)", a))
    }
  }
  inner = sub("\\)$", "", sub("sqrt\\(", "", equation))
  if (grepl("sqrt\\(-[0-9]*\\)", equation)) {
    return("square root of negative value is not allowed (inner sqrt)")
  } else if (eval(parse(text = inner)) < 0) {
    return("square root of negative value is not allowed (outer sqrt)")
  }
  result1 = eval(parse(text = equation))
  sign = ifelse(grepl("-", equation), "-", "+")
  equation = gsub("sqrt\\(", "", equation)
  equation = gsub("\\)", "", equation)
  # print(equation)
  parsed = as.integer(unlist(strsplit(equation, "[-+*]")))
  # print(parsed)
  if (is.na(parsed[1])) {
    return("can't simplify")
  }
  if (length(parsed) == 3) {
    parsed[2] = parsed[2] ^ 2 * parsed[3]
  }
  parsed[2] = parsed[2] / 4
  if (parsed[2] != trunc(parsed[2])) {
    return("can't simplify")
  }
  # cat(parsed[1], parsed[2], "\n")
  res = solve(parsed[1], parsed[2])
  if (length(res) != 2) {
    return("can't simplify")
  }
  res.str = paste(simplify(res$a), sign, simplify(res$b))
  result2 = eval(parse(text = res.str))
  if (all.equal(result1, result2)) {
      if (sign == "+") {
        res.str = as.character(sqrt(res$a) + sqrt(res$b))
      } else {
        res.str = as.character(sqrt(res$a) - sqrt(res$b))
      }
    result = parse(text = res.str)
    return(list(
      equation = res.str,
      expression = result,
      eval = eval(result)
    ))
  }
  return(sprintf("error: original = %g, simplified = %g", result1, result2))
}