Python и C/C++: ускоряем код в 10-100 раз с динамическими библиотеками

// factorial.c
#include <stdio.h>

// Экспортируемая функция для вычисления факториала
int factorial(int n) {
if (n <= 1) return 1;
return n * factorial(n – 1);
}

# test_factorial.py
import ctypes
import os

# Определяем путь к библиотеке в зависимости от ОС
if os.name == 'posix':
lib = ctypes.CDLL('./libfactorial.so')
else:
lib = ctypes.CDLL('./factorial.dll')

# Устанавливаем типы аргументов и возвращаемого значения
lib.factorial.argtypes = [ctypes.c_int]
lib.factorial.restype = ctypes.c_int

# Вызываем функцию
result = lib.factorial(10)
print(f"Factorial of 10 is: {result}")

// vector_ops.c
#include <stdlib.h>
#include <math.h>

typedef struct {
double x;
double y;
double z;
} Vector3D;

// Создание нового вектора
Vector3D* create_vector(double x, double y, double z) {
Vector3D* vec = (Vector3D*)malloc(sizeof(Vector3D));
vec->x = x;
vec->y = y;
vec->z = z;
return vec;
}

// Освобождение памяти
void free_vector(Vector3D* vec) {
free(vec);
}

// Вычисление длины вектора
double vector_length(Vector3D* vec) {
return sqrt(vec->x * vec->x + vec->y * vec->y + vec->z * vec->z);
}

// Скалярное произведение
double dot_product(Vector3D* vec1, Vector3D* vec2) {
return vec1->x * vec2->x + vec1->y * vec2->y + vec1->z * vec2->z;
}

import ctypes
import os
import math

# Загружаем библиотеку
if os.name == 'posix':
lib = ctypes.CDLL('./libvector_ops.so')
else:
lib = ctypes.CDLL('./vector_ops.dll')

# Определяем класс для Vector3D структуры
class Vector3D(ctypes.Structure):
_fields_ = [
("x", ctypes.c_double),
("y", ctypes.c_double),
("z", ctypes.c_double)
]

# Настраиваем типы функций
lib.create_vector.argtypes = [ctypes.c_double, ctypes.c_double, ctypes.c_double]
lib.create_vector.restype = ctypes.POINTER(Vector3D)

lib.free_vector.argtypes = [ctypes.POINTER(Vector3D)]
lib.free_vector.restype = None

lib.vector_length.argtypes = [ctypes.POINTER(Vector3D)]
lib.vector_length.restype = ctypes.c_double

lib.dot_product.argtypes = [ctypes.POINTER(Vector3D), ctypes.POINTER(Vector3D)]
lib.dot_product.restype = ctypes.c_double

# Использование библиотеки
v1 = lib.create_vector(1.0, 2.0, 3.0)
v2 = lib.create_vector(4.0, 5.0, 6.0)

length1 = lib.vector_length(v1)
dot = lib.dot_product(v1, v2)

print(f"Vector length: {length1}")
print(f"Dot product: {dot}")

# Не забываем освободить память
lib.free_vector(v1)
lib.free_vector(v2)

// array_ops.c
#include <stdlib.h>

// Функция для поэлементного умножения массивов
void multiply_arrays(double* arr1, double* arr2, double* result, int size) {
for (int i = 0; i < size; i++) {
result[i] = arr1[i] * arr2[i];
}
}

// Функция для поиска максимального элемента
double find_max(double* arr, int size) {
double max = arr[0];
for (int i = 1; i < size; i++) {
if (arr[i] > max) {
max = arr[i];
}
}
return max;
}

import ctypes
import os
import numpy as np

# Загружаем библиотеку
if os.name == 'posix':
lib = ctypes.CDLL('./libarray_ops.so')
else:
lib = ctypes.CDLL('./array_ops.dll')

# Настраиваем типы функций
lib.multiply_arrays.argtypes = [
ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double),
ctypes.c_int
]
lib.multiply_arrays.restype = None

lib.find_max.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.c_int]
lib.find_max.restype = ctypes.c_double

