DouKaku?beta

[topic] 分数の発見

schemeにはrationalizeという手続きがありますが、これは結果を一つしか返しません。そして、複数の結果が欲しい場合も稀にあります。

そこで、非負の実数aが一つ与えられたときに、以下の条件を満たす分数b/cをaに近い順に全て表示する手続きを考えてみてください。条件は、 1. 分数b/cよりaに近い分数d/cは存在しない 2. 分数b/cは既約 3. cは1桁の整数 です。

例をいくつかあげます(あってると思うけど...)。

a = 1.732051 12/7 7/4 16/9 5/3 9/5 3/2 2/1

a = 3.141593 22/7 25/8 19/6 28/9 16/5 13/4 3/1

a = 1920 / 1080 16/9 9/5 7/4 11/6 12/7 5/3 2/1

post comment

Posted feedbacks - Nested

Flatten Hidden

この場合はあまり役に立ってませんが、分数クラスを定義してみました。 こんな感じでどうでしょうか？

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Sample176 {

    public static List<Rational> rationalizes(final double d) {
        final Map<Rational, Double> map = new HashMap<Rational, Double>();
        for (int index = 1; index < 10; index++) {
            int nume = (int) (d * index);
            Rational r1 = new Rational(nume, index);
            Rational r2 = new Rational(nume + 1, index);
            if ((r2.doubleValue() - d) - (d - r1.doubleValue()) > 0) {
                map.put(r1, Math.abs(r1.doubleValue() - d));
            } else {
                map.put(r2, Math.abs(r2.doubleValue() - d));
            }
        }
        List<Rational> result = new ArrayList<Rational>(map.keySet());
        Collections.sort(result, new Comparator<Rational>() {
            public int compare(Rational o1, Rational o2) {
                return Double.compare(map.get(o1), map.get(o2));
            }
        });
        return result;
    }

    public static void main(String[] args) {
        System.out.println(rationalizes(1.732051));
        System.out.println(rationalizes(3.141593));
        System.out.println(rationalizes((double) 1920 / 1080));
    }

}

class Rational extends Number {
    public final int numerator;
    public final int denominator;

    public Rational(int p, int q) {
        int nume = p;
        int deno = q;
        if (deno == 0) {
            throw new IllegalArgumentException("#N/A");
        } else if (nume == 0) {
            deno = 1;
        } else if (deno < 0) {
            nume *= -1;
            deno *= -1;
        }
        int r = gcd(nume, deno);
        numerator = nume / r;
        denominator = deno / r;
    }

    private int gcd(int n, int m) {
        int p = n;
        int q = m;
        int r = p % q;
        while (r > 0) {
            p = q;
            q = r;
            r = p % q;
        }
        return q;
    }


    @Override
    public int intValue() {
        return numerator / denominator;
    }
    @Override
    public long longValue() {
        return this.intValue();
    }
    @Override
    public float floatValue() {
        return (float) this.doubleValue();
    }
    @Override
    public double doubleValue() {
        return (double) numerator / denominator;
    }


    @Override
    public int hashCode() {
        return numerator + denominator * 7;
    }
    @Override
    public boolean equals(Object obj) {
        if (this == obj) return true;
        if (!(obj instanceof Rational)) return false;
        Rational other = (Rational) obj;
        return this.denominator == other.denominator
            && this.numerator == other.numerator;
    }

    @Override
    public String toString() {
        return String.format("%d/%d", numerator, denominator);
    }
}

表示が不完全だけど、リストを取り出すだけならば
>(rationalize-list 1.732051)
(12/7 7/4 16/9 5/3 9/5 3/2 2)

分数の時が不完全な表示になります。
>(print-rationalize-list 1920/1080 )
16/9 => 16/9 9/5 7/4 11/6 12/7 5/3 2/1 
NIL

単にニュートン法の応用です。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
(defun check (given fraction min max)
  (let ((center (floor (/ (+ min max) 2))))
    (if (or (equal center min) (equal center max))
    (/ (check-near given fraction (cons min max)) fraction)
    (if (> (* given fraction) center)
        (check given fraction center max)
        (check given fraction min center)))))

(defun check-near (given fraction list)
  (if (> (- given (/ (car list) fraction))
     (- (/ (cdr list) fraction) given))
      (cdr list)
      (car list)))

