libmath.c

#include "libmath.h"





/*               ********** Number Theory **********                */

int gcd(int num1, int num2) {
    if (num2 == 0)
        return num1;
    return gcd(num2, num1 % num2); //Euclidean algorithm
}

int lcm(int num1, int num2) {
    return (num1 / gcd(num1, num2)) * num2; 
}

bool is_prime(int num) {
    if (num < 2)
        return false;
    for (int i=2; i*i <= num; i++) { //number num will have divisors till sqrt of num.
        if (num % i == 0) 
            return false;
    }
    return true;
}

int prime_factorize(int num, int *primes) {
    int idx = 0;
    for (int i=2; i * i <= num; i++) {
        while (num % i == 0) {
            primes[idx++] = i;
            num /= i; //we found the prime i, now lets get it out of the multiples by dividing
        }
    }
    if (num > 1) 
        primes[idx++] = num;
    return idx;
}

int factorial( int n){
    return n<0 ? -1 : n == 0 ? 1 : n * factorial(n-1);
}

unsigned int divisor_sum(unsigned int n) { //unsigned int for a broader range of positive integers
    unsigned int sum = 1, pow = 2;
    unsigned int prime;
    for (; (n & 1) == 0; pow <<= 1, n >>= 1) { //another reason for unsigned int, better handling of bitwise operations, we directly work with binary representation.
        sum += pow;
    }
    for (prime = 3; prime * prime <= n; prime += 2) {
        unsigned int sum_in = 1;
        for (pow = prime; n % prime == 0; pow *= prime, n /= prime) {
            sum_in += pow;
        }
        sum *= sum_in;
    }
    if (n > 1) {
        sum *= n + 1;
    }
    return sum;
}

unsigned int phi(unsigned int n) {
    unsigned int count = 0;
    for (unsigned int i = 1; i <= n; i++) {
        if (gcd(i, n) == 1) {
            count++;
        }
    }
    return count;
}


/*               ********** Matrices **********                */

int **create_matrix(int row, int col) {
    int **mat = (int **)malloc(row * sizeof(int *));
    for (int i=0; i < row; i++) {
        mat[i] = (int *)malloc(col * sizeof(int)); 
    }
    return mat;
}

void free_matrix(int **mat, int row) {
    for (int i=0; i < row; i++) {
        free(mat[i]);
    }
    free(mat);
}

void matprint(int **mat, int row, int col) {
    for (int i=0; i < row; i++) {
        for (int j=0; j < col; j++) {
            printf("%d ", mat[i][j]);
        }
        printf("\n");
    }
}

void matrix_add(int **mat1, int **mat2, int **result, int row, int col){
    int *row_1, *row_2, *row_r; //a row of mat1,mat2 and reusult matrices(see pointer arithmetics)
    for(int i=0; i<row; i++){
        row_1 = mat1[i];
        row_2 = mat2[i];    //take current rows. So now ur problem is simplified to 3 1D arrays
        row_r = result[i];
        for(int j=0; j<col; j++){
            row_r[j] = row_1[j]+row_2[j]; 
        }
    }
}

void matmul(int **mat1, int **mat2, int **result, int row_1, int col_1, int col_2){
    for(int i=0; i < row_1; i++){
        for(int j=0; j < col_2; j++){
            result[i][j]=0;
        }
    }
    for(int i=0; i < row_1; i++){
        for(int j=0; j < col_2; j++){
            int sum = 0;
            for(int k=0; k < col_1; k+=4){ //4 iterations at time, faster computation, less iterations
                sum+= mat1[i][k] * mat2[i][k];
                if(k+1 < col_1) //!!!Since we are increasing the value of current index, we have to always make sure it is in the boundary.
                    sum+= mat1[i][k+1]*mat2[k+1][j];
                if(k+2 < col_1)
                    sum+= mat1[i][k+2] * mat2[k+2][j];
                if(k+3 < col_1)
                    sum+= mat1[i][k+3] * mat2[k+3][j];
            }
            result[i][j] = sum;
        }
    }
}

void transpose(void *mat_transposed, void *mat_original, int row, int col){
	double (*dest)[row] = mat_transposed;
    double (*src)[col] = mat_original;
	for (int i=0; i < row; i++)
		for (int j=0; j < col; j++)
			dest[j][i] = src[i][j];
}
double det_in(double **mat, int size, int perm){
	if (size == 1)
        return mat[0][0];

	double sum = 0, *m[--size];
	for (int i=0; i < size; i++)
		m[i] = mat[i + 1] + 1;
	for (int i=0, sign = 1; i <= size; i++) {
		sum += sign * (mat[i][0] * det_in(m, size, perm));
		if (i == size)
            break;
		m[i] = mat[i] + 1;
		if (!perm)
            sign = -sign;
	}
	return sum;
}

double det(double *mat, int size, int perm){ //wrapper function for the det_in
	double *m[size];
	for (int i=0; i < size; i++)
		m[i] = mat + (size * i);
	return det_in(m, size, perm);
}

void ring(unsigned int n){
    for(int i=0; i < n; i++){
        for(int j=0; j < n; j++){
            if(i ==0 || i == n-1 || j==0 || j == n - 1){
                printf("%d", 1);
            }
            else{
                printf("%d", 0);
            }
        }
        printf("\n");
    }
}
/*               ********** Statistics and Probability **********                */


