[Russian][Math] Как я снимал квартиру при помощи мат.анализа

Съем квартиры

Во многих крупных городах Украины, при съеме квартиры, нужно заплатить за первый месяц, последний месяц (залог) и половину отдать риелтору. Т.е., если вы видите в объявлении квартиру за 12000 гривен в месяц, при заключении договора нужно будет отдать 12000*2 + 12000/2 = 30000 гривен. Это: $2x + \frac{x}{2}$, где x = цена за месяц.

Пока я листал объявления, я подсчитал, что сразу нужно будет заплатить 30k UAH за квартиру, предлагаемую за 12k UAH в месяц. Потом нашел другое объявление, где квартира стоила 12500 в месяц. Сколько надо будет заплатить в начале? Считать в уме было лень (особенно, разделить 12500 на 2) и я задался вопросом, можно ли подсчитать это используя информацию о предыдущей квартире? Мы ведь понимаем, что синхронно с ростом цены "за месяц" возрастает и "итоговая" цена. Если цена (в месяц) возросла на 500 гривен, то на сколько нужно увеличить "итоговую" цену? Нетрудно догадаться, что с каждой гривной цены "в месяц", "итоговая" цена возрастает на 2.5 гривны. Если цена "в месяц" возросла на 500 гривен, значит "итоговая" возрастет на 1250 гривен. Если вместо 12000 у нас теперь 12500, то "итоговая" будет не 30000 (старый результат), а 31250 гривен (старый результат плюс 1250).

Чтобы подсчитать "итоговую" цену заново, нужны три действия - умножение, деление, сложение. Чтобы подсчитать в моем случае, нужно разницу умножить на 2.5 и прибавить её к старой "итоговой" цене. Если лениться запускать калькулятор, а вместо этого считать в уме, то подсчитать во втором случае может быть легче.

Так вот, 2.5 это и есть производная от ф-ции $2x + \frac{x}{2}$. Она показывает, как растет результат ф-ции при росте $x$.

$2x + \frac{x}{2} = 2.5x$, а производная ф-ции $nx$ это просто $n$.

Постройка графика ф-ции (или табулирование)

Нужно построить график ф-ции $7x$, либо же просто построить таблицу её значений (это и называется табулированием). Мы пишем что-то такое:

#!/usr/bin/env python3

for x in range(15):
    print (x, x*7)

И это работает. А как бы сделать быстрее? Если подумать? Мы можем вместо вычисления $7x$ каждый раз, вычислять разницу, а это 7. Ведь мы считаем с шагом 1 от 1 до 14, т.е., последовательно.

$7$ это производная от $7x$.

#!/usr/bin/env python3

a=0
for x in range(15):
    print (x, a)
    a = a + 7

Конечно, сложение работает быстрее умножения.

Теперь кое-что труднее: $x^2$:

#!/usr/bin/env python3

for x in range(15):
    print (x, x*x)

Это может быть неплохим вопросом на собеседовании - попробуйте заставить эту ф-цию работать быстрее. Возможно ли это вообще? Нам нужно не пересчитывать $x^2$ каждый раз, а на каждом шаге вычислять разницу.

В Wolfram Mathematica находим ф-цию разницы:

In[]:= f[x_] = x^2

In[]:= ExpandAll[f[x + 1] - f[x]]
Out[]:= 1 + 2x

(ExpandAll сокращает выражение, делая его как можно проще/короче.)

Кстати, ничего не напоминает? $2x$ это первая производная от $x^2$.

Пишем такое:

#!/usr/bin/env python3

a=0
for x in range(15):
    print (x, a)
    a = a + 2*x + 1

Работает так же. А вот быстрее ли? Перепишем на чистый Си. Для теста, будем считать сумму этих всех значений: $\sum_{i=1}^{...} i^2$

#include <stdint.h>
#include <stdio.h>
#include <time.h>

int main()
{
        uint64_t _max=0x1ffffffff;

        uint64_t sum=0;
        time_t start=time(NULL);

        for (uint64_t i=0; i<_max; i++)
        {
                sum=sum+i*i;
        };
        printf ("sum=%lx\n", sum);
        printf ("time=%ld\n", time(NULL)-start);

        sum=0;
        start=time(NULL);
        uint64_t prev=0;
        for (uint64_t i=0; i<_max; i++)
        {
                sum=sum+prev;
                prev=prev+(2*i+1);
        };
        printf ("sum=%lx\n", sum);
        printf ("time=%ld\n", time(NULL)-start);
};

В итоге:

