Programowanie funkcyjne/Strumienie

From Studia Informatyczne

Spis treści

Wstęp

Przedstawiając konstrukcje imperatywne, powiedzieliśmy, że w programowaniu imperatywnym upływ czasu w świecie rzeczywistym jest modelowany przez upływ czasu w programie komputerowym. W "czystym" programowaniu funkcyjnym nie ma zmiennych, ani zmian stanów. Trudno też mówić o upływie czasu w programie. Wszak program składa się z implementacji pojęć matematycznych. Jedyna zależność związana z czasem polega na tym, że pojęcia wykorzystywane muszą być zdefiniowane przed pojęciami korzystającymi z nich.

Strumienie to po prostu ciągi wartości, być może nieskończone. Programowanie strumieniowe pozwala jednak na utworzenie związku między czasem i obliczeniami.

W fizyce mamy pojęcie czasoprzestrzeni. Sama czasoprzestrzeń nie zmienia się w czasie -- jest to twór matematyczny ujmujący upływ czasu i historię całej przestrzeni. Historia pojedynczej cząstki jest reprezentowana przez jej trajektorię w czasoprzestrzeni, nazywaną linią świata cząstki.

Jak zobaczymy, strumienie mogą być użyte do przedstawiania historii obliczeń iteracyjnych. Strumień może reprezentować całą, być może nieskończoną historię pewnych obliczeń. Pozwala to "jednym rzutem oka" spojrzeć na całą historię obliczenia. Jak zobaczymy, pozwoli to również na implementowanie pewnych ciekawych operacji na "liniach świata obliczeń".

Wiele z podanych w tym wykładzie przykładów pochodzi z [AS], przy czym zostały przeniesione z języka Scheme do Ocamla.

Definiowanie form specjalnych

Jak już pisaliśmy wcześniej, w Ocamlu argumenty procedur są obliczane gorliwie, tzn. najpierw są obliczane ich wartości, a dopiero potem przystępujemy do obliczania wartości procedury. Wyjątkiem są formy specjalne, np. instrukcja warunkowa if-then-else --- w zależności od wartości warunku, tylko jeden z pozostałych członów jest obliczany.

Czy można takie formy specjalne zaimplementować samemu? Tak. Potrzeba do tego dwóch rzeczy. Musimy umieć odroczyć moment obliczania wartości i wprowadzać własne makrodefinicje.

O uleniwianiu i odraczaniu obliczeń mówiliśmy w poprzednim wykładzie. W Ocamlu możemy nie tylko definiować własne makrodefinicje, ale wręcz zmieniać składnię języka (sic!). Za analizę składniową odpowiedzialny jest preprocesor P4. Nie będziemy tu mówić o nim, gdyż wykracza to poza zakres tego kursu. Jednak raz użyjemy definicji wprowadzającej prostą makrodefinicję, która realizuje formę specjalną.

Implementacja strumieni

Formalnie strumień to ciąg elementów --- być może nieskończony. Pomysł polega jednak na tym, żeby był on "leniwy", tzn. obliczane były tylko potrzebne wartości. Przyjmujemy przy tym, że wartość stają się potrzebne zgodnie z ich kolejnością występowania w strumieniu. Tak więc w każdej chwili strumień składa się ze skończonego ciągu obliczonych (potrzebnych) wartości oraz odroczonego (być może nieskończonego) ciągu pozostałych elementów. Implementując strumienie korzystamy z uleniwiania (zaimplementowanego przez moduł Lazy). Strumień, jeśli nie jest pusty, to para: wartość pierwszego elementu i odroczony strumień pozostałych wartości.

type 'a stream = Nil | Cons of 'a * 'a stream Lazy.t;;
type 'a stream = Nil | Cons of 'a * 'a stream Lazy.t

let empty s = s = Nil;;
val empty : 'a stream -> bool = <fun>