(defun rationalize-list (given)
  (let (list)
    (loop for i from 1 to 9 do
     (let ((ans (check given i 0 (* (1+ i) given))))
       (and (not (find ans list))
        (push ans list))))
    (sort list
      #'(lambda(x y)(< (abs (- x given))
               (abs (- given y)))))))

(defun print-rationalize-list (given)
  (let ((list (rationalize-list given)))
    (format t "~s => " given)
    (loop for v in list do
     (multiple-value-bind (f r) (floor v)
       (if (zerop r)
           (format t "~d/1 " v)
           (format t "~d " v))))
       (format t "~%")))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def gcd(m, n):
    assert m > 0
    assert n > 0
    while n:
        m, n = n, m % n
    assert m > 0
    return m

class Rational:
    def __init__(self, b, c):
        g = gcd(b, c)
        self.b = b // g
        self.c = c // g

    def __repr__(self):
        return "%d/%d" % (self.b, self.c)

    def __cmp__(self, that):
        return self.b * that.c - self.c * that.b

    def value(self):
        return float(self.b) / self.c

def rationalize(a):
    rationals = []
    for c in xrange(1, 10):
        b = int(a * c + 0.5)
        r = Rational(b, c)
        if r.b == b:
            rationals.append(r)
    rationals.sort(key=lambda r: abs(a - r.value()))
    return rationals

def solve(a):
    print "a = %f" % a, " ".join(repr(r) for r in rationalize(a))

def main():
    solve(1.732051)
    solve(3.141593)
    solve(1920.0 / 1080)

if __name__ == '__main__':
    main()
boost::math::gcd で既約分数のチェックは楽できました。 それ以外は普通？ 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#include <iostream>
#include <vector>
#include <functional>
#include <cmath>
#include <iterator>
#include <string>
#include <stdexcept>

#include <boost/math/common_factor_rt.hpp>
#include <boost/lexical_cast.hpp>

class rational_t
{
  int n_, d_;
  double r_;

  public:
  rational_t(int numerator, int denominator)
    : n_(numerator), d_(denominator),
    r_(static_cast<double>(n_)/static_cast<double>(d_)) {}

  double real() const
  { return r_; }
  bool isnan() const
  { return std::isnan(r_); }

  int numerator() const
  { return n_; }
  int denominator() const
  { return d_; }
};
std::ostream&
operator<<(std::ostream& os, const rational_t& r)
{
  return (os << r.numerator() << "/" << r.denominator());
}

class compare_near_to
: public std::binary_function<bool, rational_t, rational_t>
{
  double g_;
  public:
  compare_near_to(double goal) : g_(goal) {}

  bool operator()(
      const rational_t& lhs,
      const rational_t& rhs
      )
  {
    return std::fabs(g_-lhs.real()) < std::fabs(g_-rhs.real());
  }
};

rational_t
find_rational(
    double  goal,
    int     denominator
    )
{
  double min_diff = goal;

  rational_t rat(0, denominator);
  for ( int n = 1; ; ++n ) {
    rational_t rat_(n, denominator);
    double diff = std::fabs(goal - rat_.real());

    if ( diff > min_diff ) break;

    min_diff = diff;
    rat = rat_;
  }
  if ( boost::math::gcd(rat.numerator(), rat.denominator()) != 1 )
    return rational_t(0,0); // causes real is NaN
  return rat;
}

int main(int c, char** v)
{
  if ( c != 2 ) {
    std::cout << v[0] << " <real number|rational number>\n";
    return 0;
  }
  double goal;
  try {
    std::string val(v[1]);
    std::string::size_type idx;
    if ( (idx = val.find("/")) != std::string::npos )
      goal = boost::lexical_cast<double>(val.substr(0,idx)) /
            boost::lexical_cast<double>(val.substr(idx+1,val.length()-(idx+1)));
    else
      goal = boost::lexical_cast<double>(v[1]);

    if ( goal <= 0.0 )
      throw std::runtime_error("input number must be positive");
  }
  catch (const std::exception& e) {
    std::cerr << e.what() << "\n";
    return -1;
  }

  std::vector<rational_t> rationals;
  for ( int i = 1; i < 10; ++i )
    rationals.push_back(find_rational(goal, i));

  rationals.erase(
      std::remove_if(
        rationals.begin(),
        rationals.end(),
        std::mem_fun_ref(&rational_t::isnan) ),
      rationals.end());
  std::sort(
      rationals.begin(),
      rationals.end(),
      compare_near_to(goal) );

  std::copy(
      rationals.begin(),
      rationals.end(),
      std::ostream_iterator<rational_t>(std::cout, ", ") );

  return 0;
}

