Анализ ГПСЧ на Haskell

а по заданию это запрещено

job

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

Нет, в данном случае как раз-таки особенности генератора использовать можно. Один из вариантов решения - «вхолостую» сгенерить 2*m^2 чисел, зафиксировать пару и генерить дальше до совпадения. Но на хаскеле нет циклов, а как более-менее красиво это сделать рекурсивно, я не знаю. Об этом и тред.

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

Обоснуй.

Не могу обосновать, я это предполагал, полагая, что это требуется сделать для абстрактного генератора, а не конкретного.

Какой-то он у тебя неканоничный. Это не LCG, а MRG.

Попробуй почитать это, что ли http://www.sciencedirect.com/science/article/pii/S0377042705003596

Стойкость таких гпсч основана как раз на том, что мы все умрем раньше.

Я знаю, что это не LCG, просто не знал, как это называется. Спасибо за линк. Как искать я в принципе знаю, написал выше, мне интересно, как это реализовать без тонны кода и извращений?

В job я это не написал, потому что «курсовая» часть работы уже сделана. У меня есть вариант, занимающийся анализом генератора на C++ и Qt. Сейчас я хочу понять, как это сделать с упором на функциональную парадигму.

Хоть на m есть ограничения? Типа степень двойки или что.

Хотя, кажется, у тебя фибоначчиев генератор (вроде непростые допускаются). Глянь еще вот: https://books.google.com.ua/books?id=wj5OyydZ5bkC&pg=PA62&lpg=PA62&am...

можешь попробовать найти хар. полином и посмотреть, что выходит.

Почему бы просто не взять любой алкоритм поиска циклов? https://en.m.wikipedia.org/wiki/Cycle_detection

Вот на русском языке более полная статья: http://mech.math.msu.su/~fpm/ps/k10/k106/k10609.pdf

Твой генератор это функция f от одного аргумента: пары чисел (x_n,x_{n-1}), действующая так:

(x_{n-1},x_{n-2}) преобразуется в (x_n, x_{n-1})

При этом применим любой алгоритм для нахождения циклов в функции f(x), описанный в статьях выше, кроме тех, которые основываются на отношении порядка/частичного порядка на множестве X={x}.

В функциональщину здесь очень простая дверь: вместо того, чтоб вычислять f(x_n) ты обращаешься к f_stream[x_n] - элементу потока псевдослучайных чисел, с индексом x_n.

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

Лень сейчас писать более подробно, но идея думаю понятна.

Наговнякал, не приходя в сознание:

type GenState = (Int, Int)

genNext :: Int -> Int -> Int -> GenState -> GenState
genNext a b m (xn, xn') = let xn'' = (a * xn + b * xn') `mod` m
                          in xn'' `seq` (xn', xn'')

floyd :: Eq a => (a -> a) -> a -> (Int, a)
floyd f x0 = z `seq` (go2 (f z) 1, z)
  where
    z = go1 (f x0) (f (f x0))
    go1 x y | x == y = x
            | otherwise = go1 (f x) (f (f y))
    go2 x n | x == z = n
            | otherwise = go2 (f x) (n + 1)

main :: IO ()
main = do
  let a = 2
      b = 3
      m = 23
      x0 = (1, 1)
  print . fst $ floyd (genNext a b m) x0 -- Output: 176

Думаю, на списках было бы идиоматичнее сделать. Но, как сделал, так сделал.

Обзор прекрасен, спасибо.

Можно через монаду Reader делать, с другой стороны, можно сделать отдельный тип данных — вычетов по модулю m (на один параметр меньше).

Э-э-э, cryptol? Топик не читал.

1. Не понял какое-такое с - в формулах его не вижу 2. Лениво разбираться как писать генератор и период каррированными - но Haskell ленивый язык, лениться можно :) Написал кортежными - на суть не влияет 3. Непонятно в чем сложности - какие-то потоки/списки/шмиски... Функциональная парадигма какая-то... Пляшем всегда от алгоритма - тупо итеративный перебор в лоб. В языке где нет циклов - через - правильно, хвостовую рекурсию, вестимо:

a = 1
b = 1
m = 10^6

