LINUX.ORG.RU

Продолжаем тему распределённого lc-генератора.

 , , ,


1

3

Кому лень читать - надо сделать(подсказать куда копать) x*3877 mod 139968, где x = 0(а может и от одного)...139967 на 32битных числах(а в идеале на float), быстро. Быстро - это максимум штук десять умножений/сложений.

Всем кто помог в прошлой теме - спасибо. В конечном итоге вышло это: http://pastebin.com/rnLwXiUv

//g++ main.cpp -Ofast -march=native -flax-vector-conversions -std=gnu++1y -pthread//пока только avx intel sb+
(0.914s)13.53tpc
(0.073s)1.08tpc

Это конечно хорошо(на самом деле оно тупо упирается в память и реально - оно раз в 20быстрее).

Но этого мало, проблема в том, что оно на даблах из-за «быстрого» «деления» и из-за того, что интел забил на целочисленную арифметику. Ну из-за пту.

Т.е. даблы только из-за:

double mod_139968(double n) {
  return n + (floor(n * 1. / 139968) * -139968);
}
//n * 1. / 139968 - конкретно этого

Я пытался смочь сделать это на 32битных интах - у меня не вышло, ибо пту. Прошу помощи.

dikiy GPFault

Ответ на: комментарий от f1u77y

Логарифм от чего? В любом случае слишком долго. Но за идею спасибо. Скорее всего для частного случая там будет не логарифм.

registrant27492
() автор топика
Ответ на: комментарий от registrant27492

Логарифм от чего?

от степени. смотри мой последний коммент в предыдущем треде. Я еще не до конца продумал, но профит может быть даже по сравнению с Эйлером. Если этот метод еще до кучи взять.


Кстати, я бы поостерегся использовать float. Там ошибки же могут быть. Используй какую-нибудь либу длинной арифметики. Я нашел GNU MP library (gmplib).

dikiy ★★☆☆☆
()
Ответ на: комментарий от f1u77y

Хотя я тут подумал и реально. (x×n) mod 139968 = (((x × (n / 2) ) mod 139968) * 2) mod 139968 - надо просто разбить первый множитель на n, чтобы результат не переполнялся от умножения во флоате.

Спасибо.

registrant27492
() автор топика

Вот зачем все эти извращения, когда

139968 = 2^6 * 3^7
Вспоминаем малую теорему Ферма.

iVS ★★★★★
()

Глянь, как реализован Exp в Go. → https://golang.org/src/math/big/int.go?s=9982:10018#L398

Пример: http://play.golang.org/p/HpoEEbUM9h

BenchmarkExp-4	10000000	       156 ns/op

А ещё глянь тут: BN_mod_exp

BN_mod_exp() computes a to the p-th power modulo m («r=a^p % m»). This function uses less time and space than BN_exp().

В общем я всё это к тому, что не надо изобретать велосипед, а стоит воспользоваться готовыми решениями.

beastie ★★★★★
()
Последнее исправление: beastie (всего исправлений: 3)
Ответ на: комментарий от dikiy

Кстати, я бы поостерегся использовать float. Там ошибки же могут быть.

Там флоаты используются как целые числа, а как флоаты они используютсят олько в модуле, чтобы можно было быстро умножить на обратное число. Ну и само-собою потому, что они быстрее целых чисел.

Используй какую-нибудь либу длинной арифметики. Я нашел GNU MP library (gmplib).

Нельзя использовать никакую длинную арифметику, ибо теряется вся суть этой штуки.

Я уже писал в предыдущей теме, что в пистоне выполняется (((3877 ** 100000) - 1) / 3876) % 139968 за время, за которое эта реализация гинерирует 250000000 целых чисел по формуле:

(3877 ^ n * 42 + ((3877 ^ n - 1) / (3877 - 1)) * 29573) (mod m)

Мне надо 3числа(для 4тредов), даже если поверить в то, что gmp быстрее питона(разве там не гмп?) в 10раз - это просто помножит на ноль весь смысл задачи.

Я еще не до конца продумал, но профит может быть даже по сравнению с Эйлером. Если этот метод еще до кучи взять.

Это может помочь тут только в плане, если кешировать результаты и не пересчитывать заново. Я подумаю над этим.

Просто пока по условию задачи так читерить нельзя.

registrant27492
() автор топика
Ответ на: комментарий от f1u77y

