Από τους Κύκλους στη Σιλικόνη: Υλοποίηση του Αλγορίθμου Κυκλικής ανάδυσης Πρώτων Αριθμών σε FPGA

Κυκλικός Αλγόριθμος Ανάδυσης Πρώτων

Ο αλγόριθμος βασίζεται στα εξής βήματα:

1. Ορισμός δύο ομόκεντρων κύκλων ακτίνας r και 2r.
2. Περιστροφή με σταθερή γραμμική ταχύτητα από κοινό σημείο εκκίνησης.
3. Κάθε φορά που ο πρώτος κύκλος ολοκληρώνει μια περιστροφή:
   1. Αυξάνεται ένας κεντρικός μετρητής n.
   2. Αν κανένας άλλος κύκλος δεν ολοκληρώνει περιστροφή τη δεδομένη χρονική στιγμή, η τιμή n του μετρητή καταγράφεται ως πρώτος αριθμός.
   3. Αν η τιμή n είναι πρώτος αριθμός, προστίθεται ένας νέος κύκλος ακτίνας nr.

H πλήρης μαθηματική απόδειξη του αλγόριθμου περιλαμβάνεται λεπτομερώς στο preprint paper της έρευνας.

Η Φιλοσοφία της Αρχιτεκτονικής: Μεταφορά σε Ψηφιακό Κύκλωμα

Σε περιβάλλον λογισμικού (software), αλγόριθμοι όπως το Κόσκινο του Ερατοσθένη απαιτούν αποθήκευση δεδομένων στη μνήμη για το μαρκάρισμα των στοιχείων, ενώ η δοκιμαστική διαίρεση (trial division) παρουσιάζει αυξανόμενη χρονική πολυπλοκότητα καθώς οι αριθμοί μεγαλώνουν.

Η σχεδίαση σε FPGA επιτρέπει την απευθείας αποτύπωση του αλγορίθμου σε ψηφιακή λογική, όπου τα βασικά γεωμετρικά στοιχεία μεταφράζονται ως εξής:

• Περιστρεφόμενοι Κύκλοι → Προγραμματιζόμενοι Oscillators: Κάθε πρώτος αριθμός p που εντοπίζεται, ανατίθεται σε έναν ανεξάρτητο μετρητή (oscillator) με περίοδο p. Στους μετρητές περιλαμβάνεται αρχικά μόνο ο 2, ενώ ο κύκλος του 1 ψηφιοποιείται από ένα γενικό counter.
• Γραμμή Αφετηρίας (Lap Line) → Ανίχνευση Μηδενισμού: Η συμπλήρωση p βημάτων από έναν μετρητή σηματοδοτεί την ολοκλήρωση μιας πλήρους περιστροφής.
• Έλεγχος Διαιρετότητας → Συνδυαστική Πύλη OR: Η εξέταση των σημάτων coincidence (at_lap) εκτελείται μέσω μιας πύλης OR.

Χρονική Πολυπλοκότητα O(1) ανά Βήμα: Εάν τουλάχιστον ένας ταλαντωτής ενεργοποιήσει το σήμα σύμπτωσης, η πύλη OR παράγει την έξοδο "1" (σήμανση σύνθετου αριθμού). Αυτή η συνδυαστική απόκριση εξαρτάται αποκλειστικά από την καθυστέρηση διάδοσης (propagation delay) του δέντρου των πυλών OR, ανεξάρτητα από το μέγεθος του n ή τον αριθμό των ενεργών ταλαντωτών. Συνεπώς, το σύστημα παρουσιάζει σταθερή καθυστέρηση αξιολόγησης O(1) ανά υποψήφιο αριθμό για μια δεδομένη υλοποίηση υλικού.

Μια αρχική προσέγγιση με πλήρως παράλληλη και ταυτόχρονη λειτουργία των oscillators απαιτεί δύο κύκλους ρολογιού ανά βήμα (έναν για την ενημέρωση και έναν για τον έλεγχο), διατηρώντας σταθερό χρόνο O(1).