# Создаем тестовые данные с использованием NumPy
arr1 = np.array([1\.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float64)
arr2 = np.array([5\.0, 4.0, 3.0, 2.0, 1.0], dtype=np.float64)
result = np.zeros_like(arr1)

# Получаем указатели на данные массивов NumPy
arr1_ptr = arr1.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
arr2_ptr = arr2.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
result_ptr = result.ctypes.data_as(ctypes.POINTER(ctypes.c_double))

# Вызываем функцию умножения массивов
lib.multiply_arrays(arr1_ptr, arr2_ptr, result_ptr, len(arr1))

print("Результат умножения массивов:", result)

# Находим максимальный элемент
max_val = lib.find_max(result_ptr, len(result))
print("Максимальный элемент:", max_val)

pip install cffi

import numpy as np
from cffi import FFI

ffi = FFI()

# Определяем интерфейс C-функций
ffi.cdef("""
void multiply_arrays(double* arr1, double* arr2, double* result, int size);
double find_max(double* arr, int size);
""")

# Загружаем библиотеку
if os.name == 'posix':
lib = ffi.dlopen('./libarray_ops.so')
else:
lib = ffi.dlopen('./array_ops.dll')

# Создаем тестовые данные
arr1 = np.array([1\.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float64)
arr2 = np.array([5\.0, 4.0, 3.0, 2.0, 1.0], dtype=np.float64)
result = np.zeros_like(arr1)

# Преобразуем NumPy массивы в формат CFFI
arr1_ptr = ffi.cast("double*", arr1.ctypes.data)
arr2_ptr = ffi.cast("double*", arr2.ctypes.data)
result_ptr = ffi.cast("double*", result.ctypes.data)

# Вызываем функции библиотеки
lib.multiply_arrays(arr1_ptr, arr2_ptr, result_ptr, len(arr1))
print("Результат умножения массивов:", result)

max_val = lib.find_max(result_ptr, len(result))
print("Максимальный элемент:", max_val)

# build_array_ops.py
from cffi import FFI

ffi = FFI()

# Определяем C-код функций и заголовков
ffi.set_source("_array_ops", """
// Функции для работы с массивами
void multiply_arrays(double* arr1, double* arr2, double* result, int size) {
for (int i = 0; i < size; i++) {
result[i] = arr1[i] * arr2[i];
}
}

double find_max(double* arr, int size) {
double max = arr[0];
for (int i = 1; i < size; i++) {
if (arr[i] > max) {
max = arr[i];
}
}
return max;
}
""")

ffi.cdef("""
void multiply_arrays(double* arr1, double* arr2, double* result, int size);
double find_max(double* arr, int size);
""")

if __name__ == "__main__":
ffi.compile(verbose=True)

import numpy as np
from _array_ops import ffi, lib

# Создаем тестовые данные
arr1 = np.array([1\.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float64)
arr2 = np.array([5\.0, 4.0, 3.0, 2.0, 1.0], dtype=np.float64)
result = np.zeros_like(arr1)

# Преобразуем NumPy массивы в формат CFFI
arr1_ptr = ffi.cast("double*", arr1.ctypes.data)
arr2_ptr = ffi.cast("double*", arr2.ctypes.data)
result_ptr = ffi.cast("double*", result.ctypes.data)

# Вызываем функции библиотеки
lib.multiply_arrays(arr1_ptr, arr2_ptr, result_ptr, len(arr1))
print("Результат умножения массивов:", result)

max_val = lib.find_max(result_ptr, len(result))
print("Максимальный элемент:", max_val)

pip install cython

# pure_python.py
def fibonacci(n):
if n <= 1:
return n
return fibonacci(n-1) + fibonacci(n-2)

# Проверка производительности
import time
start = time.time()
result = fibonacci(30)
end = time.time()
print(f"Result: {result}, Time: {end – start:.2f} seconds")

# fibonacci.pyx
def fibonacci_py(n):
if n <= 1:
return n
return fibonacci_py(n-1) + fibonacci_py(n-2)

# Типизированная версия
cpdef int fibonacci_cy(int n):
if n <= 1:
return n
return fibonacci_cy(n-1) + fibonacci_cy(n-2)

# Чистая C-версия
cdef int fibonacci_c(int n) nogil:
if n <= 1:
return n
return fibonacci_c(n-1) + fibonacci_c(n-2)

# Обертка для вызова C-версии из Python
def fibonacci_c_wrapper(n):
return fibonacci_c(n)

# setup.py
from setuptools import setup, Extension
from Cython.Build import cythonize

extensions = [
Extension(
"fibonacci",
["fibonacci.pyx"],
)
]

setup(
name="Fibonacci",
ext_modules=cythonize(extensions, compiler_directives={'language_level': "3"}),
)

python setup.py build_ext --inplace

# benchmark.py
import time
from fibonacci import fibonacci_py, fibonacci_cy, fibonacci_c_wrapper

def benchmark(func, n):
start = time.time()
result = func(n)
end = time.time()
print(f"{func.__name__}: Result = {result}, Time = {end – start:.6f} seconds")

# Тестируем все три версии
n = 30
benchmark(fibonacci_py, n)
benchmark(fibonacci_cy, n)
benchmark(fibonacci_c_wrapper, n)

# fast_array_ops.pyx
import numpy as np
cimport numpy as np

# Объявляем типы для NumPy массивов
ctypedef np.float64_t DTYPE_t

def py_multiply_arrays(np.ndarray[DTYPE_t, ndim=1] arr1, 
np.ndarray[DTYPE_t, ndim=1] arr2):
cdef int i, n = arr1.shape[0]
cdef np.ndarray[DTYPE_t, ndim=1] result = np.zeros(n, dtype=np.float64)

for i in range(n):
result[i] = arr1[i] * arr2[i]

return result

# Оптимизированная версия с использованием указателей
def fast_multiply_arrays(np.ndarray[DTYPE_t, ndim=1] arr1, 
np.ndarray[DTYPE_t, ndim=1] arr2):
cdef int i, n = arr1.shape[0]
cdef np.ndarray[DTYPE_t, ndim=1] result = np.zeros(n, dtype=np.float64)

# Получаем прямые указатели на данные
cdef DTYPE_t* arr1_ptr = &arr1[0]
cdef DTYPE_t* arr2_ptr = &arr2[0]
cdef DTYPE_t* result_ptr = &result[0]

# Используем указатели для более быстрого доступа
for i in range(n):
result_ptr[i] = arr1_ptr[i] * arr2_ptr[i]

return result

# math_wrapper.pyx
cdef extern from "math.h":
double sin(double x)
double cos(double x)
double sqrt(double x)

def fast_sin(double x):
return sin(x)

def fast_cos(double x):
return cos(x)

def hypotenuse(double a, double b):
return sqrt(a*a + b*b)

# parallel_sum.pyx
import numpy as np
cimport numpy as np
from cython.parallel import prange

ctypedef np.float64_t DTYPE_t

def parallel_sum(np.ndarray[DTYPE_t, ndim=1] arr):
cdef int i, n = arr.shape[0]
cdef double total = 0.0

# Используем prange для параллельной обработки
for i in prange(n, nogil=True):
total += arr[i]

return total

# setup.py с поддержкой OpenMP
from setuptools import setup, Extension
from Cython.Build import cythonize
import numpy as np

extensions = [
Extension(
"parallel_sum",
["parallel_sum.pyx"],
extra_compile_args=['-fopenmp'],
extra_link_args=['-fopenmp'],
)
]

setup(
ext_modules=cythonize(extensions),
include_dirs=[np.get_include()],
)

// kmeans.c
#include <stdlib.h>
#include <math.h>
#include <float.h>

// Вычисление расстояния между двумя точками
double distance(double* point1, double* point2, int dims) {
double sum = 0.0;
for (int i = 0; i < dims; i++) {
double diff = point1[i] – point2[i];
sum += diff * diff;
}
return sqrt(sum);
}

// Алгоритм k-means
void kmeans(double* data, int n_points, int dims, int k, 
double* centroids, int* labels, int max_iter) {
// Инициализация центроидов (берем первые k точек из данных)
for (int i = 0; i < k; i++) {
for (int j = 0; j < dims; j++) {
centroids[i * dims + j] = data[i * dims + j];
}
}

for (int iter = 0; iter < max_iter; iter++) {
// Назначение меток точкам
for (int i = 0; i < n_points; i++) {
double min_dist = DBL_MAX;
int min_cluster = 0;

for (int j = 0; j < k; j++) {
double dist = distance(&data[i * dims], &centroids[j * dims], dims);
if (dist < min_dist) {
min_dist = dist;
min_cluster = j;
}
}

labels[i] = min_cluster;
}

// Пересчет центроидов
double* new_centroids = (double*)calloc(k * dims, sizeof(double));
int* counts = (int*)calloc(k, sizeof(int));

for (int i = 0; i < n_points; i++) {
int cluster = labels[i];
counts[cluster]++;

for (int j = 0; j < dims; j++) {
new_centroids[cluster * dims + j] += data[i * dims + j];
}
}

for (int i = 0; i < k; i++) {
if (counts[i] > 0) {
for (int j = 0; j < dims; j++) {
centroids[i * dims + j] = new_centroids[i * dims + j] / counts[i];
}
}
}

free(new_centroids);
free(counts);
}
}

import numpy as np
from cffi import FFI

ffi = FFI()
ffi.cdef("""
void kmeans(double* data, int n_points, int dims, int k, 
double* centroids, int* labels, int max_iter);
""")

# Загружаем библиотеку
lib = ffi.dlopen('./libkmeans.so') # или kmeans.dll в Windows

# Генерируем тестовые данные
np.random.seed(42)
data = np.random.rand(1000, 2).astype(np.float64) # 1000 точек в 2D пространстве

# Подготавливаем параметры
n_points, dims = data.shape
k = 5 # Количество кластеров
centroids = np.zeros((k, dims), dtype=np.float64)
labels = np.zeros(n_points, dtype=np.int32)
max_iter = 100

# Вызываем функцию C
lib.kmeans(
ffi.cast("double*", data.ctypes.data),
n_points,
dims,
k,
ffi.cast("double*", centroids.ctypes.data),
ffi.cast("int*", labels.ctypes.data),
max_iter
)

print("Центроиды кластеров:")
print(centroids)
print("\nРаспределение точек по кластерам:")
for i in range(k):
print(f"Кластер {i}: {np.sum(labels == i)} точек")

# Визуализация результатов
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
for i in range(k):
mask = labels == i
plt.scatter(data[mask, 0], data[mask, 1], label=f'Кластер {i}')
plt.scatter(centroids[:, 0], centroids[:, 1], c='black', marker='x', s=100, label='Центроиды')
plt.legend()
plt.title('Результаты k-means кластеризации')
plt.show()

// gaussian_blur.c
#include <stdlib.h>
#include <math.h>

// Создание ядра фильтра Гаусса
void create_gaussian_kernel(double* kernel, int radius, double sigma) {
int size = 2 * radius + 1;
double sum = 0.0;

for (int y = -radius; y <= radius; y++) {
for (int x = -radius; x <= radius; x++) {
double value = exp(-(x*x + y*y) / (2 * sigma * sigma));
kernel[(y + radius) * size + (x + radius)] = value;
sum += value;
}
}

// Нормализация ядра
for (int i = 0; i < size * size; i++) {
kernel[i] /= sum;
}
}

// Применение фильтра размытия Гаусса
void apply_gaussian_blur(unsigned char* input, unsigned char* output, 
int width, int height, int channels,
double* kernel, int radius) {
int size = 2 * radius + 1;

for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
for (int c = 0; c < channels; c++) {
double sum = 0.0;

for (int ky = -radius; ky <= radius; ky++) {
for (int kx = -radius; kx <= radius; kx++) {
int img_x = x + kx;
int img_y = y + ky;

// Проверка границ
if (img_x >= 0 && img_x < width && img_y >= 0 && img_y < height) {
int img_idx = (img_y * width + img_x) * channels + c;
int kernel_idx = (ky + radius) * size + (kx + radius);

sum += input[img_idx] * kernel[kernel_idx];
}
}
}

output[(y * width + x) * channels + c] = (unsigned char)sum;
}
}
}
}

import ctypes
import numpy as np
from PIL import Image

# Загружаем библиотеку
lib = ctypes.CDLL('./libgaussian_blur.so') # или gaussian_blur.dll в Windows

# Настраиваем типы функций
lib.create_gaussian_kernel.argtypes = [
ctypes.POINTER(ctypes.c_double),
ctypes.c_int,
ctypes.c_double
]
lib.create_gaussian_kernel.restype = None

lib.apply_gaussian_blur.argtypes = [
ctypes.POINTER(ctypes.c_ubyte),
ctypes.POINTER(ctypes.c_ubyte),
ctypes.c_int,
ctypes.c_int,
ctypes.c_int,
ctypes.POINTER(ctypes.c_double),
ctypes.c_int
]
lib.apply_gaussian_blur.restype = None

# Загружаем изображение
img = Image.open('input.jpg')
img_array = np.array(img, dtype=np.uint8)
height, width, channels = img_array.shape

# Создаем выходной массив
output_array = np.zeros_like(img_array)

# Параметры фильтра
radius = 2
sigma = 1.0
kernel_size = 2 * radius + 1
kernel = np.zeros((kernel_size, kernel_size), dtype=np.float64)

# Создаем ядро фильтра
kernel_ptr = kernel.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
lib.create_gaussian_kernel(kernel_ptr, radius, sigma)

# Применяем фильтр
input_ptr = img_array.ctypes.data_as(ctypes.POINTER(ctypes.c_ubyte))
output_ptr = output_array.ctypes.data_as(ctypes.POINTER(ctypes.c_ubyte))

lib.apply_gaussian_blur(
input_ptr, output_ptr,
width, height, channels,
kernel_ptr, radius
)

# Создаем новое изображение из результата
blurred_img = Image.fromarray(output_array)
blurred_img.save('output_blurred.jpg')
blurred_img.show()

// physics_2d.cpp
#include <cmath>
#include <vector>
#include <algorithm>

extern "C" {

struct Particle {
double x, y; // Позиция
double vx, vy; // Скорость
double radius; // Радиус
double mass; // Масса
};

// Обновление позиции частиц
void update_positions(Particle* particles, int n_particles, double dt) {
for (int i = 0; i < n_particles; i++) {
particles[i].x += particles[i].vx * dt;
particles[i].y += particles[i].vy * dt;
}
}

// Обработка столкновений со стенами
void handle_wall_collisions(Particle* particles, int n_particles, 
double width, double height, double restitution) {
for (int i = 0; i < n_particles; i++) {
Particle& p = particles[i];

// Столкновение с левой и правой стенами
if (p.x – p.radius < 0) {
p.x = p.radius;
p.vx = -p.vx * restitution;
} else if (p.x + p.radius > width) {
p.x = width – p.radius;
p.vx = -p.vx * restitution;
}

// Столкновение с верхней и нижней стенами
if (p.y – p.radius < 0) {
p.y = p.radius;
p.vy = -p.vy * restitution;
} else if (p.y + p.radius > height) {
p.y = height – p.radius;
p.vy = -p.vy * restitution;
}
}
}

// Обработка столкновений между частицами
void handle_particle_collisions(Particle* particles, int n_particles, double restitution) {
for (int i = 0; i < n_particles; i++) {
for (int j = i + 1; j < n_particles; j++) {
Particle& p1 = particles[i];
Particle& p2 = particles[j];

double dx = p2.x – p1.x;
double dy = p2.y – p1.y;
double distance = sqrt(dx * dx + dy * dy);

// Проверка на столкновение
if (distance < p1.radius + p2.radius) {
// Нормализация вектора столкновения
double nx = dx / distance;
double ny = dy / distance;

// Расстояние перекрытия
double overlap = p1.radius + p2.radius – distance;

// Корректировка позиций для предотвращения перекрытия
p1.x -= nx * overlap * 0.5;
p1.y -= ny * overlap * 0.5;
p2.x += nx * overlap * 0.5;
p2.y += ny * overlap * 0.5;

// Вычисление относительной скорости
double dvx = p2.vx – p1.vx;
double dvy = p2.vy – p1.vy;

// Проекция скорости на нормаль столкновения
double vn = dvx * nx + dvy * ny;

// Если частицы удаляются друг от друга, столкновения не происходит
if (vn > 0) continue;

// Вычисление импульса
double impulse = -(1.0 + restitution) * vn;
impulse /= (1.0 / p1.mass) + (1.0 / p2.mass);

// Применение импульса
p1.vx -= impulse * nx / p1.mass;
p1.vy -= impulse * ny / p1.mass;
p2.vx += impulse * nx / p2.mass;
p2.vy += impulse * ny / p2.mass;
}
}
}
}

// Основная функция симуляции
void simulate_step(Particle* particles, int n_particles, 
double width, double height, double dt, double restitution) {
update_positions(particles, n_particles, dt);
handle_wall_collisions(particles, n_particles, width, height, restitution);
handle_particle_collisions(particles, n_particles, restitution);
}

} // extern "C"

import ctypes
import numpy as np
import pygame
import random

# Определяем структуру Particle для ctypes
class Particle(ctypes.Structure):
_fields_ = [
("x", ctypes.c_double),
("y", ctypes.c_double),
("vx", ctypes.c_double),
("vy", ctypes.c_double),
("radius", ctypes.c_double),
("mass", ctypes.c_double)
]

# Загружаем библиотеку
lib = ctypes.CDLL('./libphysics_2d.so') # или physics_2d.dll в Windows

# Настраиваем типы функций
lib.simulate_step.argtypes = [
ctypes.POINTER(Particle),
ctypes.c_int,
ctypes.c_double,
ctypes.c_double,
ctypes.c_double,
ctypes.c_double
]
lib.simulate_step.restype = None

# Инициализация pygame
pygame.init()
width, height = 800, 600
screen = pygame.display.set_mode((width, height))
pygame.display.set_caption("2D Physics Simulation")
clock = pygame.time.Clock()

# Создаем частицы
n_particles = 50
particles = (Particle * n_particles)()

for i in range(n_particles):
particles[i].radius = random.uniform(10.0, 30.0)
particles[i].ma