Простой ГСП на Си Опции

[Herz]

Предположим, требуется моргать светодиодом с псевдослучайным интервалом от 2 до 30 секунд.
Посоветуйте простую функцию на Си для генерации целых чисел в таком малом диапазоне. Пусть от 1 до 255 максимум.
Функция rand() из библиотеки HT-PICC тяжеловата как-то...

[thermit]

x(n)=(1664525*x(n-1)+32767) mod 2^32 и берем старшие N бит.

[Guest_TSerg_*]

Любой RNG на основе сдвигового регистра 8 р, но лучше - 16 р

u8 rnd8(void)
{
static u16 seed = 0;

seed = (seed << 7) - seed + 251;

return (u8)(seed + (seed>>8));
}

[Herz]

Спасибо, коллеги за скорую помощь!

Любой RNG на основе сдвигового регистра 8 р, но лучше - 16 р

u8 rnd8(void)
{
static u16 seed = 0;

seed = (seed << 7) - seed + 251;

return (u8)(seed + (seed>>8));
}

Я только не пойму, если в теле функции переменная seed каждый раз обнуляется, то какой смысл в static?

[aaarrr]

Я только не пойму, если в теле функции переменная seed каждый раз обнуляется, то какой смысл в static?

Она не обнуляется, 0 присваивается только при инициализации программы, т.е. до первого вызова.

[Herz]

Она не обнуляется, 0 присваивается только при инициализации программы, т.е. до первого вызова.

Отлично, спасибо! А можно ли немного усложнить задачу и передавать аргументами в функцию желаемый диапазон значений?

[SM]

интервал определен как [offset, offset+max)

два варианта:

1) - два параметра, max и offset

return ((rand_value) % max) + offset;

это если % достаточно легковесный.

2) - три параметра, mask, max и offset, mask ближайшее значение типа 2^n-1, большее или равное max

temp = rand_value & mask;
if (temp >= max) temp -= max;
return temp + offset;

это без операции деления, вероятно менее затратный вычислительно.

rand_value - результат вышеприведенной ф-ции.

[Herz]

Спасибо, SM!
Но, по правде говоря, ничего не понял.
Вас на затруднит немного разжевать?
Вот если у меня есть две переменные:

Код

unsigned char min = 2;
unsigned char max = 30;

и я вызываю функцию

unsigned char rand(a, b ) :

result = rand(min, max);

желая получить число в указанном диапазоне. Как должно выглядеть её тело?
Так?

Код

unsigned char rand(a, b)
{
static unsigned int rand_value
return (unsigned char)((rand_value) % b) + a;
}

[_Pasha]

Код

u8 rnd8(u8 min, u8 max)
{
   static u16 seed = 0;
   u8 tmp;

   seed =(seed << 7) - seed + 251;
   tmp = (u8)(seed+ (seed>>8));

   while( tmp > max) tmp-=max;// для данного случая это будет всегда выигрышнее, чем  целочисленное деление

return tmp + min;
}

[SM]

Вас на затруднит немного разжевать?

желая получить число в указанном диапазоне. Как должно выглядеть её тело?

Если за базу взять ф-цию, предложенную TSerg, то вот так:

медленный код (он притащит библиотечную ф-цию деления с остатком):

суть - тупое вычисление в лоб остатка от деления ограничиваемого числа на max.

Код

u8 rnd8(u8 max, u8 offset)
{
static u16 seed = 0;
u8 temp;

seed = (seed << 7) - seed + 251;

temp =  (u8)(seed + (seed>>8));

return (temp % max) + offset;

}

самый быстрый код, но требующий маски, которая задается параметром ф-ции, и должна быть числом 2^n-1, ближайшим, больше либо равным "max".

Суть кода - операция лог. И с маской, соответствующей требованию ближайшего, большего либо равного 2^n-1, к max, ограничивает значение неким числом 2^n-1, меньшим max*2, затем в одну итерацию вычитания по условию, число ограничивается уже max.

Если требуется вычисление маски во время выполнения, то есть "max" задается не константами, и требуется максимальная скорость, то можно добавить получение маски по таблице (массиву констант на 256 значений - тут уж хотите быстро, раскошеливайтесь на память для таблиц, другие способы вычисления такой маски на АВР в любом случае будут в цикле до 8 итераций в худшем случае)...

Код

u8 rnd8(u8 mask, u8 max, u8 offset)
{
static u16 seed = 0;
u8 temp;

seed = (seed << 7) - seed + 251;

temp =  ((u8)(seed + (seed>>8))) & mask;
if (temp >= max) temp -= max;

return temp + offset;

}

