Войти
ПрограммированиеФорумОбщее

Трюки с float: быстрое вычисление логарифмов. (комментарии) (4 стр)

Страницы: 13 4 5 610 Следующая »
#45
2:16, 9 янв. 2022

> Кэхэн и компания как раз ратуют за разделение этих ошибок.
Увы, обоснования довольно скупые.

https://stackoverflow.com/questions/6430448/why-doesnt-gcc-optimi… aaaaa/6453476 :

IEEE-754 is pretty clear on the point that floating-point data do represent exact values

For numerical analysis, your brain will thank you if you interpret floating point numbers not as intervals, but as exact values (which happen to be not exactly the values that you wanted).

https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf :

Very few of us learn that roundoff in these formulas can lose more than the roots’ last sig. digit or two when the roots come too close to coincidence. Then the number of sig. digits lost roughly equals the number of leading digits in which the roots agree up to a maximum loss of about half the sig. digits carried by the arithmetic. The shock we feel at learning of this loss may abate when we learn from Backward Error-Analysis that the roots suffer little worse from roundoff during their computation than if first each of the given coefficients A, B and C had been perturbed by a rounding error or two when it was generated, and then the roots had been computed exactly from the perturbed coefficients. Thus we might conclude that the computed roots are scarcely more uncertain due to roundoff during their computation than to roundoff before their computation, and therefore the loss of more than a few of their last sig. digits, when this loss occurs, is deserved by “Ill Conditioned roots”, as they are called when they are hypersensitive to perturbations of the coefficients.
That conclusion is not quite right. It seems somewhat like blaming a crime upon its victim.
Good reasons exist for software to compute the roots each correct in every bit except perhaps its last one or two, as if every coefficient’s every bit had been supplied correctly. It may have been, for all the programmer and users of that software know. Software that accurate can be costly to program in a language that lacks support for arithmetic rather more precise than the precision of the coefficients and the desired accuracy of the roots. Why is their accuracy worth that cost?

Instead numerical troubles are usually recognized too late if ever, and then misdiagnosed at least at first if not forever. We tend to blame time wasted and other inflated costs upon that data by calling it “Ill-Conditioned” whenever we have persuaded ourselves that our computational methods, having succeeded on all the other data we tested, must be “Numerically Stable”. Name-calling makes the caller feel better without enlightening him.


#46
(Правка: 8:05) 8:04, 9 янв. 2022

}:+()___ [Smile]

Нафига деление, что в одном варианте, что в другом? Не будет ли достаточно одного полинома?

Я исследовал этот вопрос. Всё дело в скорости сходимости различных рядов. А деление, особенно в формате single, не такое дорогое, - примерно как 3 операции умножения. Можно разложить в ряд Чебышёва функцию f(t)=lb(t/3+1)/t, где t=3*(x-1), аргумент x приводится промежутку [2/3; 4/3]. В этом случае lb(x)=t*f(t)+k без деления, но тогда, во-первых, для достижения одинарной точности понадобится 10 членов разложения вместо 4-х, которые в моём варианте; во-вторых, понадобится ещё константа 4/3 для определения принадлежности x промежутку (у меня сравнение с 3/2 происходит бесплатно).

> Самое главное, вам нужно выбрать хорошую аппроксимацию.
Скорее, нужно перестать пользоваться однобоким определением точности, когда рассматривается только погрешность по y, а погрешность по x почему-то считается равной нулю. Для практических целей это просто бесполезная потеря производительности на ровном месте.

Нюанс, о котором вы говорите, для функции lb(x) не актуален, он характерен для тригонометрических функций (редуцировнание больших значений аргумента), и в какой-то мере - для экспоненты и некоторых других функций. Вы правы, что нет необходимости принимать какие-то экстраординарные меры для повышения точности в тех случаях, когда точность всё равно теряется при передаче аргумента в функцию, т.е. само x не может быть задано точно. Хотя какой-то запас по точности всё равно должен быть в разумных пределах.

#47
8:32, 9 янв. 2022

