Встала задачка такая...

Есть число от 0 до 2^24-1;

Нужно придумать алгоритм, который быстро считает десятичный логорифм от этого числа с точностью до 4 знаков после запятой. Можно создавать таблицы, но как можно меньше.

Доступные операции: + — / *.

признаюсь, честно, нужно по работе...

Здравствуйте, SergeyBi, Вы писали:


SB>Есть число от 0 до 2^24-1;
Только целые числа?

SB>Нужно придумать алгоритм, который быстро считает десятичный логорифм от этого числа с точностью до 4 знаков после запятой. Можно создавать таблицы, но как можно меньше.

Ряд Маклорена ?

M>Ряд Маклорена ?

Да, думаю очень даже подойдет Спасибо.

Здравствуйте, SergeyBi, Вы писали:

M>>Ряд Маклорена ?

SB>Да, думаю очень даже подойдет Спасибо.

Напрямую похоже не подойдет.
|x| > 1 так что ряд будет расходиться.
можно посчитаться -Ln(1/x). Тут ряд быстро должен сойтись. А для маленьких чисел например до 100 сделать табличку.

Здравствуйте, SergeyBi, Вы писали:

SB>Встала задачка такая...

SB>Есть число от 0 до 2^24-1;

SB>Нужно придумать алгоритм, который быстро считает десятичный логорифм от этого числа с точностью до 4 знаков после запятой. Можно создавать таблицы, но как можно меньше.

SB>Доступные операции: + — / *.

Битовые операции доступны (сдвиг вправо, проверка младшего)?
Если да, то легко посчитать двоичный логарифм, и потом домножить на 0.30103

Здравствуйте, MBo, Вы писали:
MBo>Если да, то легко посчитать двоичный логарифм, и потом домножить на 0.30103

Мда, про легкость двоичного логарифма я соврал.
Вот сразу десятичный (не более 46 итераций цикла)

function L10(Value: Double): Double;
var
  Table: array[1..15] of Double;
  curr, mult, next: Double;
  i: Integer;
begin
 //на самом деле эта таблица считается один раз и зашивается
  for i := 1 to 15 do
   Table[i] := Log10(1 + 1.0/(1 shl i));

  Result := 0;
  curr := 1;
  mult := 0.5;
  for i := 1 to 15 do begin
    while True do begin
      next := curr + curr * mult;
      if next > Value then
         break;
      Result := Result + Table[i];
      curr := next;
    end;
    mult := mult * 0.5;
  end;
end;

Здравствуйте, MBo, Вы писали:

MBo>Вот сразу десятичный (не более 46 итераций цикла)

Ыщщо одно
То была калька с двоичного логарифма, а модификация под десятичный втрое быстрее (а среднем 14 итераций внутреннего цикла):

function L100(Value: Double): Double;
var
  p10: array[1..16] of Double;
  l10: array[1..16] of Double;
  curr, next: Double;
  i: Integer;
begin
 //таблица считается один раз и зашивается
  for i := 1 to 16 do begin
    l10[i] := 1.0/(1 shl (i - 1)); 
    p10[i] := Power(10, l10[i]);
  end;

  Result := 0;
  curr := 1;
  for i := 1 to 16 do begin
    while True do begin
      next := curr * p10[i];
      if next > Value then
         break;
      Result := Result + l10[i];
      curr := next;
    end;
  end;
end;

Здравствуйте, SergeyBi, Вы писали:

SB>Есть число от 0 до 2^24-1;

... или, приблизительно, до 10^8

SB>Нужно придумать алгоритм, который быстро считает десятичный логорифм от этого числа с точностью до 4 знаков после запятой. Можно создавать таблицы, но как можно меньше.

4 знака после запятой — это таблица Брадиса.
Почему бы и не заколбасить её?

