exp(pi * sqrt(n))が整数に近くなるnを探す DouKaku?

1以上200未満の整数nのうち、 exp(pi * sqrt(n))がほとんど整数であるようなnを求めるコードを書いてください。なお、expは底がeである指数関数 - Wikipedia、 piは円周率、sqrtは平方根です。また「ほとんど整数である」とは整数からプラスマイナス0.0001の範囲にあることとします。

Pythonで34行のスクリプトを書いて得られた出力の例が下のようになります。

37 199148647.999978
58 24591257752.000000
67 147197952743.999999
163 262537412640768744.000000

この問題は光成さんに教えて頂いた e^{π*sqrt{163}}≒26253741640768744 が元になっています。ご協力ありがとうございました。

1から199まで(Range)の整数から、選ぶ(Select)、
基準は、
「x=Exp[Pi Sqrt[n]]として(With)、Abs[x-Round@x]が10^-4以下のもの」

「10^-4を0.0001と書いてはいけない」という問題かな

Select[Range@199,
  Function[{n}, 
    With[{x = Exp[Pi Sqrt@n]},
      Abs[x - Round@x] <= 10^-4]]]

もともとかなり古いライブラリだけど。。。使い勝手は悪くないと思う。

CReal モジュールを使う

このモジュールではデフォルトで小数点以下40桁まで表示。

see: CReal

module Main (main) where

import CReal

default (CReal)

f :: CReal -> (Int, CReal)
f n = (ceiling n, exp $ pi * sqrt n)