Aenigma
> Нюанс, о котором вы говорите, для функции lb(x) не актуален, он характерен для тригонометрических функций (редуцировнание больших значений аргумента)
Да вроде погрешность, указанная в #40, была бы примерно такая же, если бы у нас была погрешность не округления (т. е. считаем промежуточные вычисления точно), а в аргументе (для примера - тоже в 4-м разряде: log2(0.999)=-0.001443..., log2(0.9991)=-0.001299...).

#48
(Правка: 8:56) 8:53, 9 янв. 2022

Вблизи 1 у логарифма действительно есть эффект неустойчивости, когда небольшая относительная погрешность x приводит к большой относительной погрешности результата. Я написал, что не актуально, потому что проблема легко решается смещением области аппроксимации в окрестность 1, а игнорирование проблемы даёт неприемлемо большие ошибки. Т.е. тут почти нет вариантов, как действовать. Вот в случае, когда редуцируется аргумент тригонометрических функций, возможны различные подходы

#49
12:28, 9 янв. 2022

Aenigma
> Всё дело в скорости сходимости различных рядов.
Не, то что у логарифма ряд хреновый, это понятно. Но можно же попробовать альтернативные делению подходы, например, возвести аргумент в квадрат или четвертую степень; попробовать полиномы от простых полиномов и т. п.

> Я написал, что не актуально, потому что проблема легко решается смещением области аппроксимации в окрестность 1, а игнорирование проблемы даёт неприемлемо большие ошибки.
Если мы считаем результат в полную точность float/double, то эта проблема нерешаема. Соответственно, логично, что аппроксимация пониженной точности должна давать более-менее равномерное понижение относительно этой предельной точности. Т. е. эта ошибка приемлема. А когда неприемлема надо использовать log1p.

#50
13:54, 9 янв. 2022

}:+()___ [Smile]

Если мы считаем результат в полную точность float/double, то эта проблема нерешаема.

Нужно уточнить, о какой проблеме мы говорим. Если, например, пользователь хочет вычислить функцию, которая сводится к ln(x+1), то сделать это хорошо с помощью функции ln(x) он никак не сможет. Эта проблема не решаемая, вернее решаемая переходом к другой функции (кстати, коэффициенты моей аппроксимации lb(x) подходят и для аппроксимации функции lb(x+1)). Но это проблема не самой функции логарифм, а проблема корректного её использования. Сам логарифм может быть посчитан без потери точности.

Но можно же попробовать альтернативные делению подходы, например, возвести аргумент в квадрат или четвертую степень; попробовать полиномы от простых полиномов и т. п.

Нужно понимать, что даже 2-3 лишних команды могут иметь то же время работы, что и команда деления, поэтому усложнение кода вряд ли будет оправданным. Попробовать варьировать функцию можно, и есть вариант, для которого ряд сходится быстрее, чем используемый у меня - это переход к обратной функции (в смысле 1/f(t)). Но здесь добавляется ещё одно деление, поэтому такой вариант хорош только для двойной или более высокой точности и совместно с методом Паде-Чебышёва, где в любом случае есть деление.

#51
21:26, 9 янв. 2022

Небольшая зарисовка на предмет точности (вдохновлённая примером из https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf ) :
https://rextester.com/ZBI21312
Т. е. результат использование ill-conditioned функции - не обязательно ill-conditioned. Однако если забить на точность в пределах погрешности аргумента - результат неточный.

#52
0:23, 10 янв. 2022

FordPerfect
В общем, от деления у меня избавиться не получилось, но, как минимум, твой многочлен можно серьезно сократить:

    x=(x-1.5f)/(x+1.5f);
    float x2=x*x;

    return y+5.84962500721660e-01f
        +x *(2.88539316635122e+00f
        +x2*(9.61185690308810e-01f
        +x2*(6.07166866084299e-01f)));

или

    x=(x-1.5f)/(x+1.5f);
    float x2=x*x;

    return y+5.84962500720858e-01f
        +x *(2.88539005028199e+00f
        +x2*(9.61807113213581e-01f
        +x2*(5.76149707019795e-01f
        +x2*(4.42589808255434e-01f))));

Но эти многочлены, скорее всего, не очень хорошие, дают ступеньки при смыкании соседних интервалов.
Стоило бы сделать подгон многочлена, который точно смыкается сам с собой.