вот еще вариант, с операцией вычисления остатка от деления, реализованной в теле ф-ции без использования библиотечных ф-ций, средне производительный.

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

Код

u8 rnd8(u8 max, u8 offset)
{
static u16 seed = 0;

u8 temp, i;
u16 sub;

seed = (seed << 7) - seed + 251;

temp =  (u8)(seed + (seed>>8));

i=8; sub = max << 8;

do {
  sub >>= 1;
  if (sub >= temp) temp -= sub;
} while (--i);

return temp + offset;

}

для данного случая это будет всегда выигрышнее, чем целочисленное деление

Да щасссс... целочисленное деление всегда 8 итераций (см. мой третий вариант). А у вас.... прикиньте, например, для частного случая max = 3 и tmp=254.

[_Pasha]

Да щасссс... целочисленное деление всегда 8 итераций (см. мой третий вариант). А у вас.... прикиньте, например, для частного случая max = 3 и tmp=254.

А ну да, маску надо задавать.
Ага, и нафига там min в параметрах? Для красоты


Сообщение отредактировал _Pasha - Jan 18 2014, 08:06

[Herz]

Спасибо большое, есть что покурить.

Ага, и нафига там min в параметрах? Для красоты

Ну, это пустяки. Если функция возвращает число в диапазоне от 0, то можно предварительно max задать уменьшенным на min, а потом к результату этот min добавить, как у Вас.

[_Pasha]

Ну, это пустяки. Если функция возвращает число в диапазоне от 0, то можно предварительно max задать уменьшенным на min, а потом к результату этот min добавить, как у Вас.

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

[SM]

что min должен быть передан в параметрах

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

[Herz]

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

Замечание правильное. В реальности такой необходимости нет пока. А условие я формулировал так для наглядности, чтобы самому было легче разобраться.

[Herz]

Странно, но вот у меня такая функция:

Код

unsigned char rnd8(unsigned char max)
{
static unsigned int seed = 0;
unsigned char temp, i;
unsigned int    sub;

seed = (seed << 7) - seed + 251;
temp =  (unsigned char)(seed + (seed>>8));
i = 8; sub = max << 8;

do     {
      sub >>= 1;
      if (sub >= temp) temp -= sub;
    } while (--i);

return temp;
}

(минимум опустил пока для простоты)

при вызове:

Код

NUM = rnd8(MAX);

где

Код

#define        MAX        30

работает неверно. Диапазон не соблюдает. Более того, старший бит всегда 1.

[SM]

Я допустил описку, сорри... Писал прямо в форум...

if (sub >= temp) должно быть наоборот, if (temp >= sub)

[Herz]

Семён Семёныч!
Да, теперь работает отлично. Ещё раз большое спасибо.
Кстати, temp не приводится к типу sub при temp -= sub; ?

[SM]

Кстати, temp не приводится к типу sub при temp -= sub; ?

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

[Herz]

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

Нет, я имею в виду, не возвращает ли функция int в этом случае вместо char после приведения? Или типом функции это задаётся строго?

[SM]

Нет, я имею в виду, не возвращает ли функция int в этом случае вместо char после приведения? Или типом функции это задаётся строго?

В контексте Вашего вопроса, приведения в "temp -= sub", это задается типом переменной temp, при присваивании в 8-битную переменную автоматически все лишнее обрезается, и оптимизатор, по идее, должен сам определить это, и использовать 8-битные операции для всего выражения (которое по ANSI должно было считаться в 16 битах, а потом обрезаться).

Касаемо значения, возвращаемого функцией, да, оно строго задано типом в описании функции, и не может быть никаким, кроме заданного.

UPD:
PS
по логике, кстати, в том выражении надо не temp приводить к типу sub, а наоборот, учитывая, что по условию там temp >= sub, что означает, что реальная разрядность числа в sub меньше либо равна 8. Т.е. для наглядности и для подсказки оптимизатору можно написать temp -= (u8) sub;

[Herz]

ОК, я понял. То есть и у TSerg смысла в явном приведении в строке:

Код

return (u8)(seed + (seed>>8));

не было?

[SM]

То есть и у TSerg смысла в явном приведении в строке:

Да. Только наглядность. Ну и, возможно, уничтожение warning-а компилятора о том, что теряется разрядность, если этот warning включен по умолчанию.

[Fat Robot]

Не нужно вам никакого деления и остатка от деления.

1. делаете генератор в полном диапазоне int16 или int8. с помощью LFSR. Вам написали.

