diff --git a/modules/dnn/src/layers/cpu_kernels/activation_kernels.simd.hpp b/modules/dnn/src/layers/cpu_kernels/activation_kernels.simd.hpp
index 161c1e0479..a28e57c441 100644
--- a/modules/dnn/src/layers/cpu_kernels/activation_kernels.simd.hpp
+++ b/modules/dnn/src/layers/cpu_kernels/activation_kernels.simd.hpp
@@ -232,10 +232,13 @@ static void activationGELUApprox(const void* input, void* output,
     v_float32 half = vx_setall_f32(0.5f), one = vx_setall_f32(1.f);
     v_float32 v_s2pi = vx_setall_f32(sqrt2_pi), v_coeff = vx_setall_f32(coeff);
     v_float32 two = vx_setall_f32(2.f);
+    // Clamp to [-9, 9] to prevent overflow in exp(2*inner); tanh saturates here anyway
+    v_float32 clamp_hi = vx_setall_f32(9.f), clamp_lo = vx_setall_f32(-9.f);
     for (; i + vlanes <= len; i += vlanes) {
         v_float32 x = vx_load(inp + i);
         // inner = sqrt(2/pi) * x + coeff * x^3 = x * (sqrt(2/pi) + coeff * x^2)
         v_float32 inner = v_mul(x, v_add(v_s2pi, v_mul(v_coeff, v_mul(x, x))));
+        inner = v_min(v_max(inner, clamp_lo), clamp_hi);
         // tanh via exp: (exp(2*inner)-1)/(exp(2*inner)+1)
         v_float32 e2 = v_exp(v_mul(two, inner));
         v_float32 t = v_div(v_sub(e2, one), v_add(e2, one));
diff --git a/modules/dnn/src/layers/cpu_kernels/softmax.cpp b/modules/dnn/src/layers/cpu_kernels/softmax.cpp
index aa5f683be6..e07be85e51 100644
--- a/modules/dnn/src/layers/cpu_kernels/softmax.cpp
+++ b/modules/dnn/src/layers/cpu_kernels/softmax.cpp
@@ -11,6 +11,7 @@
 
 #include "../../precomp.hpp"
 #include "softmax.hpp"
+#include "opencv2/core/fast_math.hpp"
 
 namespace cv { namespace dnn {
 
@@ -101,10 +102,13 @@ void softmax(Mat &dst, const Mat &src, int axis, int axisBias, int axisStep){
                 s += axisBuf[cnDim];
             }
 
-            s = 1.f / s;
-
             // copy back the result to src
             _cnDim = 0;
+            if (s == 0.f || cvIsInf(1.f / s)) {
+                for (; _cnDim < axisStep; _cnDim++)
+                    dstPtr[srcOffset + (_cnDim + axisBias) * cnStep] = 0.f;
+            } else {
+                s = 1.f / s;
 #if CV_ENABLE_UNROLLED && defined(_M_ARM64)
             for (; _cnDim + 3 < axisStep; _cnDim += 4) {
                 dstPtr[srcOffset + (_cnDim + 0 + axisBias) * cnStep] = axisBuf[_cnDim + 0] * s;
@@ -115,6 +119,7 @@ void softmax(Mat &dst, const Mat &src, int axis, int axisBias, int axisStep){
 #endif
             for (; _cnDim < axisStep; _cnDim++)
                 dstPtr[srcOffset + (_cnDim + axisBias) * cnStep] = axisBuf[_cnDim] * s;
+            }
         }
     }, nstripes);
 }
diff --git a/modules/dnn/test/test_layers.cpp b/modules/dnn/test/test_layers.cpp
index fba5f4ad78..5be5b4eee9 100644
--- a/modules/dnn/test/test_layers.cpp
+++ b/modules/dnn/test/test_layers.cpp
@@ -41,6 +41,7 @@
 
 #include "test_precomp.hpp"
 #include <opencv2/core/ocl.hpp>
+#include <opencv2/core/fast_math.hpp>
 #include "npy_blob.hpp"
 #include <opencv2/dnn/shape_utils.hpp>
 #include <opencv2/dnn/all_layers.hpp>
@@ -2959,4 +2960,59 @@ TEST(ConvolutionWinograd, Accuracy)
     normAssert(outLarge, refLarge, "Large input after small", 0.0, 0.0);
 }
 
+TEST(Layer_Test_GeluApprox, NoNaN_LargeInput)
+{
+    LayerParams lp;
+    lp.type = "GeluApproximation";
+    lp.name = "test_gelu_approx";
+    Ptr<Layer> layer = LayerFactory::createLayerInstance("GeluApproximation", lp);
+    ASSERT_TRUE(layer != nullptr);
+
+    float data[] = {-15.f, -10.f, -7.4f, -1.f, 0.f, 1.f, 5.f, 10.6f, 15.f, 20.f};
+    int dims[] = {1, 1, 10};
+    Mat inp(3, dims, CV_32F, data);
+    std::vector<Mat> inpVec = {inp};
+    std::vector<Mat> outVec;
+
+    runLayer(layer, inpVec, outVec);
+    ASSERT_EQ(outVec.size(), (size_t)1);
+
+    Mat& out = outVec[0];
+    for (int i = 0; i < 10; i++) {
+        float val = out.ptr<float>()[i];
+        EXPECT_FALSE(cvIsNaN(val)) << "NaN at index " << i << " (input=" << data[i] << ")";
+        EXPECT_FALSE(cvIsInf(val)) << "Inf at index " << i << " (input=" << data[i] << ")";
+    }
+
+    EXPECT_NEAR(out.ptr<float>()[9], 20.f, 0.01f);
+    EXPECT_NEAR(out.ptr<float>()[0], 0.f, 1e-6f);
+    EXPECT_NEAR(out.ptr<float>()[4], 0.f, 1e-6f);
+}
+
+TEST(Layer_Test_Softmax, NoNaN_AllNegInf)
+{
+    LayerParams lp;
+    lp.type = "Softmax";
+    lp.name = "test_softmax";
+    lp.set("axis", 1);
+    Ptr<Layer> layer = LayerFactory::createLayerInstance("Softmax", lp);
+    ASSERT_TRUE(layer != nullptr);
+
+    int dims[] = {1, 8};
+    Mat inp(2, dims, CV_32F, Scalar(-std::numeric_limits<float>::infinity()));
+    std::vector<Mat> inpVec = {inp};
+    std::vector<Mat> outVec;
+
+    runLayer(layer, inpVec, outVec);
+    ASSERT_EQ(outVec.size(), (size_t)1);
+
+    Mat& out = outVec[0];
+    for (int i = 0; i < 8; i++) {
+        float val = out.ptr<float>()[i];
+        EXPECT_FALSE(cvIsNaN(val)) << "NaN at index " << i;
+        EXPECT_FALSE(cvIsInf(val)) << "Inf at index " << i;
+        EXPECT_EQ(val, 0.f) << "Expected 0 at index " << i;
+    }
+}
+
 }} // namespace