pxfm/exponents/logisticf.rs
1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1. Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2. Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3. Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::exponents::core_expf;
30use crate::polyeval::f_polyeval6;
31
32/// Logistic function
33///
34/// Computes 1/(1+exp(-x))
35pub fn f_logisticf(x: f32) -> f32 {
36 let x_u = x.to_bits();
37 let x_abs = x_u & 0x7fff_ffffu32;
38 if x_abs >= 0x42b2_0000u32 || x_abs <= 0x3d00_0000u32 {
39 // |x| >= 89 or |x| <= 2^-5
40 let exp = ((x_u >> 23) & 0xFF) as i32;
41 // |x| < 2^-25
42 if exp <= 101i32 {
43 return (1. / (2. + -x as f64)) as f32;
44 }
45
46 if x_abs <= 0x3400_0000u32 {
47 // |x| < f32::EPSILON
48 // logistic(x) ~ 0.5 + x / 4 + O(x^3)
49 #[cfg(any(
50 all(
51 any(target_arch = "x86", target_arch = "x86_64"),
52 target_feature = "fma"
53 ),
54 all(target_arch = "aarch64", target_feature = "neon")
55 ))]
56 {
57 use crate::common::f_fmlaf;
58 return f_fmlaf(0.25, x, 0.5);
59 }
60 #[cfg(not(any(
61 all(
62 any(target_arch = "x86", target_arch = "x86_64"),
63 target_feature = "fma"
64 ),
65 all(target_arch = "aarch64", target_feature = "neon")
66 )))]
67 {
68 use crate::common::f_fmla;
69 return f_fmla(0.25, x as f64, 0.5) as f32;
70 }
71 }
72
73 // When x < log(2^-150) or nan
74 if x_u >= 0xc2cf_f1b5u32 {
75 // logistic(-Inf) = 0
76 if x.is_infinite() {
77 return 0.0;
78 }
79 // logistic(nan) = nan
80 if x.is_nan() {
81 return x;
82 }
83 return if x.is_sign_positive() {
84 f32::INFINITY
85 } else {
86 0.
87 };
88 }
89 if x.is_sign_positive() && x_u >= 0x42d0_0000u32 {
90 // x >= 104
91 if x.is_nan() {
92 return f32::NAN;
93 }
94 return 1.;
95 }
96 if x.is_sign_negative() && x_u >= 0xc2cb_fe00u32 {
97 // x <= -101.99609
98 return 0.;
99 }
100
101 if x_abs <= 0x3d00_0000u32 {
102 // |x| <= 2^-5
103
104 // Polynomial for logistics function generated by Sollya:
105 // d = [-2^-5, 2^-5];
106 // f_logistic = 1/(1+exp(-x));
107 // Q = fpminimax(f_logistic, 5, [|D...|], d);
108
109 let dx = x as f64;
110 let p = f_polyeval6(
111 dx,
112 f64::from_bits(0x3fe0000000000000),
113 f64::from_bits(0x3fcffffffffffcf5),
114 f64::from_bits(0x3cfbce9c01a4987c),
115 f64::from_bits(0xbf955555524c34a4),
116 f64::from_bits(0xbda6322b5fac64d4),
117 f64::from_bits(0x3f61104f3a0dedc3),
118 );
119 return p as f32;
120 }
121 }
122
123 let v_exp = core_expf(-x);
124 let logistic = 1. / (1. + v_exp);
125 logistic as f32
126}
127
128#[cfg(test)]
129mod tests {
130 use super::*;
131
132 #[test]
133 fn test_logisticf() {
134 assert_eq!(f_logisticf(-89.), 2.227364e-39);
135 assert_eq!(f_logisticf(-104.), 0.);
136 assert_eq!(f_logisticf(-103.), 0.);
137 assert_eq!(f_logisticf(-88.9), 2.461613e-39);
138 assert_eq!(f_logisticf(-88.), 6.054601e-39);
139 assert_eq!(f_logisticf(-80.), 1.8048513e-35);
140 assert_eq!(f_logisticf(-60.), 8.756511e-27);
141 assert_eq!(f_logisticf(-40.), 4.248354e-18);
142 assert_eq!(f_logisticf(-20.), 2.0611537e-9);
143 assert_eq!(f_logisticf(-1.591388e29), 0.);
144 assert_eq!(f_logisticf(-3.), 0.047425874);
145 assert_eq!(f_logisticf(3.), 0.95257413);
146 assert_eq!(f_logisticf(20.), 1.);
147 assert_eq!(f_logisticf(55.), 1.);
148 assert_eq!(f_logisticf(-104.), 0.);
149 assert_eq!(f_logisticf(-90.), 8.19401e-40);
150 assert_eq!(f_logisticf(0.00000000000524323), 0.5);
151 assert_eq!(f_logisticf(0.00000000524323), 0.5);
152 assert_eq!(f_logisticf(0.02), 0.5049998);
153 assert_eq!(f_logisticf(-0.02), 0.49500015);
154 assert_eq!(f_logisticf(90.), 1.);
155 assert_eq!(f_logisticf(105.), 1.);
156 assert_eq!(f_logisticf(150.), 1.);
157 assert_eq!(f_logisticf(f32::INFINITY), 1.0);
158 assert_eq!(f_logisticf(f32::NEG_INFINITY), 0.);
159 assert!(f_logisticf(f32::NAN).is_nan());
160 }
161}