Update rust.yml

src/bindings/python.rs

@@ -206,6 +206,9 @@ struct PyRunArgs {
     /// bool: Should the circuit use range checks for inputs and outputs (set to false if the input is a felt)
     #[pyo3(get, set)]
     pub ignore_range_check_inputs_outputs: bool,
+    /// float: epsilon used for arguments that use division
+    #[pyo3(get, set)]
+    pub epsilon: f64,
 }
 
 /// default instantiation of PyRunArgs
@@ -238,12 +241,14 @@ impl From<PyRunArgs> for RunArgs {
             decomp_base: py_run_args.decomp_base,
             decomp_legs: py_run_args.decomp_legs,
             ignore_range_check_inputs_outputs: py_run_args.ignore_range_check_inputs_outputs,
+            epsilon: Some(py_run_args.epsilon),
         }
     }
 }
 
 impl Into<PyRunArgs> for RunArgs {
     fn into(self) -> PyRunArgs {
+        let eps = self.get_epsilon();
         PyRunArgs {
             bounded_log_lookup: self.bounded_log_lookup,
             input_scale: self.input_scale,
@@ -262,6 +267,7 @@ impl Into<PyRunArgs> for RunArgs {
             decomp_base: self.decomp_base,
             decomp_legs: self.decomp_legs,
             ignore_range_check_inputs_outputs: self.ignore_range_check_inputs_outputs,
+            epsilon: eps,
         }
     }
 }

src/circuit/ops/hybrid.rs

@@ -1,7 +1,7 @@
 use super::*;
 use crate::{
     circuit::{layouts, utils},
-    fieldutils::{integer_rep_to_felt, IntegerRep},
+    fieldutils::{IntegerRep, integer_rep_to_felt},
     graph::multiplier_to_scale,
     tensor::{self, DataFormat, Tensor, TensorType, ValTensor},
 };
@@ -15,10 +15,12 @@ use serde::{Deserialize, Serialize};
 pub enum HybridOp {
     Ln {
         scale: utils::F32,
+        eps: f64,
     },
     Rsqrt {
         input_scale: utils::F32,
         output_scale: utils::F32,
+        eps: f64,
     },
     Sqrt {
         scale: utils::F32,
@@ -42,6 +44,7 @@ pub enum HybridOp {
     Recip {
         input_scale: utils::F32,
         output_scale: utils::F32,
+        eps: f64,
     },
     Div {
         denom: utils::F32,
@@ -77,6 +80,7 @@ pub enum HybridOp {
         input_scale: utils::F32,
         output_scale: utils::F32,
         axes: Vec<usize>,
+        eps: f64,
     },
     Output {
         decomp: bool,
@@ -128,12 +132,13 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
             HybridOp::Rsqrt {
                 input_scale,
                 output_scale,
+                eps,
             } => format!(
-                "RSQRT (input_scale={}, output_scale={})",
-                input_scale, output_scale
+                "RSQRT (input_scale={}, output_scale={}, eps={})",
+                input_scale, output_scale, eps
             ),
             HybridOp::Sqrt { scale } => format!("SQRT(scale={})", scale),
-            HybridOp::Ln { scale } => format!("LN(scale={})", scale),
+            HybridOp::Ln { scale, eps } => format!("LN(scale={}, eps={})", scale, eps),
             HybridOp::RoundHalfToEven { scale, legs } => {
                 format!("ROUND_HALF_TO_EVEN(scale={}, legs={})", scale, legs)
             }
@@ -146,16 +151,18 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
             HybridOp::Recip {
                 input_scale,
                 output_scale,
+                eps,
             } => format!(
-                "RECIP (input_scale={}, output_scale={})",
-                input_scale, output_scale
+                "RECIP (input_scale={}, output_scale={}, eps={})",
+                input_scale, output_scale, eps
             ),
             HybridOp::Div { denom } => format!("DIV (denom={})", denom),
             HybridOp::SumPool {
                 padding,
                 stride,
                 kernel_shape,
-                normalized, data_format
+                normalized,
+                data_format,
             } => format!(
                 "SUMPOOL (padding={:?}, stride={:?}, kernel_shape={:?}, normalized={}, data_format={:?})",
                 padding, stride, kernel_shape, normalized, data_format
@@ -177,10 +184,11 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
                 input_scale,
                 output_scale,
                 axes,
+                eps,
             } => {
                 format!(
-                    "SOFTMAX (input_scale={}, output_scale={}, axes={:?})",
-                    input_scale, output_scale, axes
+                    "SOFTMAX (input_scale={}, output_scale={}, axes={:?}, eps={})",
+                    input_scale, output_scale, axes, eps
                 )
             }
             HybridOp::Output { decomp } => {
@@ -211,17 +219,21 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
             HybridOp::Rsqrt {
                 input_scale,
                 output_scale,
+                eps,
             } => layouts::rsqrt(
                 config,
                 region,
                 values[..].try_into()?,
                 *input_scale,
                 *output_scale,
+                *eps,
             )?,
             HybridOp::Sqrt { scale } => {
                 layouts::sqrt(config, region, values[..].try_into()?, *scale)?
             }
-            HybridOp::Ln { scale } => layouts::ln(config, region, values[..].try_into()?, *scale)?,
+            HybridOp::Ln { scale, eps } => {
+                layouts::ln(config, region, values[..].try_into()?, *scale, *eps)?
+            }
             HybridOp::RoundHalfToEven { scale, legs } => {
                 layouts::round_half_to_even(config, region, values[..].try_into()?, *scale, *legs)?
             }
@@ -255,12 +267,14 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
             HybridOp::Recip {
                 input_scale,
                 output_scale,
+                eps,
             } => layouts::recip(
                 config,
                 region,
                 values[..].try_into()?,
                 integer_rep_to_felt(input_scale.0 as IntegerRep),
                 integer_rep_to_felt(output_scale.0 as IntegerRep),
+                *eps,
             )?,
             HybridOp::Div { denom, .. } => {
                 if denom.0.fract() == 0.0 {
@@ -317,6 +331,7 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
                 input_scale,
                 output_scale,
                 axes,
+                eps,
             } => layouts::softmax_axes(
                 config,
                 region,
@@ -324,6 +339,7 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
                 *input_scale,
                 *output_scale,
                 axes,
+                *eps,
             )?,
             HybridOp::Output { decomp } => {
                 layouts::output(config, region, values[..].try_into()?, *decomp)?
@@ -364,6 +380,7 @@ impl<F: PrimeField + TensorType + PartialOrd + std::hash::Hash> Op<F> for Hybrid
             } => multiplier_to_scale((output_scale.0 * input_scale.0) as f64),
             HybridOp::Ln {
                 scale: output_scale,
+                eps: _,
             } => 4 * multiplier_to_scale(output_scale.0 as f64),
             _ => in_scales[0],
         };

src/circuit/ops/layouts.rs

@@ -303,6 +303,7 @@ pub(crate) fn recip<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
     value: &[ValTensor<F>; 1],
     input_scale: F,
     output_scale: F,
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     let input = value[0].clone();
     let input_dims = input.dims();
@@ -317,6 +318,7 @@ pub(crate) fn recip<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
             &input_evals,
             felt_to_integer_rep(input_scale) as f64,
             felt_to_integer_rep(output_scale) as f64,
+            eps,
         )
         .par_iter()
         .map(|x| Value::known(integer_rep_to_felt(*x)))
@@ -335,7 +337,7 @@ pub(crate) fn recip<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
     let claimed_output = identity(config, region, &[claimed_output], true)?;
     // divide by input_scale
     let zero_inverse_val =
-        tensor::ops::nonlinearities::zero_recip(felt_to_integer_rep(output_scale) as f64)[0];
+        tensor::ops::nonlinearities::zero_recip(felt_to_integer_rep(output_scale) as f64, eps)[0];
     let zero_inverse = create_constant_tensor(integer_rep_to_felt(zero_inverse_val), 1);
 
     let equal_zero_mask = equals_zero(config, region, &[input.clone()])?;
@@ -473,7 +475,7 @@ pub fn sqrt<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
 ///    Some(&[1, 2, 3, 2, 3, 4, 3, 4, 5]),
 /// &[3, 3],
 /// ).unwrap());
-/// let result = rsqrt::<Fp>(&dummy_config, &mut dummy_region, &[x], 1.0.into(), 1.0.into()).unwrap();
+/// let result = rsqrt::<Fp>(&dummy_config, &mut dummy_region, &[x], 1.0.into(), 1.0.into(), f64::EPSILON).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 1, 1, 1, 1, 1, 1, 1, 1]), &[3, 3]).unwrap();
 /// assert_eq!(result.int_evals().unwrap(), expected);
 /// ```
@@ -483,13 +485,21 @@ pub fn rsqrt<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
     value: &[ValTensor<F>; 1],
     input_scale: utils::F32,
     output_scale: utils::F32,
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     let sqrt = sqrt(config, region, value, input_scale)?;
 
     let felt_output_scale = integer_rep_to_felt(output_scale.0 as IntegerRep);
     let felt_input_scale = integer_rep_to_felt(input_scale.0 as IntegerRep);
 
-    let recip = recip(config, region, &[sqrt], felt_input_scale, felt_output_scale)?;
+    let recip = recip(
+        config,
+        region,
+        &[sqrt],
+        felt_input_scale,
+        felt_output_scale,
+        eps,
+    )?;
 
     Ok(recip)
 }
@@ -1547,7 +1557,7 @@ pub(crate) fn dynamic_lookup<F: PrimeField + TensorType + PartialOrd + std::hash
 ///      * Creates pairs: (index_input, value_input) for original elements
 ///      * Creates pairs: (index_output, value_output) for permuted elements
 ///      * index_input is a fixed sequence 0,1,2... corresponding to input positions
-///    
+///
 ///    - Core permutation verification:
 ///      * For each (index_input, value_input), verify there exists exactly one
 ///        (index_output, value_output) such that value_input = value_output
@@ -5702,7 +5712,7 @@ pub fn ceil<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
 /// &[1, 1, 2, 2],
 /// ).unwrap());
 ///
-/// let result = ln::<Fp>(&dummy_config, &mut dummy_region, &[x], 2.0.into()).unwrap();
+/// let result = ln::<Fp>(&dummy_config, &mut dummy_region, &[x], 2.0.into(), f64::EPSILON).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[4, 0, 4, -8]), &[1, 1, 2, 2]).unwrap();
 /// assert_eq!(result.int_evals().unwrap(), expected);
 ///
@@ -5712,6 +5722,7 @@ pub fn ln<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
     region: &mut RegionCtx<F>,
     values: &[ValTensor<F>; 1],
     scale: utils::F32,
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     // first generate the claimed val
 
@@ -5882,6 +5893,7 @@ pub fn ln<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
         &[pow2_prior_to_claimed_distance],
         scale_as_felt,
         scale_as_felt * scale_as_felt,
+        eps,
     )?;
 
     let interpolated_distance = pairwise(
@@ -5910,6 +5922,7 @@ pub fn ln<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
         &[pow2_next_to_claimed_distance],
         scale_as_felt,
         scale_as_felt * scale_as_felt,
+        eps,
     )?;
 
     let interpolated_distance_next = pairwise(
@@ -6698,12 +6711,13 @@ pub(crate) fn softmax_axes<F: PrimeField + TensorType + PartialOrd + std::hash::
     input_scale: utils::F32,
     output_scale: utils::F32,
     axes: &[usize],
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     let soft_max_at_scale = move |config: &BaseConfig<F>,
                                   region: &mut RegionCtx<F>,
                                   values: &[ValTensor<F>; 1]|
           -> Result<ValTensor<F>, CircuitError> {
-        softmax(config, region, values, input_scale, output_scale)
+        softmax(config, region, values, input_scale, output_scale, eps)
     };
 
     let output = multi_dim_axes_op(config, region, values, axes, soft_max_at_scale)?;
@@ -6718,6 +6732,7 @@ pub(crate) fn percent<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>
     values: &[ValTensor<F>; 1],
     input_scale: utils::F32,
     output_scale: utils::F32,
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     let is_assigned = values[0].all_prev_assigned();
     let mut input = values[0].clone();
@@ -6736,6 +6751,7 @@ pub(crate) fn percent<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>
         &[denom],
         input_felt_scale,
         output_felt_scale,
+        eps,
     )?;
     // product of num * (1 / denom) = input_scale * output_scale
     pairwise(config, region, &[input, inv_denom], BaseOp::Mult)
@@ -6760,7 +6776,7 @@ pub(crate) fn percent<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>
 ///     Some(&[2, 2, 3, 2, 2, 0]),
 ///     &[2, 3],
 /// ).unwrap());
-/// let result = softmax::<Fp>(&dummy_config, &mut dummy_region, &[x], 128.0.into(), (128.0 * 128.0).into()).unwrap();
+/// let result = softmax::<Fp>(&dummy_config, &mut dummy_region, &[x], 128.0.into(), (128.0 * 128.0).into(), f64::EPSILON).unwrap();
 /// // doubles the scale of the input
 /// let expected = Tensor::<IntegerRep>::new(Some(&[350012, 350012, 352768, 350012, 350012, 344500]), &[2, 3]).unwrap();
 /// assert_eq!(result.int_evals().unwrap(), expected);
@@ -6771,6 +6787,7 @@ pub fn softmax<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
     values: &[ValTensor<F>; 1],
     input_scale: utils::F32,
     output_scale: utils::F32,
+    eps: f64,
 ) -> Result<ValTensor<F>, CircuitError> {
     // get the max then subtract it
     let max_val = max(config, region, values)?;
@@ -6787,7 +6804,14 @@ pub fn softmax<F: PrimeField + TensorType + PartialOrd + std::hash::Hash>(
         },
     )?;
 
-    percent(config, region, &[ex.clone()], input_scale, output_scale)
+    percent(
+        config,
+        region,
+        &[ex.clone()],
+        input_scale,
+        output_scale,
+        eps,
+    )
 }
 
 /// Checks that the percent error between the expected public output and the actual output value

src/graph/utilities.rs

@@ -858,6 +858,7 @@ pub fn new_op_from_onnx(
             SupportedOp::Hybrid(HybridOp::Recip {
                 input_scale: (scale_to_multiplier(in_scale) as f32).into(),
                 output_scale: (scale_to_multiplier(max_scale) as f32).into(),
+                eps: run_args.get_epsilon(),
             })
         }
 
@@ -903,6 +904,7 @@ pub fn new_op_from_onnx(
             SupportedOp::Hybrid(HybridOp::Rsqrt {
                 input_scale: (scale_to_multiplier(in_scale) as f32).into(),
                 output_scale: (scale_to_multiplier(max_scale) as f32).into(),
+                eps: run_args.get_epsilon(),
             })
         }
         "Exp" => SupportedOp::Nonlinear(LookupOp::Exp {
@@ -913,6 +915,7 @@ pub fn new_op_from_onnx(
             if run_args.bounded_log_lookup {
                 SupportedOp::Hybrid(HybridOp::Ln {
                     scale: scale_to_multiplier(input_scales[0]).into(),
+                    eps: run_args.get_epsilon(),
                 })
             } else {
                 SupportedOp::Nonlinear(LookupOp::Ln {
@@ -1131,6 +1134,7 @@ pub fn new_op_from_onnx(
                 input_scale: scale_to_multiplier(in_scale).into(),
                 output_scale: scale_to_multiplier(max_scale).into(),
                 axes: softmax_op.axes.to_vec(),
+                eps: run_args.get_epsilon(),
             })
         }
         "MaxPool" => {

src/lib.rs

@@ -97,11 +97,11 @@ impl From<String> for EZKLError {
 
 use std::str::FromStr;
 
-use circuit::{table::Range, CheckMode};
+use circuit::{CheckMode, table::Range};
 #[cfg(all(feature = "ezkl", not(target_arch = "wasm32")))]
 use clap::Args;
 use fieldutils::IntegerRep;
-use graph::{Visibility, MAX_PUBLIC_SRS};
+use graph::{MAX_PUBLIC_SRS, Visibility};
 use halo2_proofs::poly::{
     ipa::commitment::IPACommitmentScheme, kzg::commitment::KZGCommitmentScheme,
 };
@@ -350,6 +350,16 @@ pub struct RunArgs {
         arg(long, default_value = "false")
     )]
     pub ignore_range_check_inputs_outputs: bool,
+    /// Optional override for epsilon value
+    #[cfg_attr(all(feature = "ezkl", not(target_arch = "wasm32")), arg(long))]
+    pub epsilon: Option<f64>,
+}
+
+impl RunArgs {
+    /// Returns the epsilon value
+    pub fn get_epsilon(&self) -> f64 {
+        self.epsilon.unwrap_or(f64::EPSILON)
+    }
 }
 
 impl Default for RunArgs {
@@ -376,6 +386,7 @@ impl Default for RunArgs {
             decomp_base: 16384,
             decomp_legs: 2,
             ignore_range_check_inputs_outputs: false,
+            epsilon: None,
         }
     }
 }

src/tensor/ops.rs

@@ -160,7 +160,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, 0, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 0, 3, 4, 5, 6]), &[1, 3, 2]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 ///
 /// let result = trilu(&a, -1, true).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 2, 3, 4, 0, 6]), &[1, 3, 2]).unwrap();
@@ -168,7 +168,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, -1, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[0, 0, 3, 0, 5, 6]), &[1, 3, 2]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 ///
 /// let a = Tensor::<IntegerRep>::new(
 ///   Some(&[1, 2, 3, 4, 5, 6]),
@@ -188,7 +188,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, 0, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 0, 0, 4, 5, 0]), &[1, 2, 3]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 ///
 /// let result = trilu(&a, -1, true).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 2, 3, 4, 5, 6]), &[1, 2, 3]).unwrap();