Идея такая: пусть p(k) = [f * 10^(k*0.0001)] — целочисленная таблица экспонент, k=[0..9999] — искомый диапазон логарифмов
f — коэффициент, при котором все экспоненты принимают целые значения
f >= 1/(10^0.0001-1) = 4342...., для красоты и удобства возьмём f = 10^4

Пусть исходное целое число X.
Найдём масштаб 10^b, при котором Y = X*10^b попадает в [p0;p10000)

Дальше — двоичным поиском — найдём позицию k числа Y в p[]
То есть, 10000*10^(k*0.0001) ~ X*10^b

10^(4-b+k*0.0001) ~ X
lg(X) = 4-b+k*0.0001

Собственно, код на питоне

from math import *

p = [int(10000 * 10**(k*0.0001)) for k in xrange(10001)]

def look(y) :
    if y < p[0] : raise Exception("too small")
    if y > p[-1] : raise Exception("too large")
    # TODO сделать двоичный, а не линейный поиск lower_bound
    for k in xrange(10000) :
        if p[k+1] > y :
            return k
    return 10000

def scale(x) :
    y,b = x,0
    while y < p[0] :
        y,b = y*10, b+1
    while x > p[0]*10 :
        y,b = y/10, b-1
    return y,b

def logo(x) :
    y,b = scale(x)
    k = look(y)
    return ((4-b)*10000+k) # остались в целочисленной арифметике

def lg(x) :
    return log(x)/log(10)

def delta(x) :
    return logo(x)-lg(x)*10000

Здравствуйте, SergeyBi, Вы писали:

SB>Есть число от 0 до 2^24-1;

SB>Нужно придумать алгоритм, который быстро считает десятичный логорифм от этого числа с точностью до 4 знаков после запятой. Можно создавать таблицы, но как можно меньше.

Делов-то! гуглить "fast logarithm", первая же ссылка дает годную пдф-ку.

Как-то давно я делал, даю как есть. Идея такова — вытаскиваем экспоненту и сразу получаем log2. Но крайне грубый. Для уточнения берем старшую честь мантиссы, столько бит, сколько у нас размер таблицы — для таблицы 256 — 8 бит. Чем больше таблица, тем точнее. А таблица вычисляется тоже как-то очень просто, если немного подумать. Но вот прямо сейчас не помню.

#include "stdio.h"
#include "math.h"
#include "time.h"


float Log2_Lut[64] = 
{
    0.0000000000f,0.0227201595f,0.0450880530f,0.0671144373f,
    0.0888095841f,0.1101833083f,0.1312449951f,0.1520036246f,
    0.1724677944f,0.1926457418f,0.2125453629f,0.2321742315f,
    0.2515396167f,0.2706484983f,0.2895075827f,0.3081233165f,
    0.3265019001f,0.3446492997f,0.3625712593f,0.3802733114f,
    0.3977607872f,0.4150388266f,0.4321123868f,0.4489862511f,
    0.4656650369f,0.4821532033f,0.4984550583f,0.5145747654f,
    0.5305163501f,0.5462837058f,0.5618805999f,0.5773106787f,
    0.5925774728f,0.6076844018f,0.6226347791f,0.6374318162f,
    0.6520786266f,0.6665782300f,0.6809335561f,0.6951474477f,
    0.7092226647f,0.7231618868f,0.7369677169f,0.7506426838f,
    0.7641892451f,0.7776097897f,0.7909066407f,0.8040820573f,
    0.8171382374f,0.8300773200f,0.8429013867f,0.8556124644f,
    0.8682125266f,0.8807034958f,0.8930872448f,0.9053655988f,
    0.9175403365f,0.9296131921f,0.9415858568f,0.9534599797f,
    0.9652371697f,0.9769189966f,0.9885069924f,1.0000000000f
};

