Проверка что число - куб

[Dmitriy40]

wrest в сообщении #1331292 писал(а):

и найти не могу

Ну, реализация ispower начинается с arith1.c:1115 (строка 1115, все файлы из каталога \src\basemath\), быстро уходит в Z_ispowerall (arith1.c:1034), которая для малых степеней уходит в is_357_power (ifactor1.c:1908), которая в строках после 1934 с помощью uis_357_powermod (ifactor1.c:1835, кстати познавательно взглянуть как сделано) для многих простых вариантов не кубов вылетает уже здесь, а в строке 1954 извлекает корень с заданной погрешностью и проверяет полученный результат возведением обратно в степень и сравнением. Корень извлекается функцией sqrtnr, вызывающей sqrtnr_abs из trans1.c:1768, для которой добрые дяденьки написали коммент "cubic Newton iteration, |a|^(1/n)", что позволяет на этом остановиться и разбирательство засчитать успешным. Ну а кто хочет может сделать больше.

[mihaild]

Набросал небольшой бенчмарк: https://github.com/mihaild/dxdy/tree/ma ... _cube_root
Второй миллиард проверяется за 12 секунд бинпоиском и за 50-60 секунд более интеллектуальными способами. Использование проверки по модулю ускоряет вектор примерно в 10 раз, интеллектуальные методы - в 30.

[kotenok gav]

А почему не извлечь кубический корень и проверить его на целость? Я всегда так делаю.

-- 09 авг 2018, 10:51 --

mihaild в сообщении #1331307 писал(а):

Набросал небольшой бенчмарк: https://github.com/mihaild/dxdy/tree/ma ... _cube_root

А можно на питоне тоже?

[Dmitriy40]

mihaild
Как думаете, почему CheckIcbrt (да и прочие) в порядке 575757,88236 выполняется дольше чем в обратном? Итоговое количество проходов сквозь оба сита одинаково, но второе сито в таком порядке должно выполняться меньшее количество раз, чем второе в обратном.

kotenok gav в сообщении #1331311 писал(а):

А почему не извлечь кубический корень и проверить его на целость? Я всегда так делаю.

До-о-о-олго! А тут идёт борьба за скорость. Потому и в бенчмарке на питоне смысла мало, и так понятно что медленно. Плюс возможные ошибки округления при вычислении корня. Тут в коде почти такой метод есть, CheckPow, вычисляет степень 1/3 и проверяет соседние целые числа на точное равенство после возведения в куб (не совсем именно так, но смысл примерно таков).

[FomaNeverov]

Вычислить куб числа не так просто и быстро. Но возвести число в 3-ю степень (три раза умножить само на себя) очевидно намного быстрее. Предварительно вычислить таблицу кубов чисел интересующего диапазона чисел и сравнивать уже с этой таблицей.

Пример, кубы 3-х чисел:

(123456789) 1881676371789154860897069
(123456790) 1881676417513891481839000
(123456791) 1881676463238628843521671

[wrest]

Dmitriy40 в сообщении #1331305 писал(а):

которая для малых степеней уходит в is_357_power (ifactor1.c:1908), которая в строках после 1934 с помощью uis_357_powermod (ifactor1.c:1835, кстати познавательно взглянуть как сделано)

Спасибо! Ну что тут скажешь? Все украдено до нас!
Проверка по модулю для куба в функции ispower PARI/GP тоже проводится, потому и ускорения у нас не вышло.
То есть покашта, для использования с PARI/GP, встроенная ispower выходит быстрее чем поиск поиск в векторе остатков и даже чем поиск в битовом массиве предвычисленном для всех возможных кубов?

-- 09.08.2018, 10:25 --

(Dmitry40)

Dmitriy40 в сообщении #1331248 писал(а):

Что тут непонятного?

Ну вы так уж не бухтите, я С++ плохо понимать, в школе не учить, чем отличается ++с от с++ сходу не знать :) Спасибо, что растолковываете.

[Dmitriy40]

FomaNeverov в сообщении #1331319 писал(а):

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

Это уже предложено. И реализовано в бенчмарке как метод CheckVector/CheckBinary.

wrest в сообщении #1331329 писал(а):