sum=aaaaaaaeffffffff
time=15
sum=aaaaaaaeffffffff
time=12

На моем старинном Intel Core 2 Duo CPU T9400 @ 2.53GHz зажатом на 2.13 GHz, первая версия выполняется за 15 секунд, вторая за 12. Может это и не впечатляющий результат, но, как говорят, "байт там, байт здесь и оппа, в итоге сэкономили существенно"...

Конечно, это работает быстрее потому что вычисление $2x$ это просто битовый сдвиг влево, затем инкремент. Это работает быстрее, чем умножать $x$ на $x$.

На самом деле, точный результат - 0x2aaaaaaa4aaaaaaaeffffffff, но наша программа вычисляет сумму по модулю $2^{64}$ (ширина регистра). В любом случае, корректный результат нам не нужен, а нужно только удостовериться, что оба способа выдают одно и тоже.

Теперь если мы будем считать $x^3$?

В Wolfram Mathematica находим ф-цию разницы:

In[]:= f[x_] = x^3

In[]:= ExpandAll[f[x + 1] - f[x]]
Out[]= 1 + 3x + 3x^2

И снова что-то знакомое. Производная от $x^3$ это $3x^2$.

Это работает, конечно:

#!/usr/bin/env python3

for x in range(0, 15):
    print (x, x*x*x)

print ("")

a=0
for x in range(0, 15):
    print (x, a)
    a = a + 3*x*x + 1 + x*3

И опять измерим скорость на чистом Си:

#include <stdint.h>
#include <stdio.h>
#include <time.h>

int main()
{
        uint64_t _max=0x1ffffffff;
        uint64_t sum=0;
        time_t start=time(NULL);

        for (uint64_t i=0; i<_max; i++)
        {
                sum=sum+i*i*i;
        };
        printf ("sum=%lx\n", sum);
        printf ("time=%ld\n", time(NULL)-start);

        sum=0;
        start=time(NULL);
        uint64_t prev=0;
        for (uint64_t i=0; i<_max; i++)
        {
                sum=sum+prev;
                uint64_t tmp=3*i;
                prev=prev+(tmp*i + tmp + 1);
        };
        printf ("sum=%lx\n", sum);
        printf ("time=%ld\n", time(NULL)-start);
};
sum=fffffffa00000001
time=21
sum=fffffffa00000001
time=17

Немного, но всё же кое-что. Во втором случае снова 2 операции умножения, но работает немного быстрее. Наверное, потому что $3x$ вычислять быстрее чем умножать $x$ на произвольную величину - $3x = 2x + x$, а $2x$ это битовый сдвиг влево, обычно компиляторы это оптимизируют так.

А для чисел с плавающей точкой?

Но рано радоваться, всё это работает хорошо только для целых чисел. Для чисел с плавающей точкой всё может быть иначе.

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

Это можно сравнить с распиливанием доски. Вам заказчик приносит доску и просит отпилить ему такую же. Но измерить вы её почему-то не можете, рулетку у вас кто-то спёр. Вы прикладываете доску заказчика к своей доске и отпиливаете кусок. Заказчик просит еще одну, а свою, "эталонную", он уносит, ему она срочно нужна. Вы свою "копию" прикладываете к другой доске и повторяете операцию. Заказчик прибегает, спешно забирает у вас "копию", оставляя вас со "второй копией". Легко заметить, что если с таким заказчиком вы так будете отрезать доски и дальше, их длина будет постепенно увеличиваться.

Если ошибка каждый раз - 1-2 мм (еще незаметно), то через 10 досок, ошибка может быть 1-2 см. Через 100 досок ошибка может быть 1-2 дм, это уже совсем неприемлемо.

Похожее происходит у нас - мы постоянно что-то прибавляем к аккумулятору, в итоге получается очень длинная цепь сложений, а для чисел с плавающей точкой это очень плохо.

#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <time.h>

int main()
{
        uint64_t _max=500000;
        float sum1=0;
        time_t start=time(NULL);

        for (uint64_t i=0; i<_max; i++)
        {
                float _i=(float)i;
                sum1=sum1+_i*_i*_i;
        };
        printf ("sum1=%f\n", sum1);
        printf ("time=%ld\n", time(NULL)-start);

        float sum2=0;
        start=time(NULL);
        float prev=0;
        for (uint64_t i=0; i<_max; i++)
        {
                float _i=(float)i;
                sum2=sum2+prev;
                float tmp=3*_i;
                prev=prev+(tmp*_i + tmp + 1);
        };
        printf ("sum2=%f\n", sum2);
        printf ("time=%ld\n", time(NULL)-start);

        double true_error=fabs((double)sum1-(double)sum2);
        printf ("true error=%f\n", true_error);
        printf ("relative true error=%f\n", (true_error / (double)sum1));
};