unsigned char Log2_Lut_Int[256] = 
{
      0,  1,  3,  4,  6,  7,  9, 10, 11, 13, 14, 16, 17, 18, 20, 21,
     22, 24, 25, 26, 28, 29, 30, 32, 33, 34, 36, 37, 38, 40, 41, 42,
     43, 45, 46, 47, 49, 50, 51, 52, 54, 55, 56, 57, 59, 60, 61, 62,
     63, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 77, 78, 79, 80, 81,
     82, 84, 85, 86, 87, 88, 89, 90, 91, 93, 94, 95, 96, 97, 98, 99,
    100,101,103,104,105,106,107,108,109,110,111,112,113,114,115,117,
    118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,
    134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,
    150,151,152,153,153,154,155,156,157,158,159,160,161,162,163,164,
    165,166,167,167,168,169,170,171,172,173,174,175,176,176,177,178,
    179,180,181,182,183,184,184,185,186,187,188,189,190,191,191,192,
    193,194,195,196,196,197,198,199,200,201,202,202,203,204,205,206,
    206,207,208,209,210,211,211,212,213,214,215,215,216,217,218,219,
    219,220,221,222,223,223,224,225,226,227,227,228,229,230,230,231,
    232,233,233,234,235,236,237,237,238,239,240,240,241,242,243,243,
    244,245,246,246,247,248,248,249,250,251,251,252,253,254,254,255
};





enum { log2_shift = 8 };
unsigned short lookup_table[1 << log2_shift];

void fill_log_table()
{
    unsigned n = 1 << log2_shift;

    float f = 1.0f;
    float d = 1.0f / float(n - 1);
    for(unsigned i = 0; i < n; ++i)
    {
        lookup_table[i] = unsigned(logf(f) * 1.4426950408889634073599246810019f * 255.0 + 0.5); 
        if(i % 4 == 0)
            printf("%3d %f %3f\n", int(pow(1.0370, i) * 10), pow(1.0370, i), pow(1.0370, lookup_table[i]));
        //printf("%3d,\n", lookup_table[i]);
        f += d;
    }
}



float fast_log2(float v)
{
    union { float f; unsigned i; } fi;
    fi.f = v;
    int    exponent = ((fi.i >> 23) & 0xFF) - 127;
    float  mantissa = Log2_Lut[(fi.i & 0x7FFFFF) >> (23 - 6)];
    return mantissa + float(exponent);
}



int fast_log2_int(float v)
{
    union { float f; unsigned i; } fi;
    fi.f = v;
    int    exponent = (((fi.i >> 23) & 0xFF) - 127) << 8;
    int    mantissa = Log2_Lut_Int[(fi.i & 0x7FFFFF) >> (23 - 8)];
    return exponent + mantissa;
}


int main()
{
    //fill_log_table();
    float s = 1/256.0;// / 256.0f;

    printf("#dim 0.01 0.01 100000 100000\n");
    printf("#title Logarithmic scale\n");
    printf("Natural\n");
    printf("Approximate\n");
    printf("Error, %%\n");

/*
    while(s <= 65536.0f)
    {
        float ls = fast_log2(s);// * 5.7707801635558536294396987240209f;
        printf("%f %f %f %f\n", s, s, pow(2, ls), fabs(s - pow(2, ls)) / s * 100);
        //s *= 1.1892071150027210667174999705605f;
        s *= 1.0905077326652576592070106557607;
    }
*/


    while(s <= 65536.0f)
    {
        float ls = fast_log2_int(s) / 256.0f;// * 5.7707801635558536294396987240209f;
        printf("%f %f %f %f\n", s, s, pow(2, ls), fabsf(s - pow(2, ls)) / s * 100);
        s *= 1.1892071150027210667174999705605f;
        //s *= 1.0905077326652576592070106557607f;
        //s *= 1.009;
    }


/*
    unsigned i;
    s = 0;
    float d = 1.1892071150027210667174999705605f;
    for(i = 0; i < 1000000000; ++i)
    {
        s += fast_log2_int(d);// * (1.0f / 256.0f);
        d += 1.1892071150027210667174999705605f;
    }
    printf("%f\n", s);
    printf("%d\n", clock());
*/

}

