Как найти (оптимальное) целочисленное соотношение разной точности?
Если у меня есть переменная m типа uint32 а также r типа uint16 а также постоянная float64 со значением т.е. f=0.5820766091346741, Как мне найти m,r которые удовлетворяют f=r/m?
Похож на Fraction.limit_denominator из python.
Этот репозиторий Github содержит различные алгоритмы наилучшего рационального приближения, но ограничивает только знаменатель.
4 ответа
Как мне найти
m,rкоторые удовлетворяютf=r/m?
= подразумевает точное.
Чтобы сделать это точно, если это возможно, см. Ниже. Этот подход не пытается наилучшим образом соответствовать, если точное решение не возможно, поскольку это не будет удовлетворять f=r/m,
Все конечные значения с плавающей точкой являются точными. "0.5820766091346741" сохранен в f может дать f близкое значение, но значение в f это точно.
Учитывая основание числа с плавающей запятой (очень часто 2), все они могут быть представлены точно с помощью: "целое число /(базовыйпоказатель)".
С двоичным 64 самый большой необходимый показатель составляет около (1023 + 53).
Поскольку OP хочет, чтобы результат соответствовал 32-битному r и 16-битный mвполне понятно, что большинство float64 (64-бит) не будет иметь точного решения - просто недостаточно комбинаций для сохранения результата.
Алгоритм ниже в прокомментированном C, предполагая основание 2.
// return true on success
bool fraction(double d, uint32_t *r, uint16_t *m) {
if (d < 0.0 || isnan(d) || d > UINT32_MAX) {
return false;
}
// Scale d to extract, hopefully a 32+15 bit integer
uint16_t power_of_2 = 32768; // largest power-of-2 in m
d *= power_of_2;
uint64_t ipart = (uint64_t) d;
// Even after scaling, `d` has a fractional part.
if (d != ipart) {
return false; // value has unrepresentable precision.
}
// while big and even, reduce the fraction
while (ipart > UINT32_MAX && (ipart % 2 == 0)) {
power_of_2 /= 2;
ipart /= 2;
}
// If reduction was insufficient ...
if (ipart > UINT32_MAX) {
return false; // value has unrepresentable precision.
}
*r = (uint32_t) ipart;
*m = power_of_2;
return true; // Success!
}
Простой ответ будет:
ROUND(f * 10^8)
f = ----------------
10^8
Затем вы можете реализовать небольшой цикл, который пытается разделить числитель и знаменатель на простые числа (начиная с 2 и выше). Что-то вроде (конечно, код не проверен):
var m = f * 10^8 ;
var r = 10^8 ;
var Prime_Numbers = [2,3,5,7,11,13,17,19,....] ;
for (var I = 0 ; I < Prime_Numbers.length ; I++) {
if ((Prime_Numbers[I] > m) ||
(Prime_Numbers[I] > r) ) {
break;
}
if (((m % Prime_Numbers[I]) == 0) &&
(r % Prime_Numbers[I]) == 0) ) {
m = m / Prime_Numbers[I] ;
r = r / Prime_Numbers[I] ;
}
console.log("Best m is: " + m) ;
console.log("Best r is: " + r) ;
:
:
}
Теперь вопрос состоит в том, сколько первичных номеров я должен включить в список?
Трудно сказать, но интуитивно не слишком много... Я бы сказал, что это будет зависеть от того, насколько вы строги в отношении ОПТИМАЛЬНО.
Надеюсь, что это дает вам некоторое направление.
Ура!!
РЕДАКТИРОВАТЬ:
Подумав немного дальше, чтобы всегда получать значения ABSOLUTE OPTIMAL, вам нужно включить все первичные числа вплоть до половины максимального значения, которое вы хотите, в качестве точности. Например, если точность тура должна быть 8 цифр (99999999), вам необходимо включить все первичные числа до (99999999/2).
РЕДАКТИРОВАТЬ 2:
Добавлено условие выхода в цикле.
Есть статья Дэвида Т. Эшли и др. который предлагает алгоритм для нахождения рационального приближения двумя целыми числами с разной точностью.
Я реализовал базовую версию, которая не содержит всей сложности упомянутой статьи 1.
Основная идея состоит в том, чтобы преобразовать число с плавающей точкой в непрерывную дробь, а затем искать конвергент высшего порядка, который удовлетворяет ограничениям. Смотрите вики для ознакомления с конвергентами.
Однако упомянутая статья описывает более сложный подход к применению ограничений на целочисленные соотношения (см. Раздел 5), который использует аналогию с решеточными структурами 1.
Я не даю вам алгоритм, потому что, IMO, продолжение дроби - верный путь.
Но я хотел показать, насколько хорошо это представление с плавающей запятой подходит для 64-битного IEEE754. Так что я немного поиграл с концепцией в Smalltalk (Squeak 64 bit).
Существует только 48 бит для представления r/m, причем многие комбинации представляют одно и то же значение (1/1=2/2=... 1/2=2/4=3/6=...), хотя уже есть 2^53 различных 64бит плавают в интервале [0,5,1,0). Таким образом, мы можем сказать, что большую часть времени мы не собираемся точно соответствовать f. Задача состоит в том, чтобы найти пару (r/m), которая округляет ближайшую к f.
Я не могу разумно играть с 48 битами, но могу с половиной точности и собрать все комбинации uint8/uint16:
v := Array new: 1<<24.
0 to: 1<<8-1 do: [:r |
0 to: 1<<16-1 do: [:m |
v at: (m<<8+r+1) put: ([r asFloat/m asFloat]
on: ZeroDivide do: [:exc | exc return: Float infinity])]].
s := v asSet sorted.
s size-2.
За исключением 0 и inf, это примерно 10 173 377 различных комбинаций из 16 777 216.
Меня интересует разрыв между двумя последовательными представимыми числами:
x := s copyFrom: 2 to: s size - 1.
y := (2 to: s size-1) collect: [:i | (s at: i) - (s at: i-1) / (s at: i) ulp].
минимум
u := y detectMin: #yourself.
около 2.71618435e8 ulp.
Давайте посмотрим, как формируются числитель и знаменатель:
p := y indexOf: u.
{((v indexOf: (x at: p)) - 1) hex.
((v indexOf: (x at: p-1)) - 1) hex}.
результат в #('16rFDFFFE' '16rFEFFFF') первые 4 цифры кодируют den (m), последние две цифры (r).
Таким образом, минимальный разрыв получается для
s1 := (1<<8-1) / (1<<8-1<<8-1).
s2 := (1<<8-2) / (1<<8-2<<8-1).
s2 asFloat - s1 asFloat / s2 asFloat ulp = u.
Это около значения 1/256 (или где-то рядом).
Мы можем предположить, что минимальный разрыв для 48-битной реперентации
s1 := (1<<16-1) / (1<<16-1<<16-1).
s2 := (1<<16-2) / (1<<16-2<<16-1).
s2 asFloat - s1 asFloat / s2 asFloat ulp.
Это около 16 ulp, не так уж и плохо, а максимальная плотность около 1/65536 (или где-то рядом).
Какова будет плотность около 0,5, как в вашем примере? Для 24-битного представления:
h := x indexOf: 0.5.
10133738. Давайте проверим точность в окрестности:
k := (h to: h +512) detectMin: [:i | (y at: i)].
u2 := y at: k.
Это 3.4903102168e10 ulp (примерно в 128 раз меньше плотности). Получается для:
s1 := (1<<8-1) / (1<<8-1<<1-1).
s2 := (1<<8-2) / (1<<8-2<<1-1).
s2 asFloat- s1 asFloat / s2 asFloat ulp = u2.
Таким образом, с 48 битами, мы можем ожидать плотность около
s1 := (1<<16-1) / (1<<16-1<<1-1).
s2 := (1<<16-2) / (1<<16-2<<1-1).
s2 asFloat- s1 asFloat / s2 asFloat ulp.
то есть 524320 ulp, или точность приблизительно 5.821121362714621e-11.
Редактировать: Как насчет наихудшей точности?
В зоне лучшей плотности:
q := (p-512 to:p+512) detectMax: [:i | y at: i].
{((v indexOf: (x at: q)) - 1) hex.
((v indexOf: (x at: q-1)) - 1) hex.}.
То есть #('16rFEFFFF' '16r10001')или, другими словами, непосредственно перед наилучшей точностью мы имеем локально худшее: w := y at: q. что составляет 6,8990021713e10 ulp для этих чисел:
s2 := (1<<8-1) / (1<<8-1<<8-1).
s1 := (1) / (1<<8).
s2 asFloat - s1 asFloat / s2 asFloat ulp = w.
В переводе на 48 бит, то есть около 1,048592e6 ulp:
s2 := (1<<16-1) / (1<<16-1<<16-1).
s1 := (1) / (1<<16).
s2 asFloat - s1 asFloat / s2 asFloat ulp.
И около 0,5 худшее составляет около 8,847936399549e12 ulp для 24 битов:
j := (h-512 to: h +512) detectMax: [:i | (y at: i)].
w2 := y at: j.
s2 := (1<<8-1) / (1<<8-1<<1-1).
s1 := (1) / (1<<1).
s2 asFloat- s1 asFloat / s2 asFloat ulp = w2.
или переведено в 48 бит, 3.4360524818e10 ulp:
s2 := (1<<16-1) / (1<<16-1<<1-1).
s1 := (1) / (1<<1).
s2 asFloat- s1 asFloat / s2 asFloat ulp.
Это примерно 3.814784579114772e-6 абсолютной точности, не очень хорошо.
Перед принятием такого представления было бы хорошо узнать, что является областью f, и знать о средней точности и точности наихудшего случая, достижимой в этой области.