g :: CReal -> Bool
g n = 0.0001 > if 0.5 > n' - n then n' - n else 1 - (n' - n)
  where n' = fromInteger $ ceiling n

main :: IO ()
main = mapM_  (putStrLn . show) $ filter (g . snd) $ map f $ [1..200]

{-
*Main> :main
(37,199148647.9999780465518567665009238753359004336659)
(58,24591257751.9999998222132414695761923552658122276102)
(67,147197952743.9999986624542245068292613125786285081833)
(163,262537412640768743.9999999999992500725971981856888793538563)
-}

あっ。200未満か。まぁいいか。

GNU拡張使用。
精度の問題か163はひっかからず。

実行結果は
37 199148647.999978
58 24591257752.000000
67 147197952743.999999

#include <stdio.h>
#define __USE_GNU
#include <math.h>

void func(int n)
{
    long double d;
    int n_normal, n_plus, n_minus;

    /* exp(pi * sqrt(n)) を求める */
    d = expl(M_PIl * sqrtl((double)n));

    /* 精度倍する。1を加減算して桁上がりするか */
    n_normal = fmodl(d*10000.0, 10000.0);
    n_plus   = n_normal + 1;
    n_minus  = n_normal - 1;

    if ((n_normal % 10000) == 0
     || (n_plus   % 10000) == 0
     || (n_minus  % 10000) == 0) {
        printf("%d %Lf\n", n, d); 
    }   
}

int main()
{
    int n;
                                                                        
    for (n = 1; n < 200; n++) {                                         
        func(n);
    }

    return 0;
}

>「10^-4を0.0001と書いてはいけない」

そうなんですか？

ええと、long doubleは18桁くらいの精度ですよね。
67の計算には足りて、163の計算には足りないのでしょうね、きっと。

Apfloatを使いました。

それだけじゃ何も芸がない感じだったので、implicitでApfloatをあまり意識せずに書ける様にしました。

see: Apfloat

import org.apfloat._
import scala.runtime._

object findAlmostInt {
  implicit def toApfloat[A <: AnyVal](i:A):Apfloat = {
    new Apfloat(i.toString, 32)
  }
  implicit def toDouble(i:Apfloat):RichDouble = i.doubleValue.toDouble
  val PI = ApfloatMath.pi(32)

  def apply(i:int, j:int) = {
    (i to j-1).filter{n =>
      val x = ApfloatMath.exp(PI.multiply(ApfloatMath.sqrt(n)))
      val y = x.ceil.subtract(x)
      (if(y<0.5){y}else{1.subtract(y)}).abs < 0.0001
    }.toList
  }
}

println(findAlmostInt(1,200))

Decimalを使って多倍長計算。sqrtは組み込みのものを使い、expは安直にマクローリン展開で収束するまで回しています。psycoを使ってますが、それでもかなり遅いです。

see: Pythonで10進数計算したり数値を丸めたりする

from decimal import *
import itertools

def decimal_exp(x):
    curr, bunsi, bunbo = Decimal("1.0"), Decimal("1.0"), 1
    for k in itertools.count(1):
        bunsi *= x
        bunbo *= k
        next = curr + bunsi / bunbo
        if curr == next:
            return curr
        curr = next

def main():
    pi = Decimal("3.141592653589793238462643383279")
    for n in xrange(1, 200):
        x = decimal_exp(pi * Decimal(n).sqrt()) # result
        i = x.quantize(Decimal("1.0"), rounding=ROUND_HALF_UP) # nearest integer
        if (x - i).quantize(Decimal("0.00001"), rounding=ROUND_DOWN) == 0:
            print n, x

if __name__ == '__main__':
    try:
        import psyco
    except ImportError:
        pass
    else:
        psyco.full()
    main()

あああ、すみません。組み終わった後に ".00..." という表記を自分好みの "0.00.." という表記に書き換えたときにエンバグしちゃったみたいで、37 が表示されなくなってました。プラス評価もらったのに申し訳ない。

--- a.orig	Sat Sep 15 15:12:33 2007
+++ a.py	Sat Sep 15 15:12:50 2007
@@ -16,7 +16,7 @@
     for n in xrange(1, 200):
         x = decimal_exp(pi * Decimal(n).sqrt()) # result
         i = x.quantize(Decimal("1.0"), rounding=ROUND_HALF_UP) # nearest integer
-        if (x - i).quantize(Decimal("0.00001"), rounding=ROUND_DOWN) == 0:
+        if (x - i).quantize(Decimal("0.0001"), rounding=ROUND_DOWN) == 0:
             print n, x
 
 if __name__ == '__main__':

f の定義がちょっとアレげなので、リファクタリング

f :: Integer -> (Integer, CReal)
f n = (n, exp $ pi * sqrt (fromInteger n))

素直にbigdecimalを使って。

see: Rubyリファレンスマニュアル - bigdecimal

require "bigdecimal"
require "bigdecimal/math.rb"

include BigMath

def exp_pi_sqrt(x)
  tmp = exp(PI(30) * sqrt(BigDecimal(x.to_s), 30), 30)
  err = BigDecimal 0.0001.to_s
  if (tmp > (tmp.floor - err)) and (tmp < (tmp.floor + err))
    print "#{x} #{tmp.to_s("F")}\n"
  elsif (tmp > (tmp.ceil - err)) and (tmp < (tmp.ceil + err))
    print "#{x} #{tmp.to_s("F")}\n"
  end
end

=begin
for i in 1 .. 200
  exp_pi_sqrt i
end
=end

強引に解いてみました。BigDecimalを使用していますがライブラリがないので自前で計算しています。expの計算方法が非効率的なので非常に時間がかかります（苦笑）。

import java.math.BigDecimal;
import java.math.MathContext;

public class Sample {
    private static final double EPS = 0.000001;
    private static final double EPS2 = 1e-25;
    private static final BigDecimal TWO = new BigDecimal(2);
    private static final BigDecimal PI = new BigDecimal
        ("3.1415926535897932384626433");
    public static BigDecimal exp(BigDecimal x) {
        int i = 1;
        BigDecimal b = new BigDecimal(i++);
        BigDecimal c = x;
        BigDecimal a = new BigDecimal(1.0);
        BigDecimal d;
        while ((d = c.divide(b, MathContext.DECIMAL128)).doubleValue() > EPS) {
            a = a.add(d);
            b = b.multiply(new BigDecimal(i++));
            c = c.multiply(x);
        }
        return a;
    }

    public static BigDecimal sqrt(BigDecimal a) {
        BigDecimal x = new BigDecimal(10);
        BigDecimal delta;
        do {
            delta = (x.multiply(x).subtract(a)).divide(TWO.multiply(x), 
                                                       MathContext.DECIMAL128);
            x = x.subtract(delta);
        } while (delta.abs().doubleValue() > EPS2);
        return x;
    }

    public static void main(String[] args) {
        for (int i = 1; i < 200; i++) {
            BigDecimal er = exp(PI.multiply(sqrt(new BigDecimal(i))));
            BigDecimal rer = new BigDecimal(er.toBigInteger());
            double a = er.subtract(rer).doubleValue();
            if (a < 0.0001 || a > 0.9999) {
                System.out.printf("%d: %f%n", i, er);
            }
        }
    }
}

ふつうに Math::BigFloat を使ってます。

実行結果：
37 => 199148647.99997804655185676650092387533591914368128
58 => 24591257751.999999822213241469576192355288556885837
67 => 147197952743.99999866245422450682926131235350749310
163 => 262537412640768743.99999999999925007259524857371873

use strict;
use warnings;
use Math::BigFloat;

my $diff = 0.0001;
my $pi   = Math::BigFloat->bpi(50);

sub exp_pi_sqrt ($) {
    return exp( $pi * Math::BigFloat->new( $_[0] )->bsqrt );
}

print( $_->[0], " => ", $_->[1], "\n" ) for 
    grep {
           abs( $_->[1] - $_->[1]->copy->bfloor )  <= $diff
        || abs( $_->[1] - $_->[1]->copy->bceil  )  <= $diff
    }
    map { [ $_ => exp_pi_sqrt($_) ] } 1 .. 200;

どう書く.orgさんの趣旨に沿っているか分かりませんが、プログラマブル電卓を使用しました☆
使い方は、コードをコピペして
　編集→クリップボードを実行
です☆

see: 多倍長電卓LM

1
2
3

for(n=1;n<200;n++)
  if((exp(pi*sqrt(n))*10000+1)%10000<2)
    printf("%d %.99f\n",n,exp(pi*sqrt(n)));

Gaucheだと任意精度の小数が扱えないので基本的に有理数で計算しています。
ただしまともにやるととんでもない時間がかかるので、以下のように計算をはしょっています。
(1) 問題は「整数に近いかどうか？」というものなので、exp(pi*sqrt(n))を級数展開した
ときの各項では小数部分のみ計算している。これで、各項の最終的な精度はdoubleで十分になる。
(2) 分母、分子の桁数が大きくなるのを防ぐため、偶数項と奇数項で計算を分けている
(偶数項ではnの平方根はいらなくなるから)。

とはいえ、これでも結構時間がかかっていて、MacPro 3GHzで1500秒程度でした。

(use srfi-1)

(define (calc-pi k)
  (do ((n 1 (+ n 1))
       (a 1 (* a (- (* 4 n n) (* 2 n))))
       (b 2 (* b 16))
       (c 1 (* c n))
       (p 0 (+ p (/ a (* b c c (- (* 2 n) 1))))))
      ((<= k n) (* 6 p))))

(define (%sqrt* a n)
  (let ((root (sqrt a)))
    (if (= (* root root) a)
        (inexact->exact root)
        (do ((i 1 (+ i 1))
             (xn a (/ (+ xn (/ a xn)) 2)))
            ((<= n i) xn)))))

(define sqrt* (cut %sqrt* <> 10))

(define pi* (calc-pi 50))

(define (fract/ n d)
  (define (str->dbl str)
    (string->number
     (string-append "0." (substring str 0 (min 15 (string-length str))))))
  (let ((dstr (number->string d))
        (nstr (number->string (modulo n d))))
    (let ((d* (str->dbl dstr))
          (n* (str->dbl nstr)))
      (/. n* d* (expt 10 (- (string-length dstr) (string-length nstr)))))))
    
(define (%fract-exp-pi* n epsilon)
  (let ((q1 (sqrt* n)))
    (let ((pin (numerator pi*))
          (pid (denominator pi*))
          (q1n (numerator q1))
          (q1d (denominator q1)))
      (define (odd-term i pn pd fract)
        (let* ((new-pn (* pn pin))
               (new-pd (* pd pid i))
               (af (fract/ (* new-pn q1n) (* new-pd q1d)))
               (new-fract (+ fract af)))
          (if (< 0.0 af epsilon)
              new-fract
              (even-term (+ i 1) new-pn new-pd new-fract))))
      (define (even-term i pn pd fract)
        (let* ((new-pn (* pn pin n))
               (new-pd (* pd pid i))
               (af (fract/ new-pn new-pd))
               (new-fract (+ fract af)))
          (if (< 0.0 af epsilon)
              new-fract
              (odd-term (+ i 1) new-pn new-pd new-fract))))
      (odd-term 1 1 1 0))))

(define fract-exp-pi* (cut %fract-exp-pi* <> 0.00001))
  
(print (filter (lambda (n)
                 (let1 x (fract-exp-pi* n)
                   (<= (abs (- x (round x))) 0.0001)))
               (iota 199 1)))

すみません、さっきのコードだと(pi * sqrt(n))^k / k! がたまたま整数に近かったときに
誤動作する可能性がありました(今回の解の探索範囲では問題ないのですが)。

--- eps.scm.orig        2007-09-16 00:31:43.000000000 +0900
+++ eps.scm     2007-09-16 00:32:29.000000000 +0900
@@ -39,9 +39,11 @@
       (define (odd-term i pn pd fract)
         (let* ((new-pn (* pn pin))
                (new-pd (* pd pid i))
-               (af (fract/ (* new-pn q1n) (* new-pd q1d)))
+               (an (* new-pn q1n))
+               (ad (* new-pd q1d))
+               (af (fract/ an ad))
                (new-fract (+ fract af)))
-          (if (< 0.0 af epsilon)
+          (if (and (< 0.0 af epsilon) (< an ad))
               new-fract
               (even-term (+ i 1) new-pn new-pd new-fract))))
       (define (even-term i pn pd fract)
@@ -49,7 +51,7 @@
                (new-pd (* pd pid i))
                (af (fract/ new-pn new-pd))
                (new-fract (+ fract af)))
-          (if (< 0.0 af epsilon)
+          (if (and (< 0.0 af epsilon) (< new-pn new-pd))
               new-fract
               (odd-term (+ i 1) new-pn new-pd new-fract))))
       (odd-term 1 1 1 0))))