Για τη βελτιστοποίηση της κατανάλωσης πόρων στο FPGA, επιλέχθηκε μια εναλλακτική αρχιτεκτονική ιεραρχικής οργάνωσης με τα εξής χαρακτηριστικά:

• Οι oscillators χωρίζονται σε clusters των 8 στοιχείων, όπου κάθε cluster ενημερώνει σειριακά τους oscillators του.
• Τα clusters λειτουργούν παράλληλα μεταξύ τους.
• Οι έξοδοι των clusters ομαδοποιούνται σε καταχωρητές (registers) των 8 bit.
• Οι καταχωρητές αυτοί οδηγούνται στην τελική πύλη OR.

Η συγκεκριμένη υλοποίηση ολοκληρώνει τον έλεγχο του συνόλου των oscillators σε 11 κύκλους ρολογιού, με τη χρονική πολυπλοκότητα να παραμένει O(1).

Πίνακας Προδιαγραφών για τον Εντοπισμό 1024 Πρώτων Αριθμών

Για την εύρεση των πρώτων 1024 πρώτων αριθμών, οι παράμετροι της αρχιτεκτονικής διαμορφώνονται ως εξής:

Αποτελέσματα σε ένα FPGA Cyclone IV E (EP4CE15E22C6)

Η αξιολόγηση του κυκλώματος πραγματοποιήθηκε μέσω του Quartus Prime σε FPGA Cyclone IV E (EP4CE15E22C6) για τον εντοπισμό των πρώτων 1024 πρώτων αριθμών (P = 1024, N_BITS = 13). Τα δεδομένα από τα reports σύνθεσης (Flow Summary) και χρονισμού (TimeQuest Timing Analyzer) συνοψίζονται παρακάτω:

Η ανάλυση του κρίσιμου μονοπατιού (Worst-Case Critical Path) έδωσε θετικό slack 3.303 ns, με Data Delay 6.618 ns. Το μονοπάτι εντοπίζεται στη διαδρομή από το bus φόρτωσης load_val[8] προς τον καταχωρητή counters[0][8] ενός εκ των clusters, λόγω του φαινομένου fan-out κατά τη διακλάδωση του σήματος στα 128 clusters.

Η αρχιτεκτονική καταλαμβάνει το 67% των διαθέσιμων LEs ενός entry-level FPGA παλαιότερης γενιάς. Σε σύγχρονες διατάξεις FPGA (π.χ. τεχνολογίας 28nm), η συχνότητα λειτουργίας αναμένεται να αυξηθεί βάσει των προδιαγραφών της σιλικόνης, χωρίς τροποποίηση του κώδικα VHDL.

Απόλυτη Χρονική Προβλεψιμότητα (Deterministic Timing)

Σε περιβάλλον λογισμικού, ο χρόνος εκτέλεσης παρουσιάζει διακυμάνσεις λόγω cache misses, διακλαδώσεων κώδικα και προγραμματισμού του λειτουργικού συστήματος. Στην εξεταζόμενη hardware αρχιτεκτονική, ο χρόνος εκτέλεσης είναι απολύτως ντετερμινιστικός.

Καθώς το σύστημα αξιολογεί κάθε υποψήφιο αριθμό n σε ένα σταθερό παράθυρο 11 κύκλων ρολογιού (8 κύκλοι για τη σειριακή ενημέρωση των slots και 3 κύκλοι για το pipeline της πύλης OR και τη λήψη απόφασης), ο απαιτούμενος αριθμός κύκλων για την εύρεση του P-οστού πρώτου αριθμού υπολογίζεται με ακρίβεια από τη σχέση:

(όπου p_P είναι η αριθμητική αξία του P-οστού πρώτου αριθμού).

Οι cycle-accurate behavioral προσομοιώσεις στο περιβάλλον ModelSim επαλήθευσαν την παραπάνω εξίσωση:

• Για P = 256 (p256 = 1619): Απαιτήθηκαν ακριβώς 17.789 κύκλοι ρολογιού
• Για P = 1024 (p1024 = 8161): Η αναζήτηση ολοκληρώθηκε σε ακριβώς 89.751 κύκλους ρολογιού.

Ο σταθερός αυτός χρονισμός καθιστά την αρχιτεκτονική κατάλληλη για συστήματα σκληρού πραγματικού χρόνου (Hard Real-Time Systems), όπου απαιτείται εγγυημένο άνω όριο καθυστέρησης (bounded latency).

Μελλοντικές Βελτιώσεις & Οριζόντια Κλιμάκωση (Scalability)

Η δομή των clusters επιτρέπει την εφαρμογή συγκεκριμένων βελτιστοποιήσεων:

• Μείωση του CHECK_PERIOD: Η τρέχουσα υλοποίηση του FSM απαιτεί 11 κύκλους για τη σάρωση. Με την τεχνική του interleaving (επικάλυψη ελέγχων), ο χρόνος αυτός δύναται να μειωθεί στους 8 κύκλους ρολογιού, επιφέροντας μείωση της χρονικής καθυστέρησης κατά ~27%.
• Wheel Factorization (Παράκαμψη Υποψηφίων): Ο κεντρικός μετρητής αυξάνεται κατά 1 σε κάθε βήμα (counter + 1). Τροποποιώντας τον ελεγκτή ώστε να παρακάμπτει τα πολλαπλασία του 2 και του 3), ο αριθμός των προς εξέταση ψηφίων μειώνεται στο 33%, με ανάλογη επίδραση στην ταχύτητα εύρεσης.
• Μεταφορά σε Block RAM (BRAM): Για μεγάλες τιμές του P (π.χ. P > 10.000), οι πίνακες counters και primes των clusters μπορούν να υλοποιηθούν στις ενσωματωμένες μνήμες (M9K) του FPGA αντί για registers. Η σειριακή προσπέλαση των slots (μέσω του σήματος slot) είναι συμβατή με τη λειτουργία dual-port των BRAMs, περιορίζοντας τη χρήση Logic Elements.

Πλεονεκτήματα & Μειονεκτήματα

Πλεονεκτήματα

• Επαλήθευση O(1) ανά Βήμα: Το συνδυαστικό δέντρο μείωσης της πύλης OR αποκρίνεται σε σταθερό χρόνο ανά υποψήφιο αριθμό, πλήρως ανεξάρτητα από το εύρος bit του n ή το πλήθος των ενεργών ταλαντωτών. Αυτή η συνδυαστική απόκριση καθορίζεται από την καθυστέρηση διάδοσης του δέντρου μείωσης (OR-reduction tree), το οποίο περιλαμβάνει βαθμίδες pipeline. Για οποιαδήποτε σταθερή παραμετροποίηση P ενεργών ταλαντωτών, αυτή η καθυστέρηση αξιολόγησης απαιτεί έναν αυστηρά σταθερό αριθμό κύκλων ρολογιού ανά υποψήφιο αριθμό, εξαλείφοντας τις διακυμάνσεις χρόνου και εξασφαλίζοντας ένα προσδιοριστικό προφίλ πολυπλοκότητας O(1) ανά βήμα.
• Απουσία χρήσης RAM: Δεν απαιτείται πίνακας καταγραφής σύνθετων αριθμών, καθώς η κατάσταση αποθηκεύεται αποκλειστικά στα flip-flops των oscillators.
• Παράλληλη λειτουργία: Τα clusters των oscillators λειτουργούν ταυτόχρονα και ανεξάρτητα, χωρίς sequential εξαρτήσεις μεταξύ τους.
• Deterministic timing: Ο ακριβής χρόνος εκτέλεσης είναι προβλέψιμος μέσω της μαθηματικής σχέσης 11 * (pP - 2) + 2.
• Δυνατότητα κλιμάκωσης: Η αρχιτεκτονική επιτρέπει την επέκταση μέσω daisy-chain σύνδεσης πολλαπλών FPGA.