2. рассматриваете выход генератора как fractional число Q1.15, если выход int16, или Q1.7, если выход int8. Таким образом у вас есть RNG с fractional выходом в диапазоне [-1; 1)

3. умножаем (fractional), прибавляем, сдвигаем, приводим типы (сделаете в удобной последовательности), чтобы получить нужный диапазон

можно обнулять старший бит и тогда на выходе RNG получать [0; 1). Дальнейшая арифметика в этом случае видится немного попроще.
Будет, правда, слегка перекошенное распределение (для ненулевого числа x != 0 будет вероятность выпадения p(x); если же x = 0, то p(0) ~= 1.5*p(x) ), но для мигалки это не принципиально, мне кажется.

Итого: сдвиги, суммы, одно умножение

[SM]

одно умножение

так вот в нем и зло. если нету ни аппаратного умножителя, ни аппаратного делителя, то и на умножение получим цикл на 8 итераций, и на нахождение остатка от деления тоже 8 итераций. При этом умножение, скорее, проиграет на 8-битном RNG и 8-битном процессоре по причине 16-битного суммирования в цикле против 8-битного вычитания, но выиграет - очень врядли (хотя, конечно, этот вопрос надо отдельно изучить, что оптимальнее на конкретной архитектуре). И при этом нахождение остатка не испортит распределение.

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

[Victor®]

Предположим, требуется моргать светодиодом с псевдослучайным интервалом от 2 до 30 секунд.
Посоветуйте простую функцию на Си для генерации целых чисел в таком малом диапазоне. Пусть от 1 до 255 максимум.
Функция rand() из библиотеки HT-PICC тяжеловата как-то...

Что-то все все усложняют...
На светодиод наверное человек смотреть будет?
Так вот через несколько он забудет, что видел при этом моргании.
Просто сделайте таблицу на пару минут и все...

Или игровой автомат? Тогда запись на видео лишит всего дохода...

П.С.
Или Вам надо офигенно нормальное распределение?
Увы - ничего не выйдет....

[Herz]

Не нужно вам никакого деления и остатка от деления.

Сразили! Я и слов-то таких не знаю: fractional. Буде Ваша воля пояснить свою идею попроще, постараюсь разобраться...

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

На данном этапе имеется PIC18, в котором аппаратное умножение присутствует. Но хочется перенести проект на PIC16, где на эту задачу желательно потратить минимум ресурсов.

Что-то все все усложняют...
На светодиод наверное человек смотреть будет?
Так вот через несколько он забудет, что видел при этом моргании.
Просто сделайте таблицу на пару минут и все...

Или игровой автомат? Тогда запись на видео лишит всего дохода...

П.С.
Или Вам надо офигенно нормальное распределение?
Увы - ничего не выйдет....

Да нет, светодиод - это пример, упрощение. На самом деле нужно посылать некий сигнал, но нерегулярно. Чтобы не привыкнуть. Таблица отпадает - памяти маловато. Хотя...
Но, от добра добра не ищут: пример от ув. SM очень даже ничего.

[Fat Robot]

Объясню предельно просто, не жонглируя терминологией

Пусть генератор после всех преобразований имеет выход X в диапазоне [0; (2^n)-1]
нужно и можно найти такие постоянные значения Y и Z, чтобы значение числа round((X*Y) >> Z) лежало в диапазоне [0; max-min]

Все, что для этого нужно это (((max-min) * X ) + (1 << (n - 1))) >> n.

Пусть n=7, a max-min=27
X в диапазоне [0; 127]
((X*27)+64)>>7 будет лежать в диапазоне [0;27]
выход операции умножения 16 бит

Если операции умножения у вас нет, но есть немного свободы в выборе значения для разности max-min, то нужно выбрать такое значение,
чтобы в нем было как можно меньше единиц в двоичных разрядах, тогда умножение из общего случая упростится в сумму сдвигов.

пусть вместо 27 из примера вас устроит 24 = (2^4) + (2^3), тогда преобразование будет таким:
(( X << 4 ) + ( X << 3 ) + 64 )>>7 и выход будет лежать в диапазоне [0;24]

Буде Ваша воля пояснить свою идею попроще, постараюсь разобраться...

[Herz]

Благодарю Вас, теперь ясно.

[Tanya]

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

? На, к примеру, 31 тоже ведь легко умножить.

[Guest_TSerg_*]

Использование чисто арифметических методов приведения диапазона ГСП к малым значениям может привести к весьма плачевным результатам в отношении качества их статистических свойств.