Я сознательно использую float (чтобы проблемы проявились как можно раньше), а для измерения относительной ошибки - double.

sum1=15625012601280660504576.000000
time=0
sum2=15625137576170320035840.000000
time=0
true error=124974889659531264.000000
relative true error=0.000008

Ошибка в шестом знаке. Впрочем, первый результат тоже далёк от корректности. Судя по Wolfram Mathematica, верное значение: 15624937500062500000000 (ошибка в пятом знаке).

Увеличиваем _max до 1000000:

sum1=249984565122584337711104.000000
time=0
sum2=249976044312089352732672.000000
time=0
true error=8520810494984978432.000000
relative true error=0.000034

Ошибка уже ощутимая - в пятом знаке. А верное значение (Wolfram Mathematica): 249999500000250000000000 (и снова ошибка в пятом знаке).

Увеличиваем _max до 10000000:

sum1=2493939293444969611153899520.000000
time=0
sum2=2487728496076280489639411712.000000
time=0
true error=6210797368689121514487808.000000
relative true error=0.002490

Тут уже ошибка в третьем знаке. Wolfram Mathematica: 2499999500000025000000000000.

Увеличиваем _max до 100000000:

sum1=20282409603651670423947251286016.000000
time=0
sum2=10141204801825835211973625643008.000000
time=1
true error=10141204801825835211973625643008.000000
relative true error=0.500000

Тут о правильности вычислений уже говорить не приходится. Wolfram Mathematica: 24999999500000002500000000000000 (вы легко можете заметить закономерности в этих числах).

Ну хоть померяем скорость при _max=10000000000:

sum1=21267647932558653966460912964485513216.000000
time=32
sum2=166153499473114484112975882535043072.000000
time=40
true error=21101494433085539482347937081950470144.000000
relative true error=0.992188

Увы, работает медленнее.

А почему "почти" производная?

А теперь, как так получается, что вот эта вот "ф-ция разницы" выдаваемая Mathematica почти как производная, но с какими-то дополнительными коэффициентами?

Находим "ф-цию разницы" для $x^2$ постепенно уменьшая шаг с 1 до 0.001:

In[]:= f[x_] = x^2

In[]:= ExpandAll[f[x + 1] - f[x]/1]
Out[]= 1 + 2x

In[]:= ExpandAll[(f[x + 0.5] - f[x])/0.5]
Out[]= 0.5 + 2x

In[]:= ExpandAll[(f[x + 0.1] - f[x])/0.1]
Out[]= 0.1 + 2x

In[]:= ExpandAll[(f[x + 0.01] - f[x])/0.01]
Out[]= 0.01 + 2x

In[]:= ExpandAll[(f[x + 0.001] - f[x])/0.001]
Out[]= 0.001 + 2x

Видите, "коэффициент" слева постепенно уменьшается, но $2x$ остается как есть.

А для $x^3$?

In[]:= f[x_] = x^3

In[]:= ExpandAll[f[x + 1] - f[x]/1]
Out[]= 1 + 3x + 3x^2

In[]:= ExpandAll[(f[x + 0.5] - f[x])/0.5]
Out[]= 0.25 + 1.5x + 3x^2

In[]:= ExpandAll[(f[x + 0.1] - f[x])/0.1]
Out[]= 0.01 + 0.3x + 3x^2

In[]:= ExpandAll[(f[x + 0.01] - f[x])/0.01]
Out[]= 0.0001 + 0.03x + 3x^2

In[]:= ExpandAll[(f[x + 0.001] - f[x])/0.001]
Out[]= 1.*10^-6 + 0.003x + 3x^2

Тут слева теперь 2 коэффициента полинома, что постепенно уменьшаются, но $3x^2$ остается как есть.

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

Производная - это предел "ф-ции разницы" при шаге стремящемся к нулю:

In[]:= f[x_] = x^2

In[]:= Limit[(f[x + step] - f[x])/step, step -> 0]
Out[]= 2x

In[]:= D[x^2]
Out[]= 2x

In[]:= f[x_] = x^3

In[]:= Limit[(f[x + step] - f[x])/step, step -> 0]
Out[]= 3x^2

In[]:= D[x^3]
Out[]= 3x^2

Предел, потому что шаг никогда не будет нулем, он только стремится к нему.


Please drop me email about bug(s) and/or suggestion(s): blog@yurichev.com.