Podstawowe "cegiełki" do budowy strumieni są dokładnie takie same jak do budowy list: pusty strumień oraz dołączenie elementu na początku strumienia. Dwie podstawowe operacje służące do dekonstrukcji strumieni, to procedury zwracające głowę i ogon danego strumienia.

let head s = 
  match s with
    Nil -> failwith "Empty" |
    Cons (h, _) -> h;;
val head : 'a stream -> 'a = <fun>

let tail s = 
  match s with
    Nil -> failwith "Empty" |
    Cons (_, t) -> Lazy.force t;;
val tail : 'a stream -> 'a stream = <fun>

Tak zdefiniowany typ danych jest trochę niewygodny w użyciu, gdyż zawsze gdy używamy Cons, powinniśmy go używać razem z lazy. Aby ułatwić sobie życie, zdefiniujemy Cons jako makro, które przed drugim swoim parametrem dodaje słowo kluczowe lazy.

\texttt{Cons}(x, y) \equiv  \texttt{Cons}(x, \texttt{lazy}\ y)

Oto niezbędne "zaklęcie", które podajemy bez wyjaśnień. Zainteresowanych odsyłamy do dokumentacji preprocesora P4.

#load "camlp4o.cma";;
#load "pa_extend.cmo";;
#load "q_MLast.cmo";;

EXTEND
  Pcaml.expr: LEVEL "simple"
  [ [ UIDENT "Cons"; param=SELF ->
        let h,l =
          match param with <:expr< ($list:[h;l]$) >> -> (h,l)
                           | _ -> raise Not_found
        in
          <:expr< Cons ($h$) (lazy $l$)>> ] ]
  ;
END;;

Operacje na strumieniach

Mając do dyspozycji takie konstruktory i selektory możemy zdefiniować kilka pomocniczych operacji na strumieniach:

let rec const_stream c = 
  Cons (c, const_stream c);; 
val const_stream : 'a -> 'a stream = <fun>

Definiuje nieskończony strumień, którego wszystkie elementy są równe c.


let ones = const_stream 1;;
val ones : int stream = Cons (1, <lazy>)

To jest nieskończony strumień jedynek.


let rec filter_stream p s = 
  if empty s then 
    Nil
  else 
    let h = head s 
    and t = tail s
    in 
      if p h then 
        Cons (h, filter_stream p t)
      else
        filter_stream p t;;