Лучше попробовать использование класса ГСП сразу генерирующего псевдо-случайные числа в желаемом амплитудном диапазоне (одновременно, к сожалению будет ограничен и период), к примеру класс LCG-генераторов:
Si = ( A*Sii + B ) mod M

К примеру, при M=83, A=2, B=0 получаем простой ГСП с периодом P=M-1=82, амплитудами 1..82:

Sii = 1..82; // Например Sii := 1;
---
Si := (Sii << 1) mod M;
Sii := Si;

Вместо mod M:
while Si<M do Si -= M;

или для данного A=2 даже:
if Si > M then Si -= M;

P.S.
Немногие удачные сочетания (A;M) c P=M-1:
(5;22), (2,29), ..., (6;41), .., (5;167), (2,5,6,8..;227)...

[SM]

Вместо mod M:
while Si<M do Si -= M;

мои глюки повторяем

while Si>=M do Si -= M;

[Guest_TSerg_*]

мои глюки повторяем

while Si>=M do Si -= M;

Ну да, очепытка, т.к. все не читал
Ниже-то - правильно.

[Herz]

Использование чисто арифметических методов приведения диапазона ГСП к малым значениям может привести к весьма плачевным результатам в отношении качества их статистических свойств.

Мне кажется, результат будет ровно тем же, если часть чисел из распределения с большей амплитудой возвращать в заданный диапазон вычитанием константы по условию. Распределение неминуемо исказится, ведь амплитуда связана с периодом жёстко. Или я не прав?

[SM]

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

Если считать саму вероятность появления того или иного кода, то тут будет выглядеть все очень просто:

допустим, изначально у любого кода вероятность 1/256.

Для случая, когда мы используем нахождение остатка, то для части кодов вероятность составит [256/N]/256, а для части ([256/N]-1)/256 ([] - ближайшее большее целое, тут я не знаю как спец-символ математический вставить), первому условию будут подчиняться первые 256 mod N кодов, 0....((256 mod N) - 1), второй части - оставшиеся N-(256 mod N) кодов (если я ничего не попутал). Причиной этому является то, что коды "заворачиваются в диапазон", перекрывая диапазон полностью, меньшее_целое(256/N) раз, а последний "заворот" происходит не полностью, перекрывая лишь первые (256 mod N) кодов, таким образом вероятность выпадения оставшихся кодов уменьшается. В худшем случае вдвое (1/256 и 2/256), ну и чем меньше N, тем меньше различие.

UPD:
Да, ну а если ГСЧ сразу строится по mod N, то у него изначально вероятности для всех генерируемых кодов равны, если, конечно, он построен правильно, то есть за весь его период каждый код выдается одинаковое кол-во раз.

[Herz]

Лучше попробовать использование класса ГСП сразу генерирующего псевдо-случайные числа в желаемом амплитудном диапазоне (одновременно, к сожалению будет ограничен и период), к примеру класс LCG-генераторов:
Si = ( A*Sii + B ) mod M

К примеру, при M=83, A=2, B=0 получаем простой ГСП с периодом P=M-1=82, амплитудами 1..82:

Sii = 1..82; // Например Sii := 1;
---

Надеюсь, для А=2, В=0 и М=83 правильно записал на Си:

Код

unsigned char rnd8(void)
{
static unsigned char seed = 1;
unsigned int temp;

temp = seed << 1;
if(temp > 83)    temp -= 83;
seed = (unsigned char)temp;

return temp;
}

P.S.
Немногие удачные сочетания (A;M) c P=M-1:
(5;22), (2,29), ..., (6;41), .., (5;167), (2,5,6,8..;227)...

Удачные в каком смысле? И как найдены?

[Guest_TSerg_*]

Сильно извиняюсь.. слегка забегаю в Инет, принимая поздравления за 58 лет

Удачные - в смысле тестирования ( автоматического подбора ) вариантов ГСП.
А вообще, есть ряд достаточно объективных тестов DIEHARD, NIST - для ГСП это норма, по другому трудно судить о качестве.

[Herz]

Сильно извиняюсь.. слегка забегаю в Инет, принимая поздравления за 58 лет

За 58 лет много уже приняли, небось? Присоединяюсь!

[Tanya]

А почему Вы не рассматриваете аппаратный генератор на реальном (природном) шуме?

[Herz]

А почему Вы не рассматриваете аппаратный генератор на реальном шуме?

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

[Tanya]

МК уже есть, а свободного места - почти нет.

И лишнего канала АЦП нет?

[Herz]

И лишнего канала АЦП нет?