Μειονεκτήματα

• Προκαθορισμένο ανώτατο όριο P: Το μέγιστο πλήθος των προς εύρεση πρώτων αριθμών περιορίζεται από τους διαθέσιμους φυσικούς πόρους του FPGA.
• Γραμμική εξάρτηση συνολικού χρόνου: Ο συνολικός χρόνος εκτέλεσης κλιμακώνεται ως O(p_P), εξαρτώμενος από την τιμή του τελευταίου πρώτου αριθμού και όχι από το πλήθος P.
• Γραμμική αύξηση πόρων: Κάθε oscillator απαιτεί σταθερό αριθμό LEs, καθιστώντας αναγκαία τη χρήση BRAM ή SRL32 για μεγάλες τιμές του P.
• Simulation overhead: Η προσομοίωση για υψηλές τιμές του P παρουσιάζει καθυστερήσεις λόγω του μεγάλου αριθμού παράλληλων διεργασιών στο λογισμικό.

Μελλοντική Έρευνα & Συγκριτική Αξιολόγηση (Benchmarking)

Σε επόμενο άρθρο, θα παρουσιάσουμε μια ολοκληρωμένη συγκριτική ανάλυση, αξιολογώντας τη σχεδίαση Cluster-Oscillator που παρουσιάστηκε εδώ έναντι των τριών κυρίαρχων αρχιτεκτονικών υλικού της βιβλιογραφίας: των υλοποιήσεων Δοκιμαστικής Διαίρεσης σε Υλικό (Hardware Trial Division), των Κοσκίνων βασισμένων σε Block RAM (BRAM-Based Sieves), καθώς και των τρεχουσών λύσεων αναφοράς (State-of-the-Art - SOTA) που χρησιμοποιούν Συστολικούς Πίνακες (Pipelined Systolic Arrays). Ιδιαίτερη έμφαση θα δοθεί στη χαρτογράφηση των ξεκάθαρων αρχιτεκτονικών ισοζυγίων (tradeoffs) ανάμεσα στον δικό μας βρόχο αξιολόγησης σταθερού χρόνου O(1) και στα όρια υψηλής συχνότητας των συστολικών πινάκων, ως προς τη ρυθμαπόδοση, την κλιμακώσιμότητα και την κατανάλωση πόρων πυριτίου.

Ο κώδικας της υλοποίησης

Ακολουθεί ο πηγαίος κώδικας σε VHDL. Το σύστημα αποτελείται από το αρχείο oscillators_cluster.vhd (υλοποίηση των 8 σειριακά ανανεώσιμων oscillators ανά cluster) και το αρχείο primes_finder.vhd (κεντρικός ελεγκτής, Pipeline Registers και FSM).

Οι παράμετροι έχουν οριστεί για την εύρεση P = 1024 πρώτων αριθμών με εύρος διαύλου N_BITS = 13. Καθώς ο 1024ος πρώτος αριθμός είναι το 8161, η επιλογή των 13 bits καλύπτει το απαιτούμενο εύρος, δεδομένου ότι 2^12 < 8161 < 2^13.

Υλοποίηση του Component oscillators_cluster

library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.NUMERIC_STD.ALL;

entity oscillators_cluster is
	generic (
		BIT_WIDTH    : integer := 10
	);
	port (
		clk          : in  std_logic;
		reset        : in  std_logic;
		run_en       : in  std_logic;
		load_en      : in  std_logic;
		load_slot    : in  std_logic_vector(2 downto 0);
		load_val     : in  std_logic_vector(BIT_WIDTH-1 downto 0);
		coincidence  : out std_logic
	);
end oscillators_cluster;

