sw: Add sample Monte Carlo application

sw/apps/montecarlo/pi_estimation/main.c

@@ -0,0 +1,55 @@
+// Copyright 2023 ETH Zurich and University of Bologna.
+// Licensed under the Apache License, Version 2.0, see LICENSE for details.
+// SPDX-License-Identifier: Apache-2.0
+
+#include "lcg.h"
+#include "math.h"
+#include "pi_estimation.h"
+#include "snrt.h"
+
+#define N_SAMPLES 1024
+
+__thread uint32_t seed0, seed1, Ap, Cp;
+
+double pi_estimate;
+
+inline void mc_init() {
+    // Double the sequences as each core produces two random numbers
+    unsigned int num_sequences = 2 * snrt_cluster_compute_core_num();
+    init_2d_lcg_params(num_sequences, 0, LCG_A, LCG_C, &seed0, &seed1, &Ap,
+                       &Cp);
+}
+
+int main() {
+    // Initialize the PRNGs for parallel Monte Carlo
+    if (snrt_is_compute_core()) mc_init();
+
+    // Store partial sum array at first free address in TCDM
+    uint32_t* reduction_array = (uint32_t*)snrt_l1_next();
+
+    // Calculate partial sums
+    uint32_t n_samples_per_core = N_SAMPLES / snrt_cluster_compute_core_num();
+    reduction_array[snrt_cluster_core_idx()] =
+        calculate_partial_sum(seed0, seed1, Ap, Cp, n_samples_per_core);
+
+    // Synchronize cores
+    snrt_cluster_hw_barrier();
+
+    // First core in cluster performs the final calculation
+    if (snrt_cluster_core_idx() == 0) {
+        // Reduce partial sums
+        uint32_t sum = 0;
+        for (int i = 0; i < snrt_cluster_compute_core_num(); i++) {
+            sum += reduction_array[i];
+        }
+
+        // Estimate pi
+        pi_estimate = estimate_pi(sum, N_SAMPLES);
+
+        // Check result
+        double err = fabs(pi_estimate - M_PI);
+        if (err > 0.05) return 1;
+    }
+
+    return 0;
+}

sw/apps/montecarlo/pi_estimation/pi_estimation.h

@@ -0,0 +1,30 @@
+// Copyright 2023 ETH Zurich and University of Bologna.
+// Licensed under the Apache License, Version 2.0, see LICENSE for details.
+// SPDX-License-Identifier: Apache-2.0
+
+__thread double one = 1.0;
+
+inline uint32_t calculate_partial_sum(uint32_t seed0, uint32_t seed1,
+                                      uint32_t Ap, uint32_t Cp,
+                                      unsigned int n_samples) {
+    uint32_t x1 = seed0;
+    uint32_t x2 = seed1;
+    double u1, u2;
+    unsigned int result = 0;
+
+    for (unsigned int i = 0; i < n_samples; i++) {
+        u1 = normalize(x1);
+        u2 = normalize(x2);
+        x1 = lcg(Ap, Cp, x1);
+        x2 = lcg(Ap, Cp, x2);
+
+        if ((u1 * u1 + u2 * u2) < one) {
+            result++;
+        }
+    }
+    return result;
+}
+
+inline double estimate_pi(uint32_t sum, uint32_t n_samples) {
+    return (double)(4 * sum) / (double)n_samples;
+}

sw/apps/prng/lcg.h

@@ -0,0 +1,59 @@
+// Copyright 2023 ETH Zurich and University of Bologna.
+// Licensed under the Apache License, Version 2.0, see LICENSE for details.
+// SPDX-License-Identifier: Apache-2.0
+
+#include "snrt.h"
+
+// Numerical Recipes from the "quick and dirty generators" list, Chapter 7.1,
+// Eq. 7.1.6 parameters from Knuth and H. W. Lewis
+#define MAX_UINT_PLUS1 4294967296.0
+#define LCG_A 1664525
+#define LCG_C 1013904223
+
+__thread double max_uint_plus_1_inverse = (double)1.0 / (double)MAX_UINT_PLUS1;
+
+// Calculate A' and C' constants
+inline void leapfrog_constants(unsigned int num_sequences, uint32_t a,
+                               uint32_t c, uint32_t* Ap, uint32_t* Cp) {
+    // Ap = a^k
+    // Cp = (a^(k-1) + a^(k-2) + ... + a^1 + a^0)*c
+    uint32_t Ap_tmp = 1;
+    uint32_t Cp_tmp = 0;
+    for (unsigned int p = 0; p < num_sequences; p++) {
+        Cp_tmp += Ap_tmp;
+        Ap_tmp *= a;
+    }
+    Cp_tmp *= c;
+
+    // Store temporary variables to outputs
+    *Ap = Ap_tmp;
+    *Cp = Cp_tmp;
+}
+
+// Calculate seed for leapfrog method's right sequence of index `sequence_idx`
+inline uint32_t right_seed(uint32_t left_seed, uint32_t a, uint32_t c,
+                           unsigned int sequence_idx) {
+    uint32_t seed = left_seed;
+    for (unsigned int p = 1; p <= sequence_idx; p++) {
+        seed = seed * a + c;
+    }
+    return seed;
+}
+
+// Generate next PRN from LCG recurrence equation
+inline uint32_t lcg(uint32_t a, uint32_t c, uint32_t previous) {
+    return previous * a + c;
+}
+
+// Normalize LCG PRN to [0, 1) range
+inline double normalize(uint32_t x) {
+    return (double)x * max_uint_plus_1_inverse;
+}
+
+inline void init_2d_lcg_params(unsigned int num_sequences, uint32_t seed,
+                               uint32_t a, uint32_t c, uint32_t* seed0,
+                               uint32_t* seed1, uint32_t* Ap, uint32_t* Cp) {
+    *seed0 = right_seed(seed, a, c, 2 * snrt_global_compute_core_idx());
+    *seed1 = right_seed(seed, a, c, 2 * snrt_global_compute_core_idx() + 1);
+    leapfrog_constants(num_sequences, a, c, Ap, Cp);
+}

target/snitch_cluster/sw/apps/Makefile

@@ -17,6 +17,7 @@ SUBDIRS += dnn/layernorm
 SUBDIRS += dnn/linear
 SUBDIRS += dnn/maxpool
 SUBDIRS += dnn/softmax
+SUBDIRS += montecarlo/pi_estimation
 
 .PHONY: all clean $(SUBDIRS)
 

target/snitch_cluster/sw/apps/montecarlo/pi_estimation/Makefile

@@ -0,0 +1,14 @@
+# Copyright 2023 ETH Zurich and University of Bologna.
+# Licensed under the Apache License, Version 2.0, see LICENSE for details.
+# SPDX-License-Identifier: Apache-2.0
+#
+# Luca Colagrande <[email protected]>
+
+PI_ESTIMATION_DIR = ../../../../../../sw/apps/montecarlo/pi_estimation
+PRNG_DIR          = ../../../../../../sw/apps/prng
+
+APP     ?= pi_estimation
+SRCS    ?= $(PI_ESTIMATION_DIR)/main.c
+INCDIRS += $(PI_ESTIMATION_DIR) $(PRNG_DIR)
+
+include ../../common.mk

target/snitch_cluster/sw/run.yaml

@@ -84,3 +84,4 @@ runs:
   #   throws illegal instruction on FDIV in simulation
   # - elf: apps/dnn/softmax/build/softmax.elf
   #   throws illegal instruction on FDIV in simulation
+  - elf: apps/montecarlo/pi_estimation/build/pi_estimation.elf