@@ -196,7 +196,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, -1, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[0, 0, 0, 4, 0, 0]), &[1, 2, 3]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 ///
 /// let a = Tensor::<IntegerRep>::new(
 ///   Some(&[1, 2, 3, 4, 5, 6, 7, 8, 9]),
@@ -216,7 +216,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, 0, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 0, 0, 4, 5, 0, 7, 8, 9]), &[1, 3, 3]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 ///
 /// let result = trilu(&a, -1, true).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 2, 3, 4, 5, 6, 0, 8, 9]), &[1, 3, 3]).unwrap();
@@ -224,7 +224,7 @@ pub fn decompose(
 ///
 /// let result = trilu(&a, -1, false).unwrap();
 /// let expected = Tensor::<IntegerRep>::new(Some(&[0, 0, 0, 4, 0, 0, 7, 8, 0]), &[1, 3, 3]).unwrap();
-/// assert_eq!(result, expected);  
+/// assert_eq!(result, expected);
 /// ```
 pub fn trilu<T: TensorType + std::marker::Send + std::marker::Sync>(
     a: &Tensor<T>,
@@ -1859,14 +1859,14 @@ pub mod nonlinearities {
     ///     Some(&[4, 25, 8, 1, 1, 1]),
     ///     &[2, 3],
     /// ).unwrap();
-    /// let result = rsqrt(&x, 1.0);
+    /// let result = rsqrt(&x, 1.0, f64::EPSILON);
     /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 0, 0, 1, 1, 1]), &[2, 3]).unwrap();
     /// assert_eq!(result, expected);
     /// ```