void permutation(mpz_t result, int n, int k){
	mpz_set_ui(result, 1); //assign value
	k = n-k;
	while (n > k)
        mpz_mul_ui(result, result, n--); //multiply
}

void combination(mpz_t result, int n, int k)
{
	permutation(result, n, k);
	while (k)
        mpz_divexact_ui(result, result, k--); /*This one line replaces the 2 lines we normally would have for normal integers (relatively small ones that we can compute with int or long):
                                                result = result/k;
                                                k--;
                                                */
}

double mean(double *arr, int len){
	double sum = 0;
	for (int i=0; i < len; i++)
		sum += arr[i];
	return sum / len;
}

double median(int* arr, int len) {
    quicksort(arr, len);
    if (len % 2 == 1)
        return arr[len / 2]; 
    return (arr[(len / 2)-1] + arr[len / 2]) / 2.0; //when no element on middle, find average of two neighbours. !!! division should be with 2.0 not 2 for not getting a int type result !!!
}

float* benford(){
    static float digits[9];
    for(int i=0; i < 10; i++){
        digits[i-1] = log10f(1 + 1.0 / i);
    }
    return digits;
}

/*               ********** Helper functions for other functions **********                */


void quicksort(int *arr, int len) {
    if (len < 2)
        return;
    int piv = arr[len / 2]; //set pivot
    int i,j; //since I am using i in line 156 157 (out of loop), let's declare it here for keeping i after the loop (error was encountered when declaring inside the loop)
    for (i = 0, j = len - 1; ; i++, j--) { //start from end and start at the same time
        while (arr[i] < piv) i++;
        while (arr[j] > piv) j--;
        if (i >= j)
            break;
        int temp = arr[i];
        arr[i]= arr[j];
        arr[j]= temp;
    }
    quicksort(arr, i);
    quicksort(arr + i, len - i);
}

int b_search(int *arr, int num, int idx_start, int idx_end) {
    if (idx_end < idx_start) {
        return -1;
    }
    int k = idx_start + ((idx_end - idx_start) / 2);
    if (arr[k] == num) {
        return k;
    }
    else if (arr[k] < num) {
        return b_search(arr, num, k + 1, idx_end);
    }
    else {
        return b_search(arr, num, idx_start, k - 1);
    }
}
double sqrt(double n) {
    if (n < 0)
        return -1;  
    double x = n;  
    double precision = 1e-6;  
    while ((x*x - n) > precision || (n - x*x) > precision) {
        x = 0.5 * (x + n / x);
    }
    return x;
}


/*               ********** Vectors **********                */



vector to_cartesian(double r,double theta, double z){
	vector A;
	A.x = r*cos(theta);
	A.y = r*sin(theta);
    A.z = z;
	return A;
}

vector add_vectors(vector A,vector B){
	vector result;
	result.x = A.x + B.x;
	result.y = A.y + B.y;
    result.z = A.z + B.z;
	return result;
}

vector subtract_vectors(vector A,vector B){
	vector result;
	result.x = A.x - B.x;
	result.y = A.y - B.y;
	result.z = A.z - B.z;
	return result;
}

vector multiply_vector(vector A,double k){
	vector result;
	result.x = k * A.x;
	result.y = k * A.y;
    result.z = k * A.z;
	return result;
}

vector divide_vector(vector A,double k){
	vector result;
	result.x = A.x/k;
	result.y = A.y/k;
    result.z = A.z/k;
	return result;
}

void print_vector(vector A){
	printf("%lf %c %c %lf %c %c %lf %c",A.x,140,(A.y>=0)?'+':'-',(A.y>=0)?A.y:fabs(A.y),150,(A.z >= 0) ? '+' : '-', fabs(A.z), 'k');
}

double dot_product(vector A, vector B){
    return A.x*B.x + A.y*B.y + A.z*B.z;
}

double vector_module(vector A){
    return sqrt(A.x * A.x + A.y * A.y + A.z * A.z);
}

double vector_angle(vector A, vector B){
    double module_A = vector_module(A);
    double module_B = vector_module(B);
    if (module_A == 0 || module_B == 0) {
        printf("Error: Zero magnitude in one of the vectors.\n");
        return -1;
    }
    double product = dot_product(A,B);
    double cos_theta = product/(module_A * module_B);
    if (cos_theta > 1.0) //adjust floating point inaccuracies
        cos_theta = 1.0;
    if (cos_theta < -1.0)
        cos_theta = -1.0;
    return acos(cos_theta);
}

double distance(vector A, vector B){
    double x = B.x - A.x;
    double y = B.y - A.y;
    double z = B.z - A.z;
    return sqrt(x*x + y*y + z*z);
}

vector reflection(vector A, vector n){
    double A_n_dot = dot_product(A, n);
    vector A_r;
    A_r.x = A.x - 2 * A_n_dot * n.x;  //See this video for the proof of the formula if you are interested: https://www.youtube.com/watch?v=BCmFsYFln2k
    A_r.y = A.y - 2 * A_n_dot * n.y;
    A_r.z = A.z - 2 * A_n_dot * n.z;

    return A_r;
}

vector lerp(vector A, vector B, double dt){
    if(dt<0 || dt>1){
        printf("Incorrect interpolation factor! It should be between 0 and 1\n");
        return A;
    }
    vector lerp_vector;
    lerp_vector.x = A.x + dt*(B.x - A.x);
    lerp_vector.y = A.y + dt*(B.y - A.y);
    lerp_vector.z = A.z + dt*(B.z - A.z);

    return lerp_vector;
}