Увы. Но, даже если бы был, без внешних компонентов ведь не обойтись. И диапазон менять сложнее.

[Guest_TSerg_*]

Спасибо за поздравление, Herz!

Однако нам не пристало даже после празднования - праздновать, работаем-да
Поэтому, "special for you" - небольшой эксклюзивный тест, предложенного мной LCG-генератора.

LCG/ГСП: M=83, A=2, B=0.

Генерация временной последовательности порождает вот такой график:
http://shot.qip.ru/00gZ9L-5OPovQkio/
(изменение числа инициализации 1..82 породит только фазовый сдвиг от числа 1)

Рассматривая график, понимаем, что явно присутствует зависимость ближайших от ближайших значений, а поэтому - устраиваем тест на автокорреляционную функцию:

Удостоверяемся, в том, что данный ГСП не может претендовать на крипто-статус.
http://shot.qip.ru/00gZ9L-6OPovQkip/

С другой стороны, разговор шел, как минимум, о желательности равновероятности чисел в избранном диапазоне.
Проверяем гистограмму.
http://shot.qip.ru/00gZ9L-5OPovQkir/

Что видим?
Идеально-равномерное распределение по числам в заданном диапазоне.
Т.е. каждое число из диапазона 1..81 встретилось только один раз за цикл.

Что и требовалось, как понимаю.

P.S.
Это всего лишь небольшой экскурс в практику RNG.
P.P.S.
Не сомневаюсь, что здесь присутствуют специалисты более высокого класса, чем я и мои объяснения они могут трактовать по своему.

[SM]

Кстати, на сколько я помню, для любого линейно-конгруэнтного генератора справедливо вот это:

для генератора X(n)=( a*X(n-1) + b ) mod m; длина последовательности P = m тогда (и только тогда), когда:

- b и m взаимно просты (их НОД = 1)
- a-1 кратно всем простым делителям m и a-1 кратно 4, если m кратно 4.

то есть, если эти условия выполнены, то a, b и m являются "удачными".

Ну и "претендовать на крипто-статус" линейно-конгруэнтный генератор не может в принципе... Это никак от АКФ не зависит, это зависит от того, на сколько сложно по какому-то количеству отсчетов генератора восстановить его функцию.

[Guest_TSerg_*]

Ну, так я и без претензий - просто сделал попытку объяснить/подобрать наипростейший RNG в диапазоне менее 8 бит с примитивным сохранением стат. значимости.

Ну и "претендовать на крипто-статус" линейно-конгруэнтный генератор не может в принципе... Это никак от АКФ не зависит,

Ишо как зависит

[SM]

Да и я без претензий... просто сделал попытку математически строго объяснить, что такое "удачные" a,b и m, и из какого множества их выбирать

Ишо как зависит

Вот у M-последовательностей, генерируемых LFSR, АКФ ну очень красивая... А на "крипто-статус" им тоже претендовать ну никак не удается

[Guest_TSerg_*]

Вот у M-последовательностей, генерируемых LFSR, АКФ ну очень красивая... А на "крипто-статус" им тоже претендовать ну никак не удается

Хм, кто бы сомневался

Лет так 40 назад, создал на основе всего лишь одного НИР большое подразделение - тестирование цифровых схем на основе сигнатурного анализа.
Как раз M-последовательности.
Как тест - работало

P.S.
Когда возникла необходимость передачи данных с 6 км глубины на поверхность океана, вот тут и встала крипто-проблема в полный рост. Есс-но - решили.
Это был 1988 г.

[Herz]

Поэтому, "special for you" - небольшой эксклюзивный тест, предложенного мной LCG-генератора.

Очень Вам признателен. На редкость познавательной получилась тема для начинающих, ИМХО.

Когда возникла необходимость передачи данных с 6 км глубины на поверхность океана, вот тут и встала крипто-проблема в полный рост. Есс-но - решили.
Это был 1988 г.

Это чтобы дельфины сообщение не перехватывали?

[Guest_TSerg_*]

Очень Вам признателен.

Да не за, что - гости разошлись, а "моcк" всегда требует корма, увы

Это чтобы дельфины сообщение не перехватывали?

Ага, такие небольшие дельфинчики от Office of Naval Intelligence ( ONI ).

[ARV]

для задач типа случайного мигания светодиодом, т.е. для тех, где критерий случайности определяет человек-наблюдатель, в сильно стесненных условиях я использую тупо содержимое flash-памяти МК в качестве псевдослучайных чисел.

Код

char rand(void){
      static __flash char *ptr;
      return ++ptr*;
}

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