#53
1:11, 10 янв. 2022

}:+()___ [Smile]
Да, твои короче и, вроде, средняя точность не хуже, но log(1)==0 хочется сохранить точно.

#54
1:48, 10 янв. 2022

Вот такое имеет на [1;2) примерно ту же точность, что мой исходный:

    x-=1.0f;
    return y+x*( 1.44269473e+00f
            +x*(-7.21308896e-01f
            +x*( 4.80051449e-01f
            +x*(-3.53295347e-01f
            +x*( 2.55612642e-01f
            +x*(-1.54604686e-01f
            +x*( 6.29557284e-02f
            +x*(-1.21056225e-02f))))))));

Можно несколько улучшить.

#55
3:28, 10 янв. 2022
float log2_v2(float x)
{
    if(!(x>0.0f)) return (x==0.0f?-INFINITY:nanf(""));
    if(x==+INFINITY) return +INFINITY;
    uint32_t bits,bias;

    memcpy(&bits,&x,sizeof(float));
    bias=(bits+(1U<<22))>>23U;
    float y=float(int32_t(bias))-127.0f;
    bits=bits-((bias-127U)<<23);
    memcpy(&x,&bits,sizeof(float));
    x=(x-1.0f)/(x+1.0f);
    float x2=x*x;
    return y+x*(2.88539008e+00f
           +x2*(9.61799077e-01f
           +x2*(5.76633549e-01f
           +x2*(4.35105166e-01f))));
}

https://rextester.com/MUKE74391 (выкладки, F8 чтобы увидеть график)
https://rextester.com/TZXM73578 (тест)

#56
3:35, 10 янв. 2022
float log2_v3(float x)
{
    if(!(x>0.0f)) return (x==0.0f?-INFINITY:nanf(""));
    if(x==+INFINITY) return +INFINITY;
    uint32_t bits,bias;

    memcpy(&bits,&x,sizeof(float));
    bias=(bits+(1U<<22))>>23U;
    float y=float(int32_t(bias))-127.0f;
    bits=bits-((bias-127U)<<23);
    memcpy(&x,&bits,sizeof(float));
    x=(x-1.0f)/(x+1.0f);
    float x2=x*x;
    return y+x*(2.88539003e+00f
           +x2*(9.61808794e-01f
           +x2*(5.76110615e-01f
           +x2*(4.42826415e-01f))));
}

RMS: 0.64 ulp, Linf: 3.0 ulp.
Привожу к [0.75;1.5).

#57
4:12, 10 янв. 2022

Оптимальный диапазон походу \((1/\sqrt{2};\sqrt{2})\), но 4 члена уже упираются в предел float. Вот для трех:

    uint32_t bias=0x004AFB0D;  // bitcast(3-sqrt(2))&0x007FFFFFUL

    uint32_t bits;
    memcpy(&bits,&x,sizeof(float));
    bits+=bias;
    float y=float(bits>>23U)-127.0f;
    bits=((bits&0x007FFFFFUL)|0x3F800000UL)-bias;
    memcpy(&x,&bits,sizeof(float));

    x=(x-1)/(x+1);
    float x2=x*x;
    
    return y
        +x *(2.88539128929759e+00f
        +x2*(9.61470821312105e-01f
        +x2*(5.98973520973232e-01f)));
#58
4:24, 10 янв. 2022

}:+()___ [Smile]
> Вот для трех:
Насколько я понимаю, RMS: 1.7 ulp, Linf: 9.0 ulp.

#59
(Правка: 7:26) 7:26, 10 янв. 2022

FordPerfect, похоже вы получили вариант, близкий к оптимальному по точности и производительности для полной точности float.

}:+()___ [Smile] Вариант с 3-мя членами разложения - хороший компромисс между точностью и производительностью, если не нужна полная точность. Она, кстати, ненамного хуже полной точности.

Для повышения точности промежуток действительно лучше выбирать [0.5^0.5; 2^0.5], но по производительности при равном числе слагаемых  [0.75; 1.5] может оказаться лучше. Надо тестировать. К сожалению, на VS компилятор генерирует не лучший код, поэтому не могу адекватно замерить производительность и сравнить со своим кодом.

Страницы: 13 4 5 610 Следующая »
ПрограммированиеФорумОбщее