@@ -60,3 +62,4 @@
                  (let1 x (fract-exp-pi* n)
                    (<= (abs (- x (round x))) 0.0001)))
                (iota 199 1)))

Squeak Smalltalk で。

組み込みの ScaledDecimal（Java でいうところの BigDecimal ）を使って、sqrt と exp を自前で用意。sqrt には、開平法を使用しています。1GHz PowerPC (OS X) で、21 分ほどかかりました(^_^;)。

| pi exp sqrt |

pi := 3.141592653589793238462643s24.

sqrt := [:x |
    | pairs digits mem rem |
    pairs := OrderedCollection new.
    [x > 0] whileTrue: [pairs addFirst: x \\ 100. x := x // 100].
    pairs addAll: (Array new: 24 withAll: 0).
    digits := rem := mem := 0.
    pairs do: [:pair |
        | next found |
        rem := rem * 100 + pair.
        (0 to: 9) findLast: [:i | rem - ((next := mem * 10 + (found := i)) * i) >= 0].
        rem := rem - (next * found).
        mem := next + found.
        digits := digits * 10 + found].
    (digits / 1e24) asScaledDecimal: 24].

exp := [:x |
    | val fac mul i delta |
    val := fac := i := 1s5.
    mul := x := x asScaledDecimal: 5.
    [(delta := mul / fac) > 0.00001s5] whileTrue: [
        val := val + delta.
        mul := mul * x.
        fac := fac * (i := i + 1)].
    val].

^(1 to: 200)
    collect: [:n | n -> (exp value: (pi * (sqrt value: n)))]
    thenSelect: [:assoc | (assoc value - (assoc value roundTo: 1)) abs <= 0.0001s4]

"=> {37 -> 199148647.99996s5. 
     58 -> 24591257751.99999s5. 
     67 -> 147197952743.99998s5. 
     163 -> 262537412640768743.99998s5} "

