diff --git a/src/blas.c b/src/blas.c
index c619256c..1f769280 100644
--- a/src/blas.c
+++ b/src/blas.c
@@ -515,3 +515,95 @@ void fix_nan_and_inf_cpu(float *input, size_t size)
             input[i] = 1.0f / i;  // pseudo random value
     }
 }
+
+
+float cosine_similarity(float *A, float *B, unsigned int feature_size)
+{
+    float mul = 0.0, d_a = 0.0, d_b = 0.0;
+
+    for (unsigned int i = 0; i < feature_size; ++i)
+    {
+        mul += A[i] * B[i];
+        d_a += A[i] * A[i];
+        d_b += B[i] * B[i];
+    }
+    float similarity;
+    float divider = sqrt(d_a) * sqrt(d_b);
+    if (divider > 0) similarity = mul / divider;
+    else similarity = 0;
+
+    return similarity;
+}
+
+// num_of_samples = 2 * loaded_images = mini_batch_size
+
+float P_constrastive(int i, int l, int num_of_samples, float **z, unsigned int feature_size, float temperature)
+{
+    if (i == l) {
+        printf(" Error: in P_constrastive must be i != l, while i = %d, l = %d \n", i, l);
+        getchar();
+    }
+
+    const float sim = cosine_similarity(z[i], z[l], feature_size);
+    const float numerator = expf(sim / temperature);
+
+    float denominator = 0;
+    int k;
+    for (k = 0; k < num_of_samples; ++k) {
+        if (k != i) {
+            const float sim_den = cosine_similarity(z[k], z[l], feature_size);
+            denominator += expf(sim_den / temperature);
+        }
+    }
+
+    float result = numerator / denominator;
+    return result;
+}
+
+// i - id of the current sample in mini_batch
+// labels[num_of_samples] - array with class_id for each sample in the current mini_batch
+// z[feature_size][num_of_samples] - array of arrays with contrastive features (output of conv-layer, f.e. 128 floats for each sample)
+// delta[feature_size] - array with deltas for backpropagation
+// temperature - scalar temperature param (temperature > 0), f.e. temperature = 0.07: Supervised Contrastive Learning
+void grad_contrastive_loss_positive(int i, int *labels, int num_of_samples, float **z, unsigned int feature_size, float temperature, float *delta)
+{
+    int j;
+    for (j = 0; j < num_of_samples; ++j) {
+        if (i != j && labels[i] == labels[j]) {
+            const double sim = cosine_similarity(z[i], z[j], feature_size);
+            const double P = P_constrastive(i, j, num_of_samples, z, feature_size, temperature);
+
+            int m;
+            for (m = 0; m < feature_size; ++m) {
+                delta[m] += (sim * z[i][m] - z[j][m]) * (1 - P);
+            }
+        }
+    }
+}
+
+// i - id of the current sample in mini_batch
+// labels[num_of_samples] - array with class_id for each sample in the current mini_batch
+// z[feature_size][num_of_samples] - array of arrays with contrastive features (output of conv-layer, f.e. 128 floats for each sample)
+// delta[feature_size] - array with deltas for backpropagation
+// temperature - scalar temperature param (temperature > 0), f.e. temperature = 0.07: Supervised Contrastive Learning
+void grad_contrastive_loss_negative(int i, int *labels, int num_of_samples, float **z, unsigned int feature_size, float temperature, float *delta)
+{
+    int j;
+    for (j = 0; j < num_of_samples; ++j) {
+        if (i != j && labels[i] == labels[j]) {
+
+            int k;
+            for (k = 0; k < num_of_samples; ++k) {
+                if (k != i && k != j) {
+                    const double sim = cosine_similarity(z[i], z[k], feature_size);
+                    const double P = P_constrastive(i, k, num_of_samples, z, feature_size, temperature);
+
+                    int m;
+                    for (m = 0; m < feature_size; ++m) {
+                        delta[m] += (z[k][m] - sim * z[i][m]) * P;
+                    }
+                }
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/blas.h b/src/blas.h
index 945d3919..473ac089 100644
--- a/src/blas.h
+++ b/src/blas.h
@@ -154,6 +154,11 @@ void rotate_weights_gpu(const float *src_weight_gpu, float *weight_deform_gpu, i
 void reduce_and_expand_array_gpu(const float *src_gpu, float *dst_gpu, int size, int groups);
 void expand_array_gpu(const float *src_gpu, float *dst_gpu, int size, int groups);
 
+float cosine_similarity(float *A, float *B, unsigned int feature_size);
+float P_constrastive(int i, int l, int num_of_samples, float **z, unsigned int feature_size, float temperature);
+void grad_contrastive_loss_positive(int i, int *labels, int num_of_samples, float **z, unsigned int feature_size, float temperature, float *delta);
+void grad_contrastive_loss_negative(int i, int *labels, int num_of_samples, float **z, unsigned int feature_size, float temperature, float *delta);
+
 #endif
 #ifdef __cplusplus
 }