Здравствуйте, мыщьхпыщьх, Вы писали:

Нуууу. Если есть плавающая арифметика, это, считай, полдела сделано. Да и логарифм наверняка найдётся — не в FPU так в библиотеке эмуляции.

А вот чисто на фиксированной арифметике (на целочисленной)?

Развитие темы с табличным методом.

Пусть у нас есть таблица экспонент для интервала [0.,1.), двоичным поиском в которой мы находим логарифм для [1.,10.)
Понятное дело, количество элементов обратно точности: для 4 знаков это 10_000, для 5 — 100_000... в общем, много.

Заметим, однако, что экспонента на этом интервале растёт очень плавно, и её можно аппроксимировать кусочно-линейной функцией.

Дальше всё просто:
Пусть powers[] — таблица экспонент: powers[k] = 10^(a + k*10^-a), где a — точность (accuracy), количество знаков
Пусть xs[],ys[],vs[] — таблица точек изломов:
— xs[j] = k,
— ys[j] = powers[k],
— vs[j] = powers[k+1]-powers[k] = j+2

lg(x) для x в [10^a,10^(a+1)] находится следующим образом
— j = lower_bound(ys,x)
— d = x-ys[j]
— k = xs[j] + d/vs[j]
— lg(x) = k

1/lg'(1.) = 2.302, 1/lg'(10.) = 23.02 — так что...
Ура, товарищи!
Мы умеем находить логарифм с помощью таблицы из 23 элементов!
(На самом деле, чем выше точность, тем больше отрезков придётся вводить: скажем, не по +1, а по +0.1 — и, соответственно, масштабировать содержимое таблицы).

Здравствуйте, Кодт, Вы писали:

К>А вот чисто на фиксированной арифметике (на целочисленной)?

А чего, можно и на фиксированной, если пораскинуть мозгом. Грубый log2 для целого числа — это просто номер его старшего бита. Вычисляется за 2 ифа с таблицей. Для уточнения можно использовать тот же принцЫп — взять старшие биты после самого старшего и корректирующую таблицу.
Для получения log10 — еще одно умножение во float, или умножение/деление в fix-point.
А, кстати, а что на выходе? — надо 4 значащих цифры — в каком виде? float или fix-point?

Здравствуйте, Кодт, Вы писали:

К>lg(x) для x в [10^a,10^(a+1)] находится следующим образом
К>- j = lower_bound(ys,x)
К>- d = x-ys[j]
К>- k = xs[j] + d/vs[j]
К>- lg(x) = k

К>1/lg'(1.) = 2.302, 1/lg'(10.) = 23.02 — так что...
К>Ура, товарищи!
К>Мы умеем находить логарифм с помощью таблицы из 23 элементов!
К>(На самом деле, чем выше точность, тем больше отрезков придётся вводить: скажем, не по +1, а по +0.1 — и, соответственно, масштабировать содержимое таблицы).

Я тут человек новый, может скажу чего не то, но по-моему, все то же самое можно получить безо всяких там lowerbound, прямой адресацией. Возможно придется задействовать 2 таблицы. По тому же принципу, что и в статье про быстрый логарифм.

Здравствуйте, мыщьхпыщьх, Вы писали:

М>Я тут человек новый, может скажу чего не то, но по-моему, все то же самое можно получить безо всяких там lowerbound, прямой адресацией. Возможно придется задействовать 2 таблицы. По тому же принципу, что и в статье про быстрый логарифм.

А ведь действительно, грубо говоря, на отрезке x = [1.000, 10.000) значение lg(x) = [0.000,1.000) с точностью до 4 разрядов определяется не более чем 4 старшими разрядами x. Отбрасываем младшие и наслаждаемся (естественно, отмасштабировав до [1000, 9999]).

Буду вкуривать, как можно получить бОльшую точность без распухания таблицы.