То есть покашта, для использования с PARI/GP, встроенная ispower выходит быстрее чем поиск поиск в векторе остатков и даже чем поиск в битовом массиве предвычисленном для всех возможных кубов?

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

[wrest]

Dmitriy40 в сообщении #1331300 писал(а):

n=15561 / 315 : 10780ms, 3994ms.

Dmitriy40 в сообщении #1331300 писал(а):

*** last result computed in 57,002 ms.

Ага, то есть скомпилированная программка С считает быстрее интерпретируемой PARI/GP...

Встаёт вопрос о GP2C: как им ользоваться и будет ли для долгих вычислений типа обсуждаемого поиска решений головоломки быстрее и насколько, и будет ли быстрее чем просто на С.

[Dmitriy40]

(mihaild)

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

   jnc   .next
   mov   EBX,EDX
   ;cmovc   EBX,EDX
.next:

как по идее и рекомендуют делать, и получаю общее замедление почти в два раза! Хотя команда cmov по доке занимает лишь два такта из почти десятка на весь цикл, и однако ж. Весь остальной код без изменений. Второй миллиард с командой cmov считается 25с, с условным переходом 15с. Парадокс.

wrest в сообщении #1331339 писал(а):

Ага, то есть скомпилированная программка С считает быстрее интерпретируемой PARI/GP...

Причина не только в этом, я же не заморачивался с кучей проверок (даже возможное переполнение не проверяю). Но сравнить чистую скорость ispower и C кода трудно, это надо пересобирать PARI с встраиванием прямо в ispower функции точного замера времени (ну или выдрать весь кусок кода ispower в свою прогу). Да и не нужно этого, как раз сравнение как было и самое близкое к реальности.

wrest в сообщении #1331339 писал(а):

и будет ли быстрее чем просто на С

Думаю более эффективно будет векторизовать (с моим любимым AVX2) проверку и распараллелить вычисления на все ядра.

[wrest]

Dmitriy40 в сообщении #1331342 писал(а):

Думаю более эффективно будет векторизовать (с моим любимым AVX2) проверку и распараллелить вычисления на все ядра.

Да, это тоже вопрос: PARI/GP на виндовс 7 отъедает 25% CPU (CPU 4-поточное).
Тама есть параллельный for который parfor но вот я попробовал, работает аномально:

(12:05) gp > parfor(i=10^15,10^15+10^8,ispower(i,3))
(12:05) gp > ##
  ***   last result computed in 18,502 ms.
(12:05) gp > for(i=10^15,10^15+10^8,ispower(i,3))
(12:06) gp > ##
  ***   last result computed in 8,238 ms.

-- 09.08.2018, 12:13 --

Dmitriy40 в сообщении #1331342 писал(а):

Причина не только в этом, я же не заморачивался с кучей проверок (даже возможное переполнение не проверяю). Но сравнить чистую скорость ispower и C кода трудно, это надо пересобирать PARI с встраиванием прямо в ispower функции точного замера времени (ну или выдрать весь кусок кода ispower в свою прогу). Да и не нужно этого, как раз сравнение как было и самое близкое к реальности.

Да не, я о другом: вот задаст Ktina очередную загадку. После отладки проверяющей функции выяснится что считать надо долго. Вопрос: лучше откомпилировать функцию через gp2c или считать в интерпретаторе GP или лучше вообще для длинных расчетов писать на c\c++ (что вообще мне кажется странным, учитывая что PARI/GP сделано специально для расчетов и само написано на с причем как мы видели, с помороченными оптимизациями)?

И главное: почему у mihaild проверка головоломки с 400 по 500 идет за одну секунду? Компьютер быстрее? Программа лучшее?

[Dmitriy40]

wrest в сообщении #1331347 писал(а):

И главное: почему у mihaild проверка головоломки с 400 по 500 идет за одну секунду?

