paul@186 | 1 | #!/usr/bin/env python2 |
paul@186 | 2 | |
paul@186 | 3 | from math import ceil, floor |
paul@186 | 4 | import sys |
paul@186 | 5 | |
paul@186 | 6 | def is_integer(x): |
paul@186 | 7 | target = round(x) * 1000 |
paul@186 | 8 | #print x, target - 100, floor(x * 1000), target + 100 |
paul@186 | 9 | return target - 100 < floor(x * 1000) < target + 100 |
paul@186 | 10 | |
paul@186 | 11 | def getscale(x): |
paul@186 | 12 | scale = getscale_part(x) |
paul@186 | 13 | if is_integer(scale): |
paul@186 | 14 | return scale |
paul@186 | 15 | return scale * getscale(scale) |
paul@186 | 16 | |
paul@186 | 17 | def getscale_part(x): |
paul@186 | 18 | part = x - int(x) |
paul@186 | 19 | if not part: |
paul@186 | 20 | return 1 |
paul@186 | 21 | elif part > 0.5: |
paul@186 | 22 | return 1 / (1 - part) |
paul@186 | 23 | else: |
paul@186 | 24 | return 1 / part |
paul@186 | 25 | |
paul@245 | 26 | def get_divider_operands(frequency, source_frequency): |
paul@245 | 27 | ratio = float(frequency) / source_frequency |
paul@245 | 28 | scale = getscale(ratio) |
paul@245 | 29 | return scale * ratio, scale |
paul@245 | 30 | |
paul@245 | 31 | def reduce_divider_operands(m, n, m_limit, n_limit): |
paul@186 | 32 | while m > m_limit and n > n_limit and m > 1 and n > 1: |
paul@186 | 33 | m /= 2 |
paul@186 | 34 | n /= 2 |
paul@186 | 35 | return m, n |
paul@186 | 36 | |
paul@186 | 37 | def show(s, m, n, f_label="f"): |
paul@186 | 38 | print "m =", m |
paul@186 | 39 | print "n =", n |
paul@186 | 40 | print "%s' =" % f_label, s * m / n |
paul@186 | 41 | print "%s' ~=" % f_label, s * round(m) / round(n) |
paul@186 | 42 | |
paul@186 | 43 | if len(sys.argv) < 4: |
paul@186 | 44 | f, s = map(float, sys.argv[1:3]) |
paul@245 | 45 | m, n = get_divider_operands(f, s) |
paul@245 | 46 | m, n = reduce_divider_operands(m, n, 0x200, 0x100000) |
paul@186 | 47 | show(s, m, n) |
paul@186 | 48 | |
paul@186 | 49 | else: |
paul@186 | 50 | f, s, i_min, i_max = map(float, sys.argv[1:5]) |
paul@186 | 51 | |
paul@186 | 52 | # Constraints: |
paul@186 | 53 | # |
paul@186 | 54 | # i_min < s * m' / di < i_max |
paul@186 | 55 | # i_min / s < m' / di < i_max / s |
paul@186 | 56 | # |
paul@186 | 57 | # f = s * m' / (di * do) |
paul@186 | 58 | # f = (s * m' / di) / do |
paul@186 | 59 | # f * do = s * m' / di |
paul@186 | 60 | # |
paul@186 | 61 | # i_min < f * do < i_max |
paul@186 | 62 | # i_min / f < do < i_max / f |
paul@186 | 63 | # |
paul@186 | 64 | # f = s * m / n |
paul@186 | 65 | # |
paul@186 | 66 | # s * m' / (di * do) = s * m / n |
paul@186 | 67 | # m' / (di * do) = m / n |
paul@186 | 68 | # m' = (m / n) * (di * do) |
paul@186 | 69 | # m' / di = (m / n) * do |
paul@186 | 70 | |
paul@186 | 71 | do_min = ceil(i_min / f) |
paul@186 | 72 | do_max = floor(i_max / f) |
paul@186 | 73 | |
paul@186 | 74 | do = do_min |
paul@186 | 75 | |
paul@186 | 76 | print do_min, "<= do <=", do_max |
paul@186 | 77 | |
paul@186 | 78 | while do <= do_max: |
paul@246 | 79 | do0 = int(do ** 0.5) |
paul@186 | 80 | do1 = int(do / do0) |
paul@186 | 81 | if do0 * do1 == int(do): |
paul@186 | 82 | i = do * f |
paul@245 | 83 | m, n = get_divider_operands(i, s) |
paul@245 | 84 | m, n = reduce_divider_operands(m, n, 0x3f, 0x3f) |
paul@186 | 85 | if m <= 0x3f and n <= 0x3f: |
paul@246 | 86 | show(s, m, n, "i") |
paul@186 | 87 | print "do =", do |
paul@186 | 88 | print "do0 =", do0 |
paul@186 | 89 | print "do1 =", do1 |
paul@186 | 90 | print "f' =", s * m / (n * do) |
paul@246 | 91 | print "m =", round(m) |
paul@246 | 92 | print "n =", round(n) |
paul@186 | 93 | print "f' ~=", s * round(m) / (round(n) * round(do)) |
paul@186 | 94 | break |
paul@186 | 95 | else: |
paul@186 | 96 | print "Ignored: m =", m, "n =", n, "do0 =", do0, "do1 =", do1 |
paul@186 | 97 | do += 1 |
paul@186 | 98 | |
paul@186 | 99 | # vim: tabstop=4 expandtab shiftwidth=4 |