gen :: (Int, Int) -> (Int, Int)
gen (x, y) = (y, (a*y + b*x) `mod` m)

period :: (Int, Int) -> Int
period p = go 1 $ gen p where
    go i n | n == p = i
           | otherwise = go (i+1) $ gen n

main = print $ period (0, 1)

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

Не понял какое-такое с - в формулах его не вижу

Есть такая контора Galois Inc. Хорошо известна хаскелевским олдфагам, так как была в своё время «одной из немногих».

Так вот, сия контора выпускает по заказу американской военщины продукт Cryptol — написанный на хаскеле DSL для исследования криптографических алгоритмов. Недавно они его заопенсорсили.

Это все очень интересно, но не проливает свет на параметр «с» в стартовом посте и его роль в формуле и мировой истории.

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

Прошу прощения, там должно быть (a*x{n-1}+b*x{n-2}+c)%m Спасибо, я как раз не мог понять, как это на хвостовой рекурсии написать

Внезапно, даже при c=0 этот код будет работать не всегда. В общем случае, начиная с пары x1,x2 мы не всегда к ней же и придем. В общем случае траектория генератора визуально будет напоминать греческую букву Ро, а не латинскую O.

Это связано с понятием «первообразного элемента по модулю». Лень гуглить, при каком соотношении между a,b,c,m любая траектория генератора будет строгим циклом.

Это решается началом не с первой пары, а с пары номером m^2

Разумеется, если взять к примеру такой волшебный генератор:

gen :: (Int, Int) -> (Int, Int)
gen (x, y) = (y, y+1)

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

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

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

Знаю неплохой математический форум для подобных вопросов - dxdy.ru Порядки там строгие, халявщиков без собственных мыслей сразу в карантин до бана, но немалая часть участников - квалифицированные специалисты (например, преподы МГУ и всяких там зарубежных университетов). Если действительно интересно и готов работать (а не как в этой теме - памагите в тривиальной задаче без собственных попыток решения) - то самое то, лучше в рунете нет (заявляю ответственно).

Спасибо за помощь :) Если всё-таки возьму на себя проект, то пойду туда, если понадобится.

памагите в тривиальной задаче без собственных попыток решения

Задача действительно тривиальная, но я пока ещё не научился мыслить в терминах ФП и хаскеля

Прокрутить генератор m^2 раз это не тоже самое, что прокрутить его m или 2m раз. Следует просто применить любой алгоритм поиска циклов типа Флойда. Не забуду если, к вечеру закину и обычную и функциональную реализации - я такое уже делал когда-то. Только не обессудь - на Scheme.

Мне, в общем-то, главное понять суть реализации, остальное я должен осилить

А почему в данном случае не подходит? Всё равно величина периода в худшем случае порядка m^2

Я её всё-таки написал, получилось не очень красиво, но всё же:

module Random where

type Generator = (Int, Int, Int, Int, Int, Int)

gen :: Generator -> (Int, Int) -> (Int, Int)
gen (a, b, c, m, x0, x1) (xp, x) = (x, (a*xp+b*x) `mod` m)

skip :: Generator -> Int -> (Int, Int)
skip (a, b, c, m, x0, x1) length  | length==1 = (x0, x1)
                                  | otherwise = skip (a, b, c, m, fst $ gen (a, b, c, m, x0, x1) (x0, x1), snd $ gen (a, b, c, m, x0, x1) (x0, x1)) (length-1)

count :: Generator -> (Int, Int) -> (Int, Int) -> Int -> Int
count gener (x0, x1) (x2, x3) cnt | (x0==x2)&&(x1==x3) = (cnt+1)
                                  | otherwise = count gener (x0, x1) (gen gener (x2, x3)) (cnt+1)

period :: Generator -> Int
period (a, b, c, m, x0, x1) = count (a, b, c, m, x0, x1) (skip (a, b, c, m, x0, x1) (m*m)) (gen (a, b, c, m, x0, x1) (skip (a, b, c, m, x0, x1) (m*m))) 0

type Val = (Int, Int)
type Generator = Val -> Val
type GenVal = (Generator, Val)

gen :: Int -> Int -> Int -> Int -> Generator
gen a b c m = go where
    go (x, y) = (y, (a*x+b*y+c) `mod` m)