А давайте грубо прикинем скорость.
Забьём на уменьшение количества проверок при MAX_Z меньше 500 и возьмём для 500. X и y будут перебираться до 2364. Для каждой пары надо проверять кубический корень. Итого $100\cdot2364\cdot2364/2\approx280\cdot10^6$. При реализации с предпросмотром по модулю 575757 benchmark выше выдаёт 1.4с на миллиард чисел или 0.4с на 280млн. Вот она, заветная цифра. Причём завышенная на четверть (реальное количество извлекаемых корней 203млн). Ну правда плюс 0.1-0.2с на инициализацию массивов.
Если скорость PARI ispower взять за 57с на миллиард, то 203млн посчитаются за 12с (проверил, да, столько). Плюс накладные расходы на возведения в степень и интерпретацию. Цикл перебора z,x,y 203млн раз, без извлечения корня, у меня выполняется 35с, плюс 12с на корни - суммарно должно быть 47с. Реально же 52с.
У Вас можете сами пересчитать по скорости ispower.

wrest в сообщении #1331347 писал(а):

учитывая что PARI/GP сделано специально для расчетов и само написано на с причем как мы видели, с помороченными оптимизациями

Вот только слово в полном названии "PARI/GP Calculator" очень жирно намекает на ориентацию на простоту вычислений, а не их скорость.

wrest в сообщении #1331347 писал(а):

Тама есть параллельный for который parfor но вот я попробовал, работает аномально:

А у Вас в 4-й строке информации о версии threading engine: single жирным тоже это написано? Что ж тогда хотите то, параллельный код в один поток плюс накладные расходы на синхронизацию потоков (лень мне искать обходятся ли они реально при исполнении в один поток) и дают видимый вами проигрыш скорости. Вполне ожидаемо, хотя и неприятно.
Вроде бы для корректной работы в несколько потоков должна быть подключена (установлена?) библиотека MPI, но как - не знаю.

[wrest]

Dmitriy40 в сообщении #1331380 писал(а):

А у Вас в 4-й строке информации о версии threading engine: single жирным тоже это написано?

А... single написано, да.

-- 09.08.2018, 14:54 --

Dmitriy40 в сообщении #1331380 писал(а):

Цикл перебора z,x,y 203млн раз, без извлечения корня, у меня выполняется 35с, плюс 12с на корни - суммарно должно быть 47с. Реально же 52с.

Ага, значит собственно определение "куб ли число" там не основная затрата по времени.

[Dmitriy40]

Значит не основная.
Код на всякий случай:

? for(z=400,500,z5=z^5;max_xy=sqrtnint(z5,4);t=0;for(x=1,max_xy,x4=x^4;for(y=x,max_xy,y4=y^4;if(x4+y4>=z5,break);if(ispower(z5-x4-y4,3,&t),printf("x=%d y=%d z=%d t=%d\n",x,y,z,t)))))
x=1360 y=1632 z=408 t=9248
x=1458 y=1458 z=486 t=26244
  ***   last result computed in 1min, 5,428 ms.
? for(z=400,500,z5=z^5;max_xy=sqrtnint(z5,4);t=0;for(x=1,max_xy,x4=x^4;for(y=x,max_xy,y4=y^4;if(x4+y4>=z5,break);if(ispower(z5-x4-y4,3),ispower(z5-x4-y4,3,&t);printf("x=%d y=%d z=%d t=%d\n",x,y,z,t)))))
  ***   last result computed in 52,229 ms.
? for(z=400,500,z5=z^5;max_xy=sqrtnint(z5,4);t=0;for(x=1,max_xy,x4=x^4;for(y=x,max_xy,t3=z5-x4-y^4;if(t3<=0,break);if(ispower(t3,3),ispower(t3,3,&t);printf("x=%d y=%d z=%d t=%d\n",x,y,z,t)))))
  ***   last result computed in 47,830 ms.
? for(z=400,500,z5=z^5;max_xy=sqrtnint(z5,4);t=0;for(x=1,max_xy,x4=x^4;z5x4=z5-x4;for(y=x,max_xy,t3=z5x4-y^4;if(t3<=0,break);if(ispower(t3,3),ispower(t3,3,&t);printf("x=%d y=%d z=%d t=%d\n",x,y,z,t)))))
  ***   last result computed in 41,591 ms.