人のふんどしで相撲を取るシリーズ。

上記のコードを任意精度演算ライブラリとしてgmpを使いCに移植しました。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>

mpf_t curr, next, bunsi, bunbo;

void mpf_exp(mpf_t x)
{
    unsigned int k;

    mpf_set_ui(curr, 1.0);
    mpf_set_ui(bunsi, 1.0);
    mpf_set_ui(bunbo, 1.0);

    for (k = 1; ;k++) {
        mpf_mul(bunsi, bunsi, x);
        mpf_mul_ui(bunbo, bunbo, k);
        mpf_div(next, bunsi, bunbo);
        mpf_add(next, next, curr);
        if (mpf_cmp(curr, next) == 0) {
            mpf_set(x, curr);
            return;
        }
        mpf_set(curr, next);
    }
}

int main(int argc, char *argv[])
{
    char buf[128], tmp[8];
    double d;
    int i;
    mpf_t x, pi;

    mpf_set_default_prec(128);

    mpf_init(curr);
    mpf_init(next);
    mpf_init(bunsi);
    mpf_init(bunbo);
    mpf_init(x);
    mpf_init(pi);

    mpf_set_str(pi, "3.1415926535897932384626433832795029", 10);

    for (i = 1; i < 200; i++) {
        mpf_set_ui(x, i);
        mpf_sqrt(x, x);
        mpf_mul(x, x, pi);
        mpf_exp(x);

        gmp_sprintf(buf, "%.*Ff", 6, x);
        d = atof(strncpy(tmp, strchr(buf, '.'), 7));
        if ((d < 0.0001) || ((1 - d) < 0.0001))
            printf("%d %s\n", i, buf);
    }

    mpf_clear(curr);
    mpf_clear(next);
    mpf_clear(bunsi);
    mpf_clear(bunbo);
    mpf_clear(x);
    mpf_clear(pi);

    return 0;
}