val filter_stream : ('a -> bool) -> 'a stream -> 'a stream = <fun>

Działa podobnie jak procedura filter dla list, tworzy strumień złożony wyłącznie z tych elementów danego strumienia, które spełniają podany predykat. Jeśli strumień będący parametrem jest nieskończony, to nieskończenie wiele jego elementów powinno spełniać predykat, tak żeby strumień wynikowy był również nieskończony.


let rec stream_ref s n = 
  if n = 0 then 
    head s
  else
    stream_ref (tail s) (n - 1);;
val stream_ref : 'a stream -> int -> 'a = <fun>

Zwraca element strumienia z podanej pozycji. Pozycje w strumieniu są indeksowane od 0.


let rec stream_map f s = 
  if empty s then 
    Nil
  else
    Cons (f (head s), stream_map f (tail s));;
val stream_map : ('a -> 'b) -> 'a stream -> 'b stream = <fun>

Działa podobnie jak procedura map dla list, stosuje dane przekształcenie f do wszystkich elementów danego strumienia i tworzy strumień wyników.

let scale_stream c s = 
  stream_map (function x -> x *. c) s;;
val scale_stream : float -> float stream -> float stream = <fun>

Procedura ta mnoży elementy danego strumienia przez zadaną stałą.

let rec stream_map2 f s1 s2 = 
  if empty s1 or empty s2 then Nil 
  else Cons (f (head s1) (head s2), stream_map2 f (tail s1) (tail s2));;
val stream_map2 : ('a -> 'b -> 'c) -> 'a stream -> 'b stream -> 'c stream = <fun>

Jest to odpowiednik procedury map2 z modułu List. Pozwala ona na proste zdefiniowanie wielu operacji łączących strumienie. Na przykład, dodawanie strumieni:

let rec add_int_streams s1 s2 = stream_map2 (+) s1 s2;;
val add_int_streams : int stream -> int stream -> int stream = <fun>
 
let sample s = 
 let rec pom i s = 
   if i = 0 || empty s then []
   else head s :: pom (i-1) (tail s)
 in 
   pom 12 s;;
val sample : 'a stream -> 'a list = <fun>

Procedury tej będziemy używać do testowania konstruowanych przez nas strumieni.

sample ones;;
- : int list = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1]
Uwaga
W tym wykladzie pokazujemy, jak zaimplementować strumienie, choć Ocaml zawiera moduł (Streams) implementujący strumienie. Jednak są to dwie różne struktury danych. Strumienie zaimplementowane w module Streams odpowiadają bardziej pojęciu strumienia w C++.

Przykłady strumieni nieskończonych

Zdefiniujmy strumień liczb naturalnych:

let rec integers_from x = 
  Cons(x, integers_from (x+1));;
val integers_from : int -> int stream = <fun>

let nats = integers_from 0;;
val nats : int stream = Cons (0, <lazy>)

sample nats;;
- : int list = [0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11]

Podobnie możemy zdefiniować strumień liczb Fibonacciego:

let fibs = 
  let rec fibo a b = 
    Cons (a, fibo b (a+b))
  in fibo 0 1;;
val fibs : int stream = Cons (0, <lazy>)

sample fibs;;
- : int list = [0; 1; 1; 2; 3; 5; 8; 13; 21; 34; 55; 89]

Jak widać, jeżeli jesteśmy w stanie skonstruować iteracyjną procedurę wyliczającą kolejne elementy strumienia, to tym samym jesteśmy w stanie zdefiniować strumień tych wartości. Procedura taka stanowi ukrytą w strumieniu "maszynerię", która na żądanie podaje kolejne elementy strumienia. Takie definicje strumieni nazywamy nie uwikłanymi (w odróżnieniu od uwikłanych, które poznamy dalej). Charakterystyczne dla definicji nie uwikłanych jest to, że bytem definiowanym rekurencyjnie jest procedura, a nie strumień. Strumień jest zdefiniowany dopiero w oparciu o rekurencyjnie zdefiniowaną procedurę.

Strumienie to tylko jeden ze sposobów konstruowania złożonych struktur danych. Można więc go łączyć z dowolnymi innymi sposobami konstruowania danych i tworzyć strumienie list, listy strumieni, czy nawet strumienie strumieni. Oto przykład, strumień strumieni powstałych z danego (nieskończonego) strumienia przez usunięcie coraz większej liczby początkowych elementów (poczynając od zera elementów):

let rec tails s = 
  if empty s then Cons (s, Nil) else Cons (s, tails (tail s));;
val tails : 'a stream -> 'a stream stream = <fun>

Strumień liczb pierwszych -- sito Eratostenesa

Spróbujmy skonstruować nieskończony strumień liczb pierwszych, generowanych metodą sita Eratostenesa. Sposób konstrukcji takiego strumienia możemy przedstawić w postaci schematu przypominającego schematy blokowe układów przetwarzających sygnały. Strumienie są na nim zaznaczone liniami ciągłymi, a pojedyncze wartości przerywanymi.

Poniższy diagram przedstawia budowę procedury sito implementującej pojedyncze sito. Zauważmy, że jest on rekurencyjny, tzn. procedura ta pojawia się również na tym schemacie. Odpowiada to temu, że liczby przesiane przez jedno sito trafiają na kolejene sito, itd.

Schemat ten przekłada się na następujący program:

let sitko p s = 
  filter_stream (function x -> (x mod p) <> 0) s;;
val sitko : int -> int stream -> int stream = <fun>

let rec sito s  = 
  Cons (head s, sito (sitko (head s) (tail s)));;
val sito : int stream -> int stream = <fun>
 
let primes = 
  sito (integers_from 2);;
val primes : int stream = Cons (2, <lazy>)

sample primes;;
- : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37]