Решения везде одинаковые, не показываю для экономии места.
Зацените разницу между первым и вторым, уже не 65с, а 52с, а ведь всего лишь убрал адрес t в ispower. А потом и ещё ускорилось, хотя казалось бы, лишь сложения/вычитания выносил за цикл ... Разница между третьим и четвёртым вариантами ровно в одной операции вычитания, повторяемой 203млн раз и занимающей 6с или 30нс/операцию (больше 100 тактов, мрак!).
И в итоге ускорилось в полтора раза! Интерпретаторы дело такое, от любого чиха колбасит скорость.

[wrest]

Dmitriy40 в сообщении #1331409 писал(а):

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

Ну в общем, собрал я из исходников на планшете с андроидом gp2c
Откомпилировал ваш последний текст при помощи gp2c, ускорение скомпилированной функции иде-то в 3 с небольшим раза.
Я правда запускал со 100 до 200 в виду хилости планшета. Интерпретация: 54,6 секунд, компилированное: 16,2 секунды.

-- 09.08.2018, 17:08 --

Dmitriy40
Вот собсно как оно в gp2c преобразовалось из

ktest()=for(z=400,500,z5=z^5;max_xy=sqrtnint(z5,4);t=0;for(x=1,max_xy,x4=x^4;z5x4=z5-x4;for(y=x,max_xy,t3=z5x4-y^4;if(t3<=0,break);if(ispower(t3,3),ispower(t3,3,&t);printf("x=%d y=%d z=%d t=%d\n",x,y,z,t)))))

В текст C:

Используется синтаксис C

/*-*- compile-command: "/data/data/com.termux/files/usr/bin/gcc -c -o ktest.gp.o -O3 -Wall -fno-strict
-aliasing -fomit-frame-pointer -fPIC -I\"/data/data/com.termux/files/usr/include\" ktest.gp.c && /data
/data/com.termux/files/usr/bin/gcc -o ktest.gp.so -shared -O3 -Wall -fno-strict-aliasing -fomit-frame-
pointer -fPIC -Wl,-shared ktest.gp.o -lc -lm -L/data/data/com.termux/files/usr/lib -lpari"; -*-*/
#include <pari/pari.h>
/*
GP;install("init_ktest","v","init_ktest","./ktest.gp.so");
GP;install("ktest","v","ktest","./ktest.gp.so");
*/
void init_ktest(void);
void ktest(void);
/*End of prototype*/

void
init_ktest(void)          /* void */
{
  pari_sp ltop = avma;
  avma = ltop;
  return;
}

void
ktest(void)       /* void */
{
  pari_sp ltop = avma;
  GEN z5 = pol_x(fetch_user_var("z5")), max_xy = pol_x(fetch_user_var("max_xy")), t = pol_x(fetch_user
_var("t")), x4 = pol_x(fetch_user_var("x4")), z5x4 = pol_x(fetch_user_var("z5x4")), t3 = pol_x(fetch_u
ser_var("t3"));
  {
    pari_sp btop = avma;
    long z;
    for (z = 100; z <= 200; ++z)
    {
      z5 = powiu(stoi(z), 5);
      max_xy = sqrtnint(z5, 4);
      t = gen_0;
{
        pari_sp btop = avma;
        GEN x = gen_0;
        for (x = gen_1; gcmp(x, max_xy) <= 0; x = gaddgs(x, 1))
        {
          x4 = gpowgs(x, 4);
          z5x4 = gsub(z5, x4);
          {
            pari_sp btop = avma;
            GEN y = gen_0;
            for (y = x; gcmp(y, max_xy) <= 0; y = gaddgs(y, 1))
            {
              t3 = gsub(z5x4, gpowgs(y, 4));
              if (gcmpgs(t3, 0) <= 0)
                break;
              if (ispower(t3, stoi(3), NULL))
              {
                ispower(t3, stoi(3), &t);
                printf0("x=%d y=%d z=%d t=%d\n", mkvec4(x, y, stoi(z), t));
              }
              if (gc_needed(btop, 1))
                gerepileall(btop, 3, &t, &y, &t3);
            }
          }
          if (gc_needed(btop, 1))
            gerepileall(btop, 5, &t, &x, &x4, &z5x4, &t3);
        }
      }
      if (gc_needed(btop, 1))
        gerepileall(btop, 6, &z5, &max_xy, &t, &x4, &z5x4, &t3);
    }
  }
  avma = ltop;
  return;
}