хотя стоп, а что не так с тем, чтобы просто посчитать в 32-битных, оно жи не переполнится

Нельзя делить.

Вот надо как-то посчитать модуль без деления не перполняя 32битные числа, а лучше флоат.

registrant27492
() автор топика
Ответ на: комментарий от registrant27492

Нельзя использовать никакую длинную арифметику, ибо теряется вся суть этой штуки.

Если решение на длинной арифметике обходит твоё на порядки, то почему бы и нет? Вон, решение на Go даёт тебе ~6500000 операций в секунду. С BigNum может быть будет ещё быстрее.

beastie ★★★★★
()
Ответ на: комментарий от registrant27492

Я уже писал в предыдущей теме, что в пистоне выполняется (((3877 ** 100000) - 1) / 3876) % 139968 за время, за которое эта реализация гинерирует 250000000 целых чисел по формуле:

(3877 ^ n * 42 + ((3877 ^ n - 1) / (3877 - 1)) * 29573) (mod m)

так сделай то же самое по этой формуле, только замени везде n на остаток от деления n % phi(139968).

dikiy ★★☆☆☆
()
Ответ на: комментарий от beastie

Если решение на длинной арифметике обходит твоё на порядки, то почему бы и нет? Вон, решение на Go даёт тебе ~6500000 операций в секунду. С BigNum может быть будет ещё быстрее.

На порядки. Как это мило.

Давай ещё раз. Как бы тебе попроще объяснить. Это и то, что ты мне предлагаешь - находятся в разных весовых категориях. Они не сравнимы.

Для начала - я просил умножение под модулем, а не степень - это раз.

А два - моя реализация даёт 5000000000 операций в секунду(рандомных чисел) из которых на каждую приходится 3 взятия модуля, т.е. у меня 15000000000 модулей в секунду в пересчёте на ведро.

Я думаю мне не надо объяснять что из себя будет представлять релизация на BigNum, но ты можешь попытаться.

registrant27492
() автор топика
Ответ на: комментарий от dikiy

так сделай то же самое по этой формуле, только замени везде n на остаток от деления n % phi(139968).

Это ничего не даст, там вычисление последовательное.

Это даст только, если кешировать вычисления, но:

http://benchmarksgame.alioth.debian.org/u64q/fasta-description.html#fasta

don't cache the random number sequence

Нельзя по условию задачи.

registrant27492
() автор топика
Ответ на: комментарий от beastie

Go

я думал, хоть в тему про алгоритмы его не получится всунуть, ан нет.

давай лучше на rust

f1u77y ★★★★
()
Ответ на: комментарий от registrant27492

А нет, я ошибся - у меня 2 взятия модуля на операцию и один на каждые 16. Но это не особо что-то меняет.

registrant27492
() автор топика

На float (суть 23-битные числа для данной задачи) можно переписать, если разбить x на старших 9 и младших 9 бит, после чего оно будет умещаться во float и после умножения на 3877.

Только операций получается дофига - 6 умножений, 3 сложения и 2 битовые - их может быть проблематично записать в sse...

У меня получилось следующее:

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

float mod_by_139968(float x)
{
	const float float_divider = 1./139968.;
	return x - 139968.f * floorf(((float)x) * float_divider);
}

float floats_nodivision(float x)
{
	float high_part = ((int)x) & ~0x1ff;
	float low_part = ((int)x) & 0x1ff;
	float high_mod = mod_by_139968(3877.f * high_part);
	return mod_by_139968(3877.f * low_part + high_mod);
}

int basic(int x)
{
	return  (x*3877)%139968;
}

int main()
{
	int x = 0;
	for(;x < 139968;++x)
	{
		if (basic(x) != floats_nodivision(x))
		{
			break;
		}

	}
	printf("sizeof(float) = %d; exiting with x = %d", sizeof(float), x);
	return 0;
}

Правда некоторая странность: от -O0 до -O3 работает, а с -Ofast выдаёт неверны результат на x = 4329, почему не разбирался, скорость не проверял.

GPFault ★★
()

А ещё можно воспользоваться хорошей кармой данных чисел: (3877 * 2^8) mod 139968 = 12736 получается довольно маленьким, что позволяет сложить в один float первые_8_бит*3877 + остальные_биты*12736

По количеству операций это выглядит более перспективно - битовые операции более простые «разделение на байты» и 4 умножения, 2 сложения

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