-    pub fn rsqrt(a: &Tensor<IntegerRep>, scale_input: f64) -> Tensor<IntegerRep> {
+    pub fn rsqrt(a: &Tensor<IntegerRep>, scale_input: f64, eps: f64) -> Tensor<IntegerRep> {
         a.par_enum_map(|_, a_i| {
             let kix = (a_i as f64) / scale_input;
-            let fout = scale_input / (kix.sqrt() + f64::EPSILON);
+            let fout = scale_input / (kix.sqrt() + eps);
             let rounded = fout.round();
             Ok::<_, TensorError>(rounded as IntegerRep)
         })
@@ -2339,15 +2339,20 @@ pub mod nonlinearities {
     ///     &[2, 3],
     /// ).unwrap();
     /// let k = 2_f64;
-    /// let result = recip(&x, 1.0, k);
+    /// let result = recip(&x, 1.0, k, f64::EPSILON);
     /// let expected = Tensor::<IntegerRep>::new(Some(&[1, 2, 1, 0, 2, 2]), &[2, 3]).unwrap();
     /// assert_eq!(result, expected);
     /// ```
-    pub fn recip(a: &Tensor<IntegerRep>, input_scale: f64, out_scale: f64) -> Tensor<IntegerRep> {
+    pub fn recip(
+        a: &Tensor<IntegerRep>,
+        input_scale: f64,
+        out_scale: f64,
+        eps: f64,
+    ) -> Tensor<IntegerRep> {
         a.par_enum_map(|_, a_i| {
             let rescaled = (a_i as f64) / input_scale;
             let denom = if rescaled == 0_f64 {
-                (1_f64) / (rescaled + f64::EPSILON)
+                (1_f64) / (rescaled + eps)
             } else {
                 (1_f64) / (rescaled)
             };
@@ -2366,16 +2371,16 @@ pub mod nonlinearities {
     /// use ezkl::fieldutils::IntegerRep;
     /// use ezkl::tensor::ops::nonlinearities::zero_recip;
     /// let k = 2_f64;
-    /// let result = zero_recip(1.0);
+    /// let result = zero_recip(1.0, f64::EPSILON);
     /// let expected = Tensor::<IntegerRep>::new(Some(&[4503599627370496]), &[1]).unwrap();
     /// assert_eq!(result, expected);
     /// ```
-    pub fn zero_recip(out_scale: f64) -> Tensor<IntegerRep> {
+    pub fn zero_recip(out_scale: f64, eps: f64) -> Tensor<IntegerRep> {
         let a = Tensor::<IntegerRep>::new(Some(&[0]), &[1]).unwrap();
 
         a.par_enum_map(|_, a_i| {
             let rescaled = a_i as f64;
-            let denom = (1_f64) / (rescaled + f64::EPSILON);
+            let denom = (1_f64) / (rescaled + eps);
             let d_inv_x = out_scale * denom;
             Ok::<_, TensorError>(d_inv_x.round() as IntegerRep)
         })