skip :: GenVal -> Int -> Val
skip (g, v) = go v where
    go v i | i <= 0    = v
           | otherwise = go (g v) (i-1)

period :: GenVal -> Int
period gv@(g, _) = go 1 $ g s where
    s = skip gv 1000 -- все равно не знаем m )
    go i n | n == s    = i
           | otherwise = go (i+1) $ g n

main = print $ period (gen 2 3 0 23, (1, 1))

ЗЫ генератор может содержать в себе стартовое значение и возвращать в качестве результата пару - новое значение и новый генератор, но мне мой вариант показался проще и универсальнее - генератор определяется параметрами, а стартовое значение задается отдельно. При больших m надо компилировать с -О2 или самому принудительно ставить строгие вычисления параметров хвостовых рекурсий. Но все равно это лажа - мы не знаем m (по условию) и не знаем сколько скипить. Короче, ТС, раздели задачи - составление алгоритма и его реализацию.

Во, золотые слова. Слушайте этого анона, и не будет skip 1000 у вас в коде.

Если кто научит, как избавиться от тучи лишних параметров, поменяю

Правильный ответ — декомпозиция. Более конкретно, в данном случае надо разделить код, который ищет цикл, от кода, который описывает генератор. И побольше пользоваться хорошей штукой под названием «полиморфизм», а так же не стесняться замыканий.

Посмотри мой код выше, может станет понятнее.

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

Нахождение периода взял из Вашего кода нахождение периода и дописал нахождение предпериода. Осталось добавить расчёт отклонения распределения от нормального. В связи с этим возник вопрос: как это написать на иммутабельных массивах? Суть алгоритма - сгенерить 400 чисел, разбить диапазон генерации на 20 частей и посчитать, сколько чисел из каждого интервала. Не могу понять, как это реализовать?

То есть задача - на основе массива создать новый, в котором будет изменён(увеличен) ровно один элемент

на основе массива создать новый, в котором будет изменён(увеличен) ровно один элемент

Соответствующие функции есть в кажом API массивов. Для примера: раз и два. Рекомендую пользоваться вторым пакетом, там много плюшек всяких хороших. Кстати, можно и Data.Map поэксплуатировать.

Благодарю

Наговнякал пример:

import System.Random
import qualified Data.Vector as V

getStats :: Int -> [Int] -> V.Vector Int
getStats len indices = V.accum (+) (V.replicate len 0) (zip indices (repeat 1))

genSmallInts :: StdGen -> Int -> [Int]
genSmallInts seed len = randomRs (0, len - 1) seed

main :: IO ()
main = do
  let len = 20
      tosses = 400
  seed <- newStdGen
  print . getStats len . take tosses $ genSmallInts seed len

Не подскажете, как использовать замыкания в монадическом коде? Весь функциональный код стал намного чище, а вот с монадической частью проблема. Как сократить предпоследнюю строку тут? http://pastebin.com/c4Uqy9L6

Прям точно так же и делать: разделять код. Что-то в стиле:

recalc :: Builder -> IO ()
recalc builder = do
  let doMyJob label = do
        entry <- builderGetObject builder castToEntry
        text <- entryGetText entry
        return $ read text
      labels = ["a", "b", "c", "m", "x1", "x2"]
  [a, b, c, m, x1, x2] <- mapM doMyJob ["entry" ++ l | l <- labels]
  area <- builderGetObject builder castToDrawingArea "drawingarea1"
  w <- widgetGetAllocatedWidth area
  h <- widgetGetAllocatedHeight area
  putStrLn $ research gener (x1, x2)
  area `on` draw $ (drawHistogram w h (gen a b c m) (x1, x2) 100 m)
  widgetQueueDraw area

Нашёл опечатку:

entry <- builderGetObject builder castToEntry label

Январь. Делать курсач. Это пять. Хоть сдал?

Вот, наговнякал, если Вам интересно: https://github.com/SeTSeR/HaskellRandom Последнее, о чём Вы сказали, я ещё доделаю

Это не курсач, я ещё школьник. Да, вот такие у нас домашние задания по инфе.

Сдавать буду во вторник