architecture Behavioral of oscillators_cluster is

	type count_array is array (0 to 7) of unsigned(BIT_WIDTH-1 downto 0);

	signal counters	: count_array := (others => to_unsigned(2, BIT_WIDTH));
	signal primes 	: count_array := (others => to_unsigned(2, BIT_WIDTH));
	signal slot : unsigned(2 downto 0) := (others => '0');
	signal slot_delayed : integer range 0 to 7 := 0;
	signal pipe_counter : unsigned(BIT_WIDTH-1 downto 0) := to_unsigned(2, BIT_WIDTH);
	signal pipe_prime   : unsigned(BIT_WIDTH-1 downto 0) := to_unsigned(2, BIT_WIDTH);
 
begin

	process(clk)
		variable load_slot_idx    : integer;
	begin
		if rising_edge(clk) then
			if reset = '1' then
				 slot         <= (others => '0');
				 slot_delayed <= 0;
				 pipe_counter <= to_unsigned(2, BIT_WIDTH);
				 pipe_prime   <= to_unsigned(2, BIT_WIDTH);
				 coincidence  <= '0';
			else
				load_slot_idx    := to_integer(unsigned(load_slot));

				if run_en = '1' then
					if slot_delayed = 0 then
						if pipe_counter = 1 then
							counters(slot_delayed) <= pipe_prime;
							coincidence <= '1';
						else
							counters(slot_delayed) <= pipe_counter - 1;
							coincidence <= '0'; -- Καθαρισμός στην αρχή του frame
						end if;
					else
						if pipe_counter = 1 then
							counters(slot_delayed) <= pipe_prime;
							coincidence <= '1';
						else
							counters(slot_delayed) <= pipe_counter - 1;
						end if;
					end if;

					pipe_counter <= counters(to_integer(slot));
					pipe_prime   <= primes(to_integer(slot));
					slot_delayed <= to_integer(slot); -- Αποθήκευση της θέσης για τον επόμενο κύκλο
					slot <= slot + 1;
				  
				end if;

				if load_en = '1' then
					primes(load_slot_idx) <= unsigned(load_val);
					counters(load_slot_idx)     <= unsigned(load_val);
				end if;

			end if;
		end if;
	end process;

end Behavioral;

Υλοποίηση του Κεντρικού Controller

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;

entity primes_finder is
    generic (
        P      : integer := 1024;
        N_BITS : integer := 13
    );
    port (
        clk         : in  std_logic;
        rst         : in  std_logic;
        counter_out : out unsigned(N_BITS-1 downto 0);
        prime_out   : out unsigned(N_BITS-1 downto 0);
        done_out        : out std_logic
    );
end entity primes_finder;

architecture rtl of primes_finder is

	constant CHECK_PERIOD : integer := 11;
	constant NUM_CLUSTERS : integer := (P + 7) / 8;
	constant NUM_GROUPS   : integer := (NUM_CLUSTERS + 7) / 8;

	signal counter					: unsigned(N_BITS-1 downto 0) := (others => '0');
	signal timer        			: integer range 0 to CHECK_PERIOD := 0; 
	signal primes_count 			: integer range 0 to P := 0;
	signal prime_found  			: unsigned(N_BITS-1 downto 0);
	signal any_lap           	: std_logic; 
	signal state_initialized 	: std_logic := '0';
	signal done 					: std_logic := '0'; 
			
	signal run_en				: std_logic;
	signal load_en_all		: std_logic := '0';
	signal load_en_vector	: std_logic_vector(NUM_CLUSTERS-1 downto 0) := (others => '0');
	signal load_slot			: std_logic_vector(2 downto 0);
	signal load_val			: std_logic_vector(N_BITS-1 downto 0);

	signal all_coincidences  : std_logic_vector(NUM_CLUSTERS - 1 downto 0) := (others => '0');
	
	type group_regs_t is array (0 to NUM_GROUPS - 1) of std_logic_vector(7 downto 0);
   signal reg_groups : group_regs_t := (others => (others => '0'));
   signal groups_has_lap : std_logic_vector(NUM_GROUPS - 1 downto 0) := (others => '0');
	 
	type state_t is (RUNNING, FINISHED);
	signal state : state_t := RUNNING;
 