同じ値の分数かどうかを「入力との差が等しいかどうか」で判定しているので、入力の値によっては誤差が出るかも……

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def rationalizes value
  # [入力との差, [分母, 分子]]からなる配列を作成する
  candidates = (1..9).map do |moderator|
    numerator  = (value*moderator).round
    difference = (value-numerator.quo(moderator)).abs
    [difference, [moderator, numerator]]
  end
  
  # 入力との差、分母の順に比較してソートすることで、
  # 入力との差が等しい場合には分母の小さい分数が先に並ぶ
  # その上で、差が等しい場合は最初の分数のみを選択することで、
  # 既約分数のみが返却される
  latest_difference = nil
  candidates.sort.inject([]) do |result, (difference, (moderator, numerator))|
    unless difference == latest_difference
      latest_difference = difference
      result << "#{numerator}/#{moderator}"
    end
    result
  end
end

p rationalizes(1.732051)      #=> ["12/7", "7/4", "16/9", "5/3", "9/5", "3/2", "2/1"]
p rationalizes(3.141593)      #=> ["22/7", "25/8", "19/6", "28/9", "16/5", "13/4", "3/1"]
p rationalizes(1920 / 1080.0) #=> ["16/9", "9/5", "7/4", "11/6", "12/7", "5/3", "2/1"]

「差が等しい」だと符号が違う値を同じ分数だとみなしてしまうので、値自体を比較に使用するように修正しました。

また、条件1には抵触しませんが「近い順に全て」という問いなので、「差が等しく」「分母が同じで」「分子が異なる」分数が検知できないのは間違いだと考えて修正しました。

※例えば0.7を入力した場合、分母が5の分数は4/5しか見つけられないが、3/5も4/5と同じだけ近いので検出されるべきではないか？　という話

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def rationalizes value
  # [浮動小数点数, [分母, 分子]]からなる配列を作成する
  candidates = (1..9).inject([]) do |result, moderator|
    # 四捨五入すれば分母の近似値が求められるが、
    # .5の場合は両方を使用する
    numerator = value * moderator
    if numerator%1.0 == 0.5
      numerators = [numerator.floor, numerator.ceil]
    else
      numerators = [numerator.round]
    end
    
    numerators.each do |numerator|
      float = numerator.quo(moderator)
      result << [float, [moderator, numerator]]
    end
    result
  end
  
  # 入力との差、分母、分子の順に比較してソートすることで、
  # 入力との差が等しい場合には分母、分子の小さい分数が先に並ぶ
  candidates = candidates.sort_by do |float, fraction|
    [(value-float).abs, fraction]
  end
  
  # 浮動小数点数が等しい場合は最初の分数のみを選択することで、
  # 既約分数のみが返却される
  latest_float = nil
  candidates.inject([]) do |result, (float, (moderator, numerator))|
    unless float == latest_float
      latest_float = float
      result << "#{numerator}/#{moderator}"
    end
    result
  end
end

p rationalizes(1.732051)      #=> ["12/7", "7/4", "16/9", "5/3", "9/5", "3/2", "2/1"]
p rationalizes(3.141593)      #=> ["22/7", "25/8", "19/6", "28/9", "16/5", "13/4", "3/1"]
p rationalizes(1920 / 1080.0) #=> ["16/9", "9/5", "7/4", "11/6", "12/7", "5/3", "2/1"]

p rationalizes(0.7)           #=> ["5/7", "2/3", "3/4", "3/5", "4/5", "1/2", "1/1"]

分数が組み込みなので簡単。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
(use srfi-1)

(define (sort-by xs cmp-fn key-fn)
  (map car (sort (map (lambda (x) (cons x (key-fn x))) xs)
                 (lambda (a b) (cmp-fn (cdr a) (cdr b))))))

(define (main args)
  (let ((n (string->number (cadr args))))
    (for-each
     print
     (sort-by (delete-duplicates
               (map (lambda (m)
                      (/ (inexact->exact (round (* m n))) m))
                    (iota 9 1)))
              <
              (lambda (x) (abs (- x n)))))))

Index

Feed

Other

Link

Pathtraq

loading...
Creative Commons License このサイトの内容は、 クリエイティブ・コモンズ・ライセンスの下でライセンスされています。 [詳細]