float mod_by_139968(float x)
{
	const float float_divider = 1./139968.;
	return x - 139968.f * floorf(((float)x) * float_divider);
}

const float low_mul = 3877;
const float high_mul = (3877 * 1<<8) % 139968;

float floats_nodivision_for_special_numbers(float x)
{
	float low_byte = ((int)x) & 0xff;
	float high_part = (((int)x) >> 8);
	return mod_by_139968(low_mul * low_byte + high_mul * high_part);
}

int basic(int x)
{
	return  (x*3877)%139968;
}

int main()
{
	assert(low_mul * 0xff + high_mul * (139967 >> 8) < 1 << 23);
	int x = 0;
	for(;x < 139968;++x)
	{
		if (basic(x) != floats_nodivision_for_special_numbers(x))
		{
			break;
		}

	}
	printf("sizeof(float) = %d; exiting with x = %d", sizeof(float), x);
	return 0;
}

это работает и с -Ofast.

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

GPFault ★★
()
Ответ на: комментарий от GPFault

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

битовые операции более простые

Сдвиг дорогой.

Грубо. сложение + 2 умножения + битоп + сдвиг. Это такт с 2-х fma + такт со сдвига. В районе 2.5 тактов, если учитывать фантомный 3-й модуль * 2 модуля полных и того 5тактов. Т.е. это выгодно. Там дальше можно выкинуть одну конвертацию.

Но с флоатом ещё одна проблема:

(3877 ^ n * 42 + ((3877 ^ n - 1) / (3877 - 1)) * 29573) (mod m)
//139968*4*29573 уже не сможет

Переходить на int то же не вариант, ибо он медленнее.

registrant27492
() автор топика
Ответ на: комментарий от GPFault

Кстати, если речь идёт о работе на x86

Да.

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

Табличкой? 139968 флотов? Будет раз в 20медленней текущей реализации.

registrant27492
() автор топика
Ответ на: комментарий от registrant27492

К сожалению я не очень в курсе актуальных сложностей операций на x86-платформе - на моих дешёвых ARMах наоборот всё на int перевожу)

По идее непосредственно сдвиг для float не нужен - он может быть заменён AND, зануляющим младшие 8 бит. Потом множитель high_mul можно взять в 256 раз меньше (без потери точности множителя, так как он во float).

Или And тоже дорогой?

GPFault ★★
()
Ответ на: комментарий от GPFault

К сожалению я не очень в курсе актуальных сложностей операций на x86-платформе - на моих дешёвых ARMах наоборот всё на int перевожу)

Ну сейчас на х86 трупут fma в 2раза выше, чем у целочисленного умножения/сдвига и такое же, как у целочисленного сложения.

Т.е. да, сейчас на х86 сделать 2сложения+2умножения стоит столько же, сколько сделать сдвиг/умножить.

Это исправили(частично) только в самом последнем штеуде.

По идее непосредственно сдвиг для float не нужен - он может быть заменён AND, зануляющим младшие 8 бит. Потом множитель high_mul можно взять в 256 раз меньше (без потери точности множителя, так как он во float).

Спасибо.

Или And тоже дорогой?

Да, флотовский and стоит столько же, сколько и сдвиг, но это быстрее, ибо сдвига для флоатов там нет и приходится кастовать, на это там есть пенальти.

И флоатовый and то же исправили.

registrant27492
() автор топика
Ответ на: комментарий от Deleted

А что такого? Я стараюсь ради пацанов. Пацан в соседнем треде расстроился, что сишка сливает расту. Не беда - я исправлю ситуацию. Но мне не интересно как пацанам на 10% - мне интересно в 20раз. Такие дела. Пойду что-ли хоть пока 10-15раз для разогрева запощу. Даже причёсывать не буду. Профиты полутают крестовики ибо прототип вырос из крестовой версии.

Надо меня направлять в нужное русло - чтоб я пилил. И пацаны молодцы - мне помогают. Я щастлив ^_^

registrant27492
() автор топика
Ответ на: комментарий от Deleted

Да и выше я пруфцанул реализацию:

uint32_t gen_random(void) {
  const int IM = 139968, IA = 3877, IC = 29573;
  static int last = 42;
  last = (last * IA + IC) % IM;
  return last;
}

Которая работает от 20раз быстрее и полностью горизонтально масштабируема(пацаны же это любят). И меня никто не похвалил, а я целую неделю карпел. Как так можно?

registrant27492
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.