便乗。このライブラリは重宝しそう。

function doukaku60(){ // 精度的に 58 が限界。
  with(Math) for(var r = [], f, n = 1; n < 200; n++)
    if(abs((f = exp(PI * sqrt(n))) - round(f)) < 1e-4) r.push([n, f]);
  return r;
}
// ↓ Apfloat を用いて書き換え //

function doukaku60_(){
  var F = Packages.org.apfloat.Apfloat;
  var M = Packages.org.apfloat.ApfloatMath;
  for(var r = [], p = 25, h = new F(0.5, p), f, n = 1; n < 200; n++){
    f = M.exp(M.pi(p).multiply(M.sqrt(new F(n, p))));
    if(M.abs(f.subtract(f.add(h).floor())).floatValue() < 1e-4) r.push([n, f.toString(true)]);
  }
  return r;
}
print(doukaku60_().join('\n'))

上記のコードをctypesを使ってpythonで実装してみました。

Cのヘッダファイル内で#ifdefやら#defineが多用されているので、もしかしたら
環境依存になっているかもしれません。

呼び出し毎にインスタンスを生成するexp関数と、静的なインスタンスを持つ
exp2関数の速度を比較すると、exp2の方が20%ほどはやいようです。

from ctypes import *

gmp = cdll.LoadLibrary('libgmp.so')

for s in ['set_default_prec', 'init', 'clear', 'set', 'set_ui', 'set_str', 'mul'
, 'mul_ui', 'div', 'add', 'cmp', 'sqrt', 'printf', 'snprintf']:
  try:
    gmp.__setattr__(s, gmp.__getattr__('__gmpf_'+s))
  except:
    gmp.__setattr__(s, gmp.__getattr__('__gmp_'+s))