begin 

	counter_out	<= counter;
	prime_out		<= prime_found;
	done_out		<= done;

	run_en <= '1' when (state = RUNNING and timer < 8) else '0';

	GEN_CLUSTERS: for i in 0 to NUM_CLUSTERS - 1 generate
		cluster_inst : entity work.oscillators_cluster
			generic map (
				BIT_WIDTH => N_BITS
			)
			port map (
				clk          => clk,
				reset        => rst,
				run_en       => run_en,
				load_en      => load_en_vector(i),
				load_slot    => load_slot,
				load_val     => load_val,
				coincidence  => all_coincidences(i)
			);
	end generate;
 
	process(clk, rst)
	begin
		if rst = '1' then
			reg_groups <= (others => (others => '0'));

		elsif rising_edge(clk) then
			reg_groups <= (others => (others => '0'));
			for i in 0 to NUM_CLUSTERS - 1 loop
				reg_groups(i / 8)(i rem 8) <= all_coincidences(i);
			end loop;
		end if;
	end process;

	GEN_GROUP_DETECTORS: for g in 0 to NUM_GROUPS - 1 generate
		groups_has_lap(g) <= '1' when reg_groups(g) /= "00000000" else '0';
	end generate;

	any_lap <= '0' when groups_has_lap = (groups_has_lap'range => '0') else '1';
	 
	process(clk, rst)
	begin
		if rst = '1' then
			counter           <= to_unsigned(2, N_BITS); 
			prime_found       <= to_unsigned(2, N_BITS);
			primes_count      <= 1;         
			done					<= '0';
			timer             <= 0;
			load_en_vector    <= (others => '0');
			state             <= RUNNING;

		elsif rising_edge(clk) then
			load_en_vector <= (others => '0');
		  
			case state is
				when RUNNING =>
			 
					if timer = 0 then
						if done /= '1' then
							  counter <= counter + 1;
						end if;
					end if;
			 
					if timer = CHECK_PERIOD - 1 then
						if any_lap = '0' then
						 
							prime_found  <= counter;
							primes_count <= primes_count + 1;

							load_en_vector(primes_count / 8) <= '1'; 
							load_slot <= std_logic_vector(to_unsigned(primes_count rem 8, 3));
							load_val  <= std_logic_vector(counter);

							if (primes_count + 1 >= P) then
								done  <= '1';
								state <= FINISHED;
							end if;
						end if;
						
						timer <= 0;
					else
						timer <= timer + 1;
					end if;

				when FINISHED =>
					null;
				end case;
			end if;
	end process;

end architecture rtl;

Open Science

Πιστεύοντας στη φιλοσοφία του open-source και της ελεύθερης διακίνησης της γνώσης, ολόκληρο το project είναι προσβάσιμο από την κοινότητα:

• Διαβάστε το Preprint Paper στο Zenodo: DOI: 10.5281/zenodo.20622711
• Κατεβάστε τον Κώδικα και τα Testbenches στο GitHub: GitHub Repository

Πρόκειται για την πρώτη hardware υλοποίηση του Αλγορίθμου Κυκλικής Ανάδυσης Πρώτων Αριθμών. Ο κώδικας είναι πλήρως παραμετροποιήσιμος και ανοιχτός για πειραματισμό — είτε σε simulator είτε σε πραγματικό FPGA. Κάθε νέος πρώτος που βρίσκεται από το σύστημα αποτελεί και εμπειρική επαλήθευση της θεωρίας.

Θα χαρώ πολύ να ακούσω τις σκέψεις σας, τις προτάσεις σας για βελτιώσεις στον κώδικα VHDL ή τις δικές σας δοκιμές σε περιβάλλοντα προσομοίωσης.