Dzięki leniwości strumieni, w zależności od tego ile liczb odczytamy ze strumienia, powstanie odpowiednia liczba sit i odpowiednia ilość liczb zostanie przez nie przesiana.

Definicje uwikłane

Definiując nieskończone strumienie nie musimy tego zawsze robić poprzez podanie odpowiedniej procedury rekurencyjnej (por. integers-from, fibo, sieve). Możemy rekurencyjnie użyć strumienia do zdefiniowania jego samego.

let rec ones =
  Cons (1, ones);;

sample ones;;
- : int list = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1]

Do budowy bardziej skomplikowanych strumieni potrzebujemy dodatkowych operacji tworzących strumienie. Na przykład, do utworzenia strumienia liczb naturalnych czy liczb Fibonacciego potrzebujemy dodawania strumieni.

let rec nats =
  Cons (0, add_int_streams ones nats);;
val nats : int stream = Cons (0, <lazy>)

sample nats;;
- : int list = [0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11]

Sposób obliczania kolejnych elementów strumieni zdefiniowanych w sposób uwikłany najlepiej przedstawić sobie w postaci tabelki:


Zamiast dodawania strumienia jedynek, możemy użyć operacji zwiększania o jeden:

let succ x = x + 1;;
val succ : int -> int = <fun>

let rec nats = 
  Cons(0, stream_map succ nats);;
val nats : int stream = Cons (0, <lazy>)

sample nats;;
- : int list = [0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11]


A oto uwikłana definicja liczb Fibonacciego:

let rec fibs = 
  Cons (0, Cons (1, add_int_streams fibs (tail fibs)));;
val fibs : int stream = Cons (0, <lazy>)

sample fibs;;
- : int list = [0; 1; 1; 2; 3; 5; 8; 13; 21; 34; 55; 89]


Możemy też w uwikłany sposób zdefiniować strumień liczb pierwszych. Użyjemy do tego predykatu prime sprawdzającego, czy liczba jest pierwsza ... a predykat prime zdefiniujemy używając strumienia liczb pierwszych.

let divisible x y = x mod y = 0;;
val divisible : int -> int -> bool = <fun>

let square x = x * x;;
val square : int -> int = <fun>

let rec primes = 
  Cons (2, filter_stream prime (integers_from 3))
and prime n = 
  let rec iter ps = 
    if head ps * head ps > n then true
    else if divisible n (head ps) then false
    else iter (tail ps)
  in 
    iter primes;;
val primes : int stream = Cons (2, <lazy>)
val prime : int -> bool = <fun>

sample primes;;
- : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37]

Całość działa poprawnie, ponieważ strumień jest leniwą strukturą danych. Pierwszy element strumienia liczb pierwszych, 2, jest dany explicite. Natomiast sprawdzając czy kolejne liczby są pierwsze, zawsze mamy już obliczone potrzebne do tego liczby pierwsze.

Konkretnie, podana explicite pierwsza liczba pierwsza, 2, pozwala sprawdzić pierwszość liczb do 4. To wystarczy, żeby stwierdzić, że 3 jest liczbą pierwszą. W każdej chwili, największa obliczona dotąd liczba pierwsza p wystarcza do przetestowania pierwszości liczb do p^2 włącznie. Natomiast kolejna liczba pierwsza p' spełnia p' < 2p \le p^2. Tak więc, każda kolejna wyznaczona liczba pierwsza zwiększa możliwy zakres testowania pierwszości o więcej niż odległość do kolejnej liczby pierwszej.


Iteracje jako strumienie

W przypadku definicji nie uwikłanych używaliśmy procedur iteracyjnych do zdefiniowania strumieni. Możemy ten mechanizm odwrócić i użyć strumieni do opisania procesów iteracyjnych -- kolejne elementy strumienia mogą reprezentować kolejne kroki iteracji.

Metoda Newtona przybliżania pierwiastków

Oto przykład, strumień kolejnych przybliżeń pierwiastka kwadratowego, metodą Newtona:

let sqrt_improve x g = (g +. x /. g) /. 2.0;;
val sqrt_improve : float -> float -> float = <fun>

let sqrt_stream x = 
  let rec guesses = 
    Cons (1.0, 
      stream_map (sqrt_improve x) guesses)
  in guesses;;
val sqrt_stream : float -> float stream = <fun>

sample (sqrt_stream 42.0);;
- : float list = 
[1.; 21.5; 11.7267441860465116; 7.65415047676148141; 
 6.57068474328365859; 6.48135630633341897; 6.48074072764349385; 6.48074069840786; 
 6.48074069840786; 6.48074069840786; 6.48074069840786; 6.48074069840786]

Przybliżenia \pi

Skonstruujemy strumień przybliżeń liczby \pi. Użyjemy do tego celu szeregu zbieżnego do \frac{\pi}{4}:

\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots

Najpierw definiujemy strumień mianowników (1, -3, 5, -7, \dots), następnie przekształcamy go w strumień sumowanych ułamków, po czym wyznaczamy strumień sum częściowych naszego szeregu.

let pi_summands = 
  let succ x = 
    if x > 0.0 then 
      -.x -. 2.0
    else
      -.x +. 2.0
  in
    let rec s = 
      Cons (1.0, stream_map succ s)
    in
      stream_map (fun x -> 1.0 /. x) s;;
val pi_summands : float stream = Cons (1., <lazy>)

sample pi_summands;;
- : float list =
[1.; -0.333333333333333315; 0.2; -0.142857142857142849; 0.111111111111111105;
 -0.0909090909090909116; 0.0769230769230769273; -0.0666666666666666657; 
 0.0588235294117647051; -0.0526315789473684181; 0.0476190476190476164;
  -0.0434782608695652162]

let rec add_float_streams s1 s2 = stream_map2 (+.) s1 s2;;
val add_float_streams : float stream -> float stream -> float stream = <fun>

let partial_sums s = 
  let rec ps = 
    Cons (head s, add_float_streams (tail s) ps)
  in ps;;
val partial_sums : float stream -> float stream = <fun>

Strumień przybliżający \pi uzyskujemy mnożąc elementy strumienia sum częściowych przez 4.

let pi_stream = scale_stream 4.0 (partial_sums (pi_summands));;
val pi_stream : float stream = Cons (4., <lazy>)

sample pi_stream;;
- : float list =
[4.; 2.66666666666666696; 3.46666666666666679; 2.89523809523809561;
 3.33968253968254025; 2.97604617604617649; 3.28373848373848443;
 3.01707181707181782; 3.25236593471887669; 3.04183961892940324;
 3.23231580940559393; 3.05840276592733318]

Taki strumień jest co prawda zbieżny, ale bardzo wolno. Po zsumowaniu pierwszych 1000 elementów ustalone są dopiero trzy pierwsze cyfry: 3.14. Jest to dosyć oczywiste, jeżeli zauważymy, że suma pierwszych n elementów jest obarczona błędem rzędu \frac{1}{n}.

Można "przyspieszyć" zbieżność tego szeregu stosując tzw. "akcelerator Eulera". Jeśli S_n jest sumą pierwszych n wyrazów szeregu, to szereg zakcelerowany ma postać:

S_{n+1} - \frac{(S_{n+1} - S_n)^2}{S_{n-1} - 2 S_n + S_{n+1}}

Działa on dobrze np. dla ciągów naprzemiennie zbieżnych, ale nie tylko. Oto jego implementacja:

let euler_transform st =
  let transform s = 
    let s0 = stream_ref s 0
    and s1 = stream_ref s 1
    and s2 = stream_ref s 2
    in (s2 -. (s2 -. s1) *. (s2 -. s1) /. (s0 -. 2.0 *. s1 +. s2))
  in
    stream_map transform (tails st);;
val euler_transform : float stream -> float stream = <fun>

sample (euler_transform pi_stream);;
- : float list =
[3.16666666666666696; 3.13333333333333375; 3.14523809523809561;
 3.13968253968254; 3.14271284271284346; 3.14088134088134163;
3.14207181707181782; 3.14125482360776553; 3.14183961892940333;
3.1414067184965031; 3.14173609926066666; 3.14147968900425623]

Taki strumień jest już zbieżny w rozsądnym czasie. Przeliczenie 1000 elementów przyspieszonego szeregu daje przybliżenie 3.141592653. Jeszcze lepsze wyniki daje wielokrotne przyspieszanie. Skonstruujmy strumień kolejnych przyspieszeń strumienia sum częściowych i wybierzemy z niego strumień pierwszych elementów:

let rec pi_table = 
  Cons (pi_stream, stream_map euler_transform pi_table);;
val pi_table : float stream stream = Cons (Cons (4., <lazy>), <lazy>)

let pi_stream_acc = stream_map head pi_table;;
val pi_stream_acc : float stream = Cons (4., <lazy>)

sample pi_stream_acc;;
- : float list =
[4.; 3.16666666666666696; 3.14210526315789496; 3.14159935731900486;
 3.14159271403377849; 3.14159265397529275; 3.14159265359117645;
 3.14159265358977802; 3.14159265358979534; 3.14159265358979489; nan; nan]

let pi = stream_ref pi_stream_acc 9;;
- : float = 3.14159265358979489
      

Jak widać, już dziewiąty element takiego strumienia stanowi dobre przybliżenie \pi, gdyż błąd jest dopiero na 15 miejscu po przecinku. Dalsze przybliżanie nie daje już rezultatów, gdyż ze względu na ograniczoną precyzję obliczeń, w mianownikach liczonych ułamków pojawiają sie zera.

Trójkąt Pascala

Strumień strumieni możemy sobie przedstawić jako nieskończoną dwuwymiarową tablicę elementów. Można więc zdefiniować trójkąt Pascala, jako strumień strumieni liczb całkowitych.

Zdefiniujemy go w sposób uwikłany. Pierwsza kolumna zawiera same jedynki. Każda kolejna kolumna zaczyna się od jedynki, a pozostałe elementy uzyskujemy sumując poprzedzający element w tej samej kolumnie oraz element na tej samej pozycji w poprzedzającej kolumnie.

W pewnym sensie, wyraziliśmy tutaj dwie zagnieżdżone nieskończone pętle: jedna obliczająca kolejne elementy strumieni, a druga obliczająca kolejne strumienie. Dzięki odroczeniu obliczania elementów strumieni, obie pętle wykonują faktycznie tylko tyle kroków, ile jest konieczne do obliczenia badanych elementów.

let rec pascal = 
  let next s = 
    let rec wyn = 
      Cons (1, add_int_streams (tail s) wyn)
    in 
      wyn
  in 
    Cons (ones, stream_map next pascal);;
val pascal : int stream stream = Cons (Cons (1, <lazy>), <lazy>)

List.map sample (sample pascal);;
- : int list list =
[[1;  1;  1;   1;    1;    1;     1;     1;     1;      1;      1;      1];
 [1;  2;  3;   4;    5;    6;     7;     8;     9;     10;     11;     12];
 [1;  3;  6;  10;   15;   21;    28;    36;    45;     55;     66;     78];
 [1;  4; 10;  20;   35;   56;    84;   120;   165;    220;    286;    364];
 [1;  5; 15;  35;   70;  126;   210;   330;   495;    715;   1001;   1365];
 [1;  6; 21;  56;  126;  252;   462;   792;  1287;   2002;   3003;   4368];
 [1;  7; 28;  84;  210;  462;   924;  1716;  3003;   5005;   8008;  12376];
 [1;  8; 36; 120;  330;  792;  1716;  3432;  6435;  11440;  19448;  31824];
 [1;  9; 45; 165;  495; 1287;  3003;  6435; 12870;  24310;  43758;  75582];
 [1; 10; 55; 220;  715; 2002;  5005; 11440; 24310;  48620;  92378; 167960];
 [1; 11; 66; 286; 1001; 3003;  8008; 19448; 43758;  92378; 184756; 352716];
 [1; 12; 78; 364; 1365; 4368; 12376; 31824; 75582; 167960; 352716; 705432]]



Scalanie strumieni

Zdefiniujmy pożyteczną procedurę, która scala dwa strumienie posortowane niemalejąco w jeden:

let rec merge s1 s2 = 
  if empty s1 then s2 else
  if empty s2 then s1 else 
    let h1 = head s1 
    and h2 = head s2
    in 
      if h1 < h2 then 
        Cons (h1, merge (tail s1) s2)
      else  
        Cons (h2, merge s1 (tail s2));;
val merge : 'a stream -> 'a stream -> 'a stream = <fun>

Procedurę tę możemy wykorzystać do zdefiniowania w sposób uwikłany strumienia postaci: jedna jedynka, dwie dwójki, trzy trójki itd.

let rec s = Cons (1, merge (stream_map (fun x -> x+1) s) (integers_from 2));;
val s : int stream = Cons (1, <lazy>)

sample s;;
- : int list = [1; 2; 2; 3; 3; 3; 4; 4; 4; 4; 5; 5]

Strumienie par - liczby Ramanujana

Skonstruujemy teraz strumień liczb Ramanujana. Liczby Ramanujana to takie liczby naturalne, które można przedstawić jako sumy sześcianów dwóch liczb naturalnych na dwa różne sposoby. Pierwsze z tych liczb, to: 1729, 4104, 13832, 20683.

Nasza konstrukcja będzie przebiegać następująco:

  • każdej parze liczb naturalnych przypiszemy wagę równą sumie ich sześcianów,
  • utworzymy strumień nieuporządkowanych par liczb naturalnych, uporządkowany wg. niemalejących wag par,
  • pary o równych wagach znajdą się obok siebie, co ułatwi wychwycenie takich wag, którym odpowiada więcej niż jedna para,
  • wynik będą stanowić takie powtarzające się wagi.

Pary nieuporządkowane będziemy reprezentować za pomocą par uporządkowanych, w których pierwsza współrzędna nie przekracza drugiej. Istotny będzie fakt, że wagi par rosną wraz ze zwiększeniem jednego z elementów pary. Wagi par definiujemy następująco:

let cube x = x * x * x;;
val cube : int -> int = <fun>

let weight (x, y) = cube x + cube y;;
val weight : int * int -> int = <fun>

Następująca procedura scala dwa różnowartościowe strumienie, uporządkowane wg. niemalejących wag elementów w jeden uporządkowany strumień różnowartościowy. Elementy powtarzające się są usuwane. Waga jest określona w postaci funkcji w. Bedziemy używać tej procedury do scalania par liczb naturalnych.

let rec merge_weighted w s1 s2 = 
  if empty s1 then s2 else
  if empty s2 then s1 else 
    let h1 = head s1 
    and h2 = head s2
    in 
      if w h1 < w h2 then 
        Cons (h1, merge_weighted w (tail s1) s2)
      else if w h1 > w h2 then 
        Cons (h2, merge_weighted w s1 (tail s2))
      else if h1 = h2 then 
        Cons (h1, merge_weighted w (tail s1) (tail s2))
      else 
        Cons (h1, Cons (h2, merge_weighted w (tail s1) (tail s2)));;
val merge_weighted : ('a -> 'b) -> 'a stream -> 'a stream -> 'a stream = <fun>

Niech s będzie rosnącym strumieniem liczb. Poniższa procedura tworzy strumień (nieuporządkowanych) par liczb z s, uporządkowany wg. niemalejących wag tych par. Istotne jest tu założenie, że wraz ze wzrostem jednej z liczb, waga pary rośnie.

Pierwszym elementem strumienia wynikowego jest oczywiście para (head s, head s). Pozostałe pary otrzymujemy scalając strumień par zawierających head s i jakąś inną liczbę, ze strumieniem par nie zawierających head s. Interesujący nas strumień par, to pairs nats.

let rec pairs s = 
  Cons ((head s, head s), 
        merge_weighted 
          weight
          (stream_map (function x -> (head s, x)) (tail s))
          (pairs (tail s)));;
val pairs : int stream -> (int * int) stream = <fun>

sample (pairs nats);;
- : (int * int) list = 
[(0, 0); (0, 1); (1, 1); (0, 2); (1, 2); (2, 2);  
 (0, 3); (1, 3); (2, 3); (3, 3); (0, 4); (1, 4)]

Poniższa procedura wychwytuje powtarzające się wagi par.

let rec non_uniq w s = 
  let rec skip we s =
    if empty s then Nil 
    else if we = w (head s) then skip we (tail s)
    else s
  in
    if empty s then Nil 
    else if empty (tail s) then Nil 
    else
      let h1 = head s
      and h2 = head (tail s)
      in 
        if w h1 = w h2 then 
          Cons (w h1, non_uniq w (skip (w h1) s))
        else
          non_uniq w (tail s);;
val non_uniq : ('a -> 'b) -> 'a stream -> 'b stream = <fun>

Strumień liczb Ramanujana możemy zdefiniować następująco:

let ramanujan = non_uniq weight (pairs nats);;
val ramanujan : int stream = Cons (1729, <lazy>)

sample ramanujan;;
- : int list = [1729; 4104; 13832; 20683; 32832; 39312; 40033; 46683; 64232; 65728; 110656; 110808]

Back-tracking + kontynuacje = strumień wyników

Pokażemy teraz jak łącząc technikę back-trackingu i kontynuacji możemy wygenerować wynik w postaci strumienia rozwiązań. Rozwiązanie to jest o tyle ciekawe, że poszukiwania kolejnych rozwiązań są prowadzone w miarę, jak ze strumienia wyników wyjmujemy kolejne elementy. Technikę tę pokażemy na dobrze znanym przykładzie problemu ośmiu hetmanów.

Predykat szach określa czy hetmany stojące na dwóch określonych pozycjach szachują się nawzajem.

let szach (x1, y1) (x2, y2) = x1=x2 or y1=y2 or x1+y1=x2+y2 or x1-y1=x2-y2;;
val szach : int * int -> int * int -> bool = <fun>

Predykat mozna_dostawic określa, czy do już stojących na szachownicy hetmanów można dostawić kolejnego.

exception Szach;;
exception Szach

let mozna_dostawic h het =
  try fold_left (fun ok x -> if szach h x then raise Szach else ok) true het
  with Szach -> false;;
val mozna_dostawic : int * int -> (int * int) list -> bool = <fun>

Właściwy back-tracking realizuje procedura dostaw, która przeszukuje wszystkie możliwe ustawienia hetmanów na planszy.

let hetman n =
  let rec dostaw het k l cont =
    if k = 0 then Cons(het, force cont) else
    if l = 0 then force cont else
    if mozna_dostawic (k, l) het then 
      dostaw ((k,l)::het) (k-1) n (lazy (dostaw het k (l-1) cont))
    else dostaw het k (l-1) cont 
  in dostaw [] n n (lazy Nil);;

Ćwiczenia

Ćwiczenia