def mpf():
  class mpf_t(Structure):
    _fields_ = [("_mp_prec", c_int), ("_mp_size", c_int), ("_mp_exp", c_int), ("
_mp_d", c_void_p)]

    def __init__(self):
      Structure.__init__(self)
      gmp.init(byref(self))

    def __del__(self):
      try:
        gmp.clear(byref(self))
      except:
        pass

  return pointer(mpf_t())

def exp(rop, op):
  curr, bunsi, bunbo, next = mpf(), mpf(), mpf(), mpf()

  gmp.set_ui(curr, c_ulong(1))
  gmp.set_ui(bunsi, c_ulong(1))
  gmp.set_ui(bunbo, c_ulong(1))

  k = 1
  while True:
    gmp.mul(bunsi, bunsi, op)
    gmp.mul_ui(bunbo, bunbo, c_ulong(k))
    gmp.div(next, bunsi, bunbo)
    gmp.add(next, next, curr)
    if gmp.cmp(curr, next) == 0:
      gmp.set(rop, curr)
      return
    gmp.set(curr, next)
    k += 1

def exp2(rop, op):
  try:
    curr, bunsi, bunbo, next = exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next
  except:
    exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next = mpf(), mpf(), mpf(), mpf()
    curr, bunsi, bunbo, next = exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next

  gmp.set_ui(curr, c_ulong(1))
  gmp.set_ui(bunsi, c_ulong(1))
  gmp.set_ui(bunbo, c_ulong(1))

  k = 1
  while True:
    gmp.mul(bunsi, bunsi, op)
    gmp.mul_ui(bunbo, bunbo, c_ulong(k))
    gmp.div(next, bunsi, bunbo)
    gmp.add(next, next, curr)
    if gmp.cmp(curr, next) == 0:
      gmp.set(rop, curr)
      return
    gmp.set(curr, next)
    k += 1

def frac_f(x):
  buf = c_buffer(256)
  gmp.snprintf(buf, len(buf)-1, "%.Ff", x)
  return float(buf.value[buf.value.index('.'):])

def main():
  from time import time
  t = time()

  gmp.set_default_prec(128)

  x, pi = mpf(), mpf()

  gmp.set_str(pi, "3.141592653589793238462643383279502884197", 10)

  for i in range(1, 200):
    gmp.set_ui(x, c_ulong(i))
    gmp.sqrt(x, x)
    gmp.mul(x, x, pi)
#    exp(x, x)
    exp2(x, x)
    f = frac_f(x)
    if f < 0.0001 or (1-f) < 0.0001:
      gmp.printf("%d %.6Ff\n", i, x)

  print time() - t

if __name__ == '__main__':
  main()

20000以下で調べてみたら，

(37,199148647.9999780465518567665009238753359004336659)
(58,24591257751.9999998222132414695761923552658122276102)
(67,147197952743.9999986624542245068292613125786285081833)
(163,262537412640768743.9999999999992500725971981856888793538563)
(232,604729957825300084759.9999921715268564302785946808125512858845)
(652,68925893036109279891085639286943768.0000000001637386442092346075723290625709)
(719,3842614373539548891490294277805829192.9999872495660121875632701836570684449713)
(1169,44555719382988281777368496770130045948309444044.9999608028638684615024311524958053653676)
(1467,18095625621654510801615355531263454706630064771074975.9999999901236936712413276522472419790897)
(2608,4750778730825177725463920948909726618214491718039471366318747406368792.0000003084643221299811801879962000157848)
(4075,1247257156019637304856107520018074552566824585862995272173368815794085495792299621093743.9999936541874689715769690801805608661462)
(5773,46309587632860353087565367053742331250287153098248715578209888177688338779879045292937243508078581367989.9999155688538141662485158867634827239495)
(5868,327451666639079200503292535866541250265248788274691526825971156747731856100971255480468836963064283775072.0000971752541625920841201776565931010652)
(14370,35853082260707987565058966569844138230073788113687426211111838655289387060751425253235617998546754095483027646820100689677619375065141669969772221012766669726610356.0000809961676073067730760452585504610564)
(19183,932865712335864748985892137487224407024879778525651187858506512365813448395651706659915858370562540302669863675912417090372354324023009941419688829317421179304232837551154837099435561184585.9999366484683843100572967742009198686654)

という結果でした．163が題意の意味で一番整数に近いようです．
所要時間は，time ではかって，

2128.54s user 8.40s system 96% cpu 36:49.42 total

でした．

module Main (main) where

import Data.List
import System.Environment
import System.IO
import CReal

default (CReal)

f :: CReal -> (Int, CReal)
f n = (ceiling n, exp $ pi * sqrt n)

g :: CReal -> Bool
g n = 0.0001 > h n

h :: CReal -> CReal
h n = if 0.5 > n' - n then n' - n else 1 - (n' - n)
  where n' = fromInteger $ ceiling n

main :: IO ()
main = do { args <- getArgs
          ; case args of
              []  -> loop 100 0
              n:_ -> loop (read n) 0
          }
loop n m 
 | n == m = return ()
 | True   = mapM_ (k . f) [(m*100)+1 .. (m+1)*100] >> loop n (m+1)

k :: (Int,CReal) -> IO ()
k n@(i,r) = if g r then hPutStrLn stdout (show n) >> hFlush stdout
	    else putStr (show i ++ "\r") >> hFlush stdout

DouKaku?

exp(pi * sqrt(n))が整数に近くなるnを探す

Posted feedbacks - Flatten