mirror of
https://hub.njuu.cf/TheAlgorithms/C-Plus-Plus.git
synced 2023-10-11 13:05:55 +08:00
Merge branch 'master' of https://github.com/ayaankhan98/C-Plus-Plus
This commit is contained in:
Ayaan Khan
commit
9217b4d8a8
@ -25,7 +25,7 @@ if(DOXYGEN_FOUND)
|
||||
set(DOXYGEN_GENERATE_MAN NO)
|
||||
set(DOXYGEN_USE_MATHJAX YES)
|
||||
set(DOXYGEN_GENERATE_HTML YES)
|
||||
set(DOXYGEN_HTML_TIMESTAMP YES)
|
||||
# set(DOXYGEN_HTML_TIMESTAMP YES)
|
||||
set(DOXYGEN_EXTRACT_STATIC YES)
|
||||
set(DOXYGEN_INLINE_SOURCES YES)
|
||||
set(DOXYGEN_CREATE_SUBDIRS YES)
|
||||
@ -68,6 +68,7 @@ endif()
|
||||
add_subdirectory(math)
|
||||
add_subdirectory(others)
|
||||
add_subdirectory(search)
|
||||
add_subdirectory(ciphers)
|
||||
add_subdirectory(strings)
|
||||
add_subdirectory(sorting)
|
||||
add_subdirectory(geometry)
|
||||
|
||||
@ -10,6 +10,9 @@
|
||||
* [Rat Maze](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/backtracking/rat_maze.cpp)
|
||||
* [Sudoku Solve](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/backtracking/sudoku_solve.cpp)
|
||||
|
||||
## Ciphers
|
||||
* [Hill Cipher](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/ciphers/hill_cipher.cpp)
|
||||
|
||||
## Data Structures
|
||||
* [Avltree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/data_structures/avltree.cpp)
|
||||
* [Binary Search Tree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/data_structures/binary_search_tree.cpp)
|
||||
@ -98,10 +101,14 @@
|
||||
* [Adaline Learning](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/adaline_learning.cpp)
|
||||
* [Kohonen Som Topology](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/kohonen_som_topology.cpp)
|
||||
* [Kohonen Som Trace](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/kohonen_som_trace.cpp)
|
||||
* [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/ordinary_least_squares_regressor.cpp)
|
||||
|
||||
## Math
|
||||
* [Armstrong Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/armstrong_number.cpp)
|
||||
* [Binary Exponent](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/binary_exponent.cpp)
|
||||
* [Check Amicable Pair](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/check_amicable_pair.cpp)
|
||||
* [Check Prime](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/check_prime.cpp)
|
||||
* [Complex Numbers](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/complex_numbers.cpp)
|
||||
* [Double Factorial](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/double_factorial.cpp)
|
||||
* [Eulers Totient Function](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/eulers_totient_function.cpp)
|
||||
* [Extended Euclid Algorithm](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/extended_euclid_algorithm.cpp)
|
||||
@ -137,11 +144,11 @@
|
||||
* [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)
|
||||
* [Golden Search Extrema](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/golden_search_extrema.cpp)
|
||||
* [Lu Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decompose.cpp)
|
||||
* [Lu Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decomposition.h)
|
||||
* [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/newton_raphson_method.cpp)
|
||||
* [Ode Forward Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_forward_euler.cpp)
|
||||
* [Ode Midpoint Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_midpoint_euler.cpp)
|
||||
* [Ode Semi Implicit Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_semi_implicit_euler.cpp)
|
||||
* [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ordinary_least_squares_regressor.cpp)
|
||||
* [Qr Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decompose.h)
|
||||
* [Qr Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decomposition.cpp)
|
||||
* [Qr Eigen Values](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_eigen_values.cpp)
|
||||
|
||||
38
README.md
38
README.md
@ -1,18 +1,36 @@
|
||||
# The Algorithms - C++ # {#mainpage}
|
||||
[](https://gitpod.io/#https://github.com/TheAlgorithms/C-Plus-Plus)
|
||||
<!-- the suffix in the above line is required for doxygen to consider this as the index page of the generated documentation site -->
|
||||
|
||||
[](https://gitpod.io/#https://github.com/TheAlgorithms/C-Plus-Plus)
|
||||
[](https://lgtm.com/projects/g/TheAlgorithms/C-Plus-Plus/context:cpp)
|
||||
[](https://gitter.im/TheAlgorithms)
|
||||
<a href="https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md"></a>
|
||||
[](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md)
|
||||
|
||||
|
||||
<a href="https://TheAlgorithms.github.io/C-Plus-Plus"></a>
|
||||
<a href="https://github.com/TheAlgorithms/C-Plus-Plus/actions?query=workflow%3A%22Awesome+CI+Workflow%22"></a>
|
||||
[](https://TheAlgorithms.github.io/C-Plus-Plus)
|
||||
[](https://github.com/TheAlgorithms/C-Plus-Plus/actions?query=workflow%3A%22Awesome+CI+Workflow%22)
|
||||
|
||||
[Online Documentation](https://TheAlgorithms.github.io/C-Plus-Plus).
|
||||
## Overview
|
||||
|
||||
The repository is a collection of open-source implementation of a variety of algorithms implemented in C++ and licensed under [MIT License](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/LICENSE). The algorithms span a variety of topics from computer science, mathematics and statistics, data science, machine learning, engineering, etc.. The implementations and the associated documentation are meant to provide a learning resource for educators and students. Hence, one may find more than one implementation for the same objective but using a different algorithm strategies and optimizations.
|
||||
|
||||
## Features
|
||||
|
||||
* The repository provides implementations of various algorithms in one of the most fundamental general purpose languages - [C++](https://en.wikipedia.org/wiki/C%2B%2B).
|
||||
* Well documented source code with detailed explanations provide a valuable resource for educators and students alike.
|
||||
* Each source code is atomic using [STL classes](https://en.wikipedia.org/wiki/Standard_Template_Library) and _no external libraries_ are required for their compilation and execution. Thus the fundamentals of the algorithms can be studied in much depth.
|
||||
* Source codes are [compiled and tested](https://github.com/TheAlgorithms/C-Plus-Plus/actions?query=workflow%3A%22Awesome+CI+Workflow%22) for every commit on the latest versions of three major operating systems viz., Windows, MacOS and Ubuntu (Linux) using MSVC 16 2019, AppleClang 11.0 and GNU 7.5.0 respectively.
|
||||
* Strict adherence to [C++11](https://en.wikipedia.org/wiki/C%2B%2B11) standard ensures portability of code to embedded systems as well like ESP32, ARM Cortex, etc. with little to no changes.
|
||||
* Self-checks within programs ensure correct implementations with confidence.
|
||||
* Modular implementations and OpenSource licensing enable the functions to be utilized conveniently in other applications.
|
||||
|
||||
## Documentation
|
||||
|
||||
[Online Documentation](https://TheAlgorithms.github.io/C-Plus-Plus) is generated from the repository source codes directly. The documentation contains all resources including source code snippets, details on execution of the programs, diagrammatic representation of program flow, and links to external resources where necessary. The documentation also introduces interactive source code with links to documentation for C++ STL library functions used.
|
||||
Click on [Files menu](https://TheAlgorithms.github.io/C-Plus-Plus/files.html) to see the list of all the files documented with the code.
|
||||
|
||||
### Algorithms implemented in C++ (for education)
|
||||
The implementations are for learning purpose. They may be less efficient than the implementation in the standard library.
|
||||
[Documentation of Algorithms in C++](https://thealgorithms.github.io/C-Plus-Plus) by [The Algorithms Contributors](https://github.com/TheAlgorithms/C-Plus-Plus/graphs/contributors) is licensed under [CC BY-SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/?ref=chooser-v1)<br/>
|
||||
<a href="https://creativecommons.org/licenses/by-sa/4.0"><img alt="Creative Commons License" style="height:22px!important;margin-left: 3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/cc.svg" /><img alt="Credit must be given to the creator" style="height:22px!important;margin-left: 3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/by.svg" /><img alt="Adaptations must be shared under the same terms" style="height:22px!important;margin-left: 3px;vertical-align:text-bottom;" src="https://mirrors.creativecommons.org/presskit/icons/sa.svg" /></a>
|
||||
|
||||
### Contribute Guidelines
|
||||
Read our [Contribution Guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md) before you contribute.
|
||||
## Contributions
|
||||
|
||||
As a community developed and maintained repository, we welcome new un-plagiarized quality contributions. Please read our [Contribution Guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md).
|
||||
|
||||
18
ciphers/CMakeLists.txt
Normal file
18
ciphers/CMakeLists.txt
Normal file
@ -0,0 +1,18 @@
|
||||
# If necessary, use the RELATIVE flag, otherwise each source file may be listed
|
||||
# with full pathname. RELATIVE may makes it easier to extract an executable name
|
||||
# automatically.
|
||||
file( GLOB APP_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp )
|
||||
# file( GLOB APP_SOURCES ${CMAKE_SOURCE_DIR}/*.c )
|
||||
# AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR} APP_SOURCES)
|
||||
foreach( testsourcefile ${APP_SOURCES} )
|
||||
# I used a simple string replace, to cut off .cpp.
|
||||
string( REPLACE ".cpp" "" testname ${testsourcefile} )
|
||||
add_executable( ${testname} ${testsourcefile} )
|
||||
|
||||
set_target_properties(${testname} PROPERTIES LINKER_LANGUAGE CXX)
|
||||
if(OpenMP_CXX_FOUND)
|
||||
target_link_libraries(${testname} OpenMP::OpenMP_CXX)
|
||||
endif()
|
||||
install(TARGETS ${testname} DESTINATION "bin/ciphers")
|
||||
|
||||
endforeach( testsourcefile ${APP_SOURCES} )
|
||||
543
ciphers/hill_cipher.cpp
Normal file
543
ciphers/hill_cipher.cpp
Normal file
@ -0,0 +1,543 @@
|
||||
/**
|
||||
* @file hill_cipher.cpp
|
||||
* @brief Implementation of [Hill
|
||||
* cipher](https://en.wikipedia.org/wiki/Hill_cipher) algorithm.
|
||||
*
|
||||
* Program to generate the encryption-decryption key and perform encryption and
|
||||
* decryption of ASCII text using the famous block cipher algorithm. This is a
|
||||
* powerful encryption algorithm that is relatively easy to implement with a
|
||||
* given key. The strength of the algorithm depends on the size of the block
|
||||
* encryption matrix key; the bigger the matrix, the stronger the encryption and
|
||||
* more difficult to break it. However, the important requirement for the matrix
|
||||
* is that:
|
||||
* 1. matrix should be invertible - all inversion conditions should be satisfied
|
||||
* and
|
||||
* 2. its determinant must not have any common factors with the length of
|
||||
* character set
|
||||
* Due to this restriction, most implementations only implement with small 3x3
|
||||
* encryption keys and a small subset of ASCII alphabets.
|
||||
*
|
||||
* In the current implementation, I present to you an implementation for
|
||||
* generating larger encryption keys (I have attempted upto 10x10) and an ASCII
|
||||
* character set of 97 printable characters. Hence, a typical ASCII text file
|
||||
* could be easily encrypted with the module. The larger character set increases
|
||||
* the modulo of cipher and hence the matrix determinants can get very large
|
||||
* very quickly rendering them ill-defined.
|
||||
*
|
||||
* \note This program uses determinant computation using LU decomposition from
|
||||
* the file lu_decomposition.h
|
||||
* \note The matrix generation algorithm is very rudimentary and does not
|
||||
* guarantee an invertible modulus matrix. \todo Better matrix generation
|
||||
* algorithm.
|
||||
*
|
||||
* @author [Krishna Vedala](https://github.com/kvedala)
|
||||
*/
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <ctime>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
#include "../numerical_methods/lu_decomposition.h"
|
||||
|
||||
/**
|
||||
* operator to print a matrix
|
||||
*/
|
||||
template <typename T>
|
||||
static std::ostream &operator<<(std::ostream &out, matrix<T> const &v) {
|
||||
const int width = 15;
|
||||
const char separator = ' ';
|
||||
|
||||
for (size_t row = 0; row < v.size(); row++) {
|
||||
for (size_t col = 0; col < v[row].size(); col++)
|
||||
out << std::left << std::setw(width) << std::setfill(separator)
|
||||
<< v[row][col];
|
||||
out << std::endl;
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
/** \namespace ciphers
|
||||
* \brief Algorithms for encryption and decryption
|
||||
*/
|
||||
namespace ciphers {
|
||||
/** dictionary of characters that can be encrypted and decrypted */
|
||||
static const char *STRKEY =
|
||||
"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789~!@#$%^&"
|
||||
"*()_+`-=[]{}|;':\",./<>?\\\r\n \0";
|
||||
|
||||
/**
|
||||
* @brief Implementation of [Hill
|
||||
* Cipher](https://en.wikipedia.org/wiki/Hill_cipher) algorithm
|
||||
*/
|
||||
class HillCipher {
|
||||
private:
|
||||
/**
|
||||
* @brief Function to generate a random integer in a given interval
|
||||
*
|
||||
* @param a lower limit of interval
|
||||
* @param b upper limit of interval
|
||||
* @tparam T type of output
|
||||
* @return random integer in the interval \f$[a,b)\f$
|
||||
*/
|
||||
template <typename T1, typename T2>
|
||||
static const T2 rand_range(T1 a, T1 b) {
|
||||
// generate random number between 0 and 1
|
||||
long double r = static_cast<long double>(std::rand()) / RAND_MAX;
|
||||
|
||||
// scale and return random number as integer
|
||||
return static_cast<T2>(r * (b - a) + a);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Function overload to fill a matrix with random integers in a given
|
||||
* interval
|
||||
*
|
||||
* @param M pointer to matrix to be filled with random numbers
|
||||
* @param a lower limit of interval
|
||||
* @param b upper limit of interval
|
||||
* @tparam T1 type of input range
|
||||
* @tparam T2 type of matrix
|
||||
* @return determinant of generated random matrix
|
||||
*
|
||||
* @warning There will need to be a balance between the matrix size and the
|
||||
* range of random numbers. If the matrix is large, the range of random
|
||||
* numbers must be small to have a well defined keys. Or if the matrix is
|
||||
* smaller, the random numbers range can be larger. For an 8x8 matrix, range
|
||||
* should be no more than \f$[0,10]\f$
|
||||
*/
|
||||
template <typename T1, typename T2>
|
||||
static double rand_range(matrix<T2> *M, T1 a, T1 b) {
|
||||
for (size_t i = 0; i < M->size(); i++) {
|
||||
for (size_t j = 0; j < M[0][0].size(); j++) {
|
||||
M[0][i][j] = rand_range<T1, T2>(a, b);
|
||||
}
|
||||
}
|
||||
|
||||
return determinant_lu(*M);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Compute
|
||||
* [GCD](https://en.wikipedia.org/wiki/Greatest_common_divisor) of two
|
||||
* integers using Euler's algorithm
|
||||
*
|
||||
* @param a first number
|
||||
* @param b second number
|
||||
* @return GCD of \f$a\f$ and \f$b\f$
|
||||
*/
|
||||
template <typename T>
|
||||
static const T gcd(T a, T b) {
|
||||
if (b > a) // ensure always a < b
|
||||
std::swap(a, b);
|
||||
|
||||
while (b != 0) {
|
||||
T tmp = b;
|
||||
b = a % b;
|
||||
a = tmp;
|
||||
}
|
||||
|
||||
return a;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief helper function to perform vector multiplication with encryption
|
||||
* or decryption matrix
|
||||
*
|
||||
* @param vector vector to multiply
|
||||
* @param key encryption or decryption key matrix
|
||||
* @return corresponding encrypted or decrypted text
|
||||
*/
|
||||
static const std::valarray<uint8_t> mat_mul(
|
||||
const std::valarray<uint8_t> &vector, const matrix<int> &key) {
|
||||
std::valarray<uint8_t> out(vector); // make a copy
|
||||
|
||||
size_t L = std::strlen(STRKEY);
|
||||
|
||||
for (size_t i = 0; i < key.size(); i++) {
|
||||
int tmp = 0;
|
||||
for (size_t j = 0; j < vector.size(); j++) {
|
||||
tmp += key[i][j] * vector[j];
|
||||
}
|
||||
out[i] = static_cast<uint8_t>(tmp % L);
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Get the character at a given index in the ::STRKEY
|
||||
*
|
||||
* @param idx index value
|
||||
* @return character at the index
|
||||
*/
|
||||
static inline char get_idx_char(const uint8_t idx) { return STRKEY[idx]; }
|
||||
|
||||
/**
|
||||
* @brief Get the index of a character in the ::STRKEY
|
||||
*
|
||||
* @param ch character to search
|
||||
* @return index of character
|
||||
*/
|
||||
static inline uint8_t get_char_idx(const char ch) {
|
||||
size_t L = std::strlen(STRKEY);
|
||||
|
||||
for (size_t idx = 0; idx <= L; idx++)
|
||||
if (STRKEY[idx] == ch)
|
||||
return idx;
|
||||
|
||||
std::cerr << __func__ << ":" << __LINE__ << ": (" << ch
|
||||
<< ") Should not reach here!\n";
|
||||
return 0;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Convenience function to perform block cipher operations. The
|
||||
* operations are identical for both encryption and decryption.
|
||||
*
|
||||
* @param text input text to encrypt or decrypt
|
||||
* @param key key for encryption or decryption
|
||||
* @return encrypted/decrypted output
|
||||
*/
|
||||
static const std::string codec(const std::string &text,
|
||||
const matrix<int> &key) {
|
||||
size_t text_len = text.length();
|
||||
size_t key_len = key.size();
|
||||
|
||||
// length of output string must be a multiple of key_len
|
||||
// create output string and initialize with '\0' character
|
||||
size_t L2 = text_len % key_len == 0
|
||||
? text_len
|
||||
: text_len + key_len - (text_len % key_len);
|
||||
std::string coded_text(L2, '\0');
|
||||
|
||||
// temporary array for batch processing
|
||||
int i;
|
||||
#ifdef _OPENMP
|
||||
#pragma parallel omp for private(i)
|
||||
#endif
|
||||
for (i = 0; i < L2 - key_len + 1; i += key_len) {
|
||||
std::valarray<uint8_t> batch_int(key_len);
|
||||
for (size_t j = 0; j < key_len; j++) {
|
||||
batch_int[j] = get_char_idx(text[i + j]);
|
||||
}
|
||||
|
||||
batch_int = mat_mul(batch_int, key);
|
||||
|
||||
for (size_t j = 0; j < key_len; j++) {
|
||||
coded_text[i + j] =
|
||||
STRKEY[batch_int[j]]; // get character at key
|
||||
}
|
||||
}
|
||||
|
||||
return coded_text;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get matrix inverse using Row-transformations. Given matrix must
|
||||
* be a square and non-singular.
|
||||
* \returns inverse matrix
|
||||
**/
|
||||
template <typename T>
|
||||
static matrix<double> get_inverse(matrix<T> const &A) {
|
||||
// Assuming A is square matrix
|
||||
size_t N = A.size();
|
||||
|
||||
matrix<double> inverse(N, std::valarray<double>(N));
|
||||
for (size_t row = 0; row < N; row++) {
|
||||
for (size_t col = 0; col < N; col++) {
|
||||
// create identity matrix
|
||||
inverse[row][col] = (row == col) ? 1.f : 0.f;
|
||||
}
|
||||
}
|
||||
|
||||
if (A.size() != A[0].size()) {
|
||||
std::cerr << "A must be a square matrix!" << std::endl;
|
||||
return inverse;
|
||||
}
|
||||
|
||||
// preallocate a temporary matrix identical to A
|
||||
matrix<double> temp(N, std::valarray<double>(N));
|
||||
for (size_t row = 0; row < N; row++) {
|
||||
for (size_t col = 0; col < N; col++)
|
||||
temp[row][col] = static_cast<double>(A[row][col]);
|
||||
}
|
||||
|
||||
// start transformations
|
||||
for (size_t row = 0; row < N; row++) {
|
||||
for (size_t row2 = row; row2 < N && temp[row][row] == 0; row2++) {
|
||||
// this to ensure diagonal elements are not 0
|
||||
temp[row] = temp[row] + temp[row2];
|
||||
inverse[row] = inverse[row] + inverse[row2];
|
||||
}
|
||||
|
||||
for (size_t col2 = row; col2 < N && temp[row][row] == 0; col2++) {
|
||||
// this to further ensure diagonal elements are not 0
|
||||
for (size_t row2 = 0; row2 < N; row2++) {
|
||||
temp[row2][row] = temp[row2][row] + temp[row2][col2];
|
||||
inverse[row2][row] =
|
||||
inverse[row2][row] + inverse[row2][col2];
|
||||
}
|
||||
}
|
||||
|
||||
if (temp[row][row] == 0) {
|
||||
// Probably a low-rank matrix and hence singular
|
||||
std::cerr << "Low-rank matrix, no inverse!" << std::endl;
|
||||
return inverse;
|
||||
}
|
||||
|
||||
// set diagonal to 1
|
||||
double divisor = temp[row][row];
|
||||
temp[row] = temp[row] / divisor;
|
||||
inverse[row] = inverse[row] / divisor;
|
||||
// Row transformations
|
||||
for (size_t row2 = 0; row2 < N; row2++) {
|
||||
if (row2 == row)
|
||||
continue;
|
||||
double factor = temp[row2][row];
|
||||
temp[row2] = temp[row2] - factor * temp[row];
|
||||
inverse[row2] = inverse[row2] - factor * inverse[row];
|
||||
}
|
||||
}
|
||||
|
||||
return inverse;
|
||||
}
|
||||
|
||||
static int modulo(int a, int b) {
|
||||
int ret = a % b;
|
||||
if (ret < 0)
|
||||
ret += b;
|
||||
return ret;
|
||||
}
|
||||
|
||||
public:
|
||||
/**
|
||||
* @brief Generate encryption matrix of a given size. Larger size matrices
|
||||
* are difficult to generate but provide more security. Important conditions
|
||||
* are:
|
||||
* 1. matrix should be invertible
|
||||
* 2. determinant must not have any common factors with the length of
|
||||
* character key
|
||||
* There is no head-fast way to generate hte matrix under the given
|
||||
* numerical restrictions of the machine but the conditions added achieve
|
||||
* the goals. Bigger the matrix, greater is the probability of the matrix
|
||||
* being ill-defined.
|
||||
*
|
||||
* @param size size of matrix (typically \f$\text{size}\le10\f$)
|
||||
* @param limit1 lower limit of range of random elements (default=0)
|
||||
* @param limit2 upper limit of range of random elements (default=10)
|
||||
* @return Encryption martix
|
||||
*/
|
||||
static matrix<int> generate_encryption_key(size_t size, int limit1 = 0,
|
||||
int limit2 = 10) {
|
||||
matrix<int> encrypt_key(size, std::valarray<int>(size));
|
||||
matrix<int> min_mat = encrypt_key;
|
||||
int mat_determinant = -1; // because matrix has only ints, the
|
||||
// determinant will also be an int
|
||||
int L = std::strlen(STRKEY);
|
||||
|
||||
double dd;
|
||||
do {
|
||||
// keeping the random number range smaller generates better
|
||||
// defined matrices with more ease of cracking
|
||||
dd = rand_range(&encrypt_key, limit1, limit2);
|
||||
mat_determinant = static_cast<int>(dd);
|
||||
|
||||
if (mat_determinant < 0)
|
||||
mat_determinant = (mat_determinant % L);
|
||||
} while (std::abs(dd) > 1e3 || // while ill-defined
|
||||
dd < 0.1 || // while singular or negative determinant
|
||||
!std::isfinite(dd) || // while determinant is not finite
|
||||
gcd(mat_determinant, L) != 1); // while no common factors
|
||||
// std::cout <<
|
||||
|
||||
return encrypt_key;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Generate decryption matrix from an encryption matrix key.
|
||||
*
|
||||
* @param encrypt_key encryption key for which to create a decrypt key
|
||||
* @return Decryption martix
|
||||
*/
|
||||
static matrix<int> generate_decryption_key(matrix<int> const &encrypt_key) {
|
||||
size_t size = encrypt_key.size();
|
||||
int L = std::strlen(STRKEY);
|
||||
|
||||
matrix<int> decrypt_key(size, std::valarray<int>(size));
|
||||
int det_encrypt = static_cast<int>(determinant_lu(encrypt_key));
|
||||
|
||||
int mat_determinant = det_encrypt < 0 ? det_encrypt % L : det_encrypt;
|
||||
|
||||
matrix<double> tmp_inverse = get_inverse(encrypt_key);
|
||||
double d2 = determinant_lu(decrypt_key);
|
||||
|
||||
// find co-prime factor for inversion
|
||||
int det_inv = -1;
|
||||
for (int i = 0; i < L; i++) {
|
||||
if (modulo(mat_determinant * i, L) == 1) {
|
||||
det_inv = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (det_inv == -1) {
|
||||
std::cerr << "Could not find a co-prime for inversion\n";
|
||||
std::exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
mat_determinant = det_inv * det_encrypt;
|
||||
|
||||
// perform modular inverse of encryption matrix
|
||||
int i;
|
||||
#ifdef _OPENMP
|
||||
#pragma parallel omp for private(i)
|
||||
#endif
|
||||
for (i = 0; i < size; i++) {
|
||||
for (int j = 0; j < size; j++) {
|
||||
int temp = std::round(tmp_inverse[i][j] * mat_determinant);
|
||||
decrypt_key[i][j] = modulo(temp, L);
|
||||
}
|
||||
}
|
||||
return decrypt_key;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Generate encryption and decryption key pair
|
||||
*
|
||||
* @param size size of matrix key (typically \f$\text{size}\le10\f$)
|
||||
* @param limit1 lower limit of range of random elements (default=0)
|
||||
* @param limit2 upper limit of range of random elements (default=10)
|
||||
* @return std::pair<matrix<int>, matrix<int>> encryption and decryption
|
||||
* keys as a pair
|
||||
*
|
||||
* @see ::generate_encryption_key
|
||||
*/
|
||||
static std::pair<matrix<int>, matrix<int>> generate_keys(size_t size,
|
||||
int limit1 = 0,
|
||||
int limit2 = 10) {
|
||||
matrix<int> encrypt_key = generate_encryption_key(size);
|
||||
matrix<int> decrypt_key = generate_decryption_key(encrypt_key);
|
||||
double det2 = determinant_lu(decrypt_key);
|
||||
while (std::abs(det2) < 0.1 || std::abs(det2) > 1e3) {
|
||||
encrypt_key = generate_encryption_key(size, limit1, limit2);
|
||||
decrypt_key = generate_decryption_key(encrypt_key);
|
||||
det2 = determinant_lu(decrypt_key);
|
||||
}
|
||||
return std::make_pair(encrypt_key, decrypt_key);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Encrypt a given text using a given key
|
||||
*
|
||||
* @param text string to encrypt
|
||||
* @param encrypt_key key for encryption
|
||||
* @return encrypted text
|
||||
*/
|
||||
static const std::string encrypt_text(const std::string &text,
|
||||
const matrix<int> &encrypt_key) {
|
||||
return codec(text, encrypt_key);
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Decrypt a given text using a given key
|
||||
*
|
||||
* @param text string to decrypt
|
||||
* @param decrypt_key key for decryption
|
||||
* @return decrypted text
|
||||
*/
|
||||
static const std::string decrypt_text(const std::string &text,
|
||||
const matrix<int> &decrypt_key) {
|
||||
return codec(text, decrypt_key);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace ciphers
|
||||
|
||||
/**
|
||||
* @brief Self test 1 - using 3x3 randomly generated key
|
||||
*
|
||||
* @param text string to encrypt and decrypt
|
||||
*/
|
||||
void test1(const std::string &text) {
|
||||
// std::string text = "Hello world!";
|
||||
std::cout << "======Test 1 (3x3 key) ======\nOriginal text:\n\t" << text
|
||||
<< std::endl;
|
||||
|
||||
std::pair<matrix<int>, matrix<int>> p =
|
||||
ciphers::HillCipher::generate_keys(3, 0, 100);
|
||||
matrix<int> ekey = p.first;
|
||||
matrix<int> dkey = p.second;
|
||||
|
||||
// matrix<int> ekey = {{22, 28, 25}, {5, 26, 15}, {14, 18, 9}};
|
||||
// std::cout << "Encryption key: \n" << ekey;
|
||||
std::string gibberish = ciphers::HillCipher::encrypt_text(text, ekey);
|
||||
std::cout << "Encrypted text:\n\t" << gibberish << std::endl;
|
||||
|
||||
// matrix<int> dkey = ciphers::HillCipher::generate_decryption_key(ekey);
|
||||
// std::cout << "Decryption key: \n" << dkey;
|
||||
std::string txt_back = ciphers::HillCipher::decrypt_text(gibberish, dkey);
|
||||
std::cout << "Reconstruct text:\n\t" << txt_back << std::endl;
|
||||
|
||||
std::ofstream out_file("hill_cipher_test1.txt");
|
||||
out_file << "Block size: " << ekey.size() << "\n";
|
||||
out_file << "Encryption Key:\n" << ekey;
|
||||
out_file << "\nDecryption Key:\n" << dkey;
|
||||
out_file.close();
|
||||
|
||||
assert(txt_back == text);
|
||||
std::cout << "Passed :)\n";
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Self test 2 - using 8x8 randomly generated key
|
||||
*
|
||||
* @param text string to encrypt and decrypt
|
||||
*/
|
||||
void test2(const std::string &text) {
|
||||
// std::string text = "Hello world!";
|
||||
std::cout << "======Test 2 (8x8 key) ======\nOriginal text:\n\t" << text
|
||||
<< std::endl;
|
||||
|
||||
std::pair<matrix<int>, matrix<int>> p =
|
||||
ciphers::HillCipher::generate_keys(8, 0, 3);
|
||||
matrix<int> ekey = p.first;
|
||||
matrix<int> dkey = p.second;
|
||||
|
||||
std::string gibberish = ciphers::HillCipher::encrypt_text(text, ekey);
|
||||
std::cout << "Encrypted text:\n\t" << gibberish << std::endl;
|
||||
|
||||
std::string txt_back = ciphers::HillCipher::decrypt_text(gibberish, dkey);
|
||||
std::cout << "Reconstruct text:\n\t" << txt_back << std::endl;
|
||||
|
||||
std::ofstream out_file("hill_cipher_test2.txt");
|
||||
out_file << "Block size: " << ekey.size() << "\n";
|
||||
out_file << "Encryption Key:\n" << ekey;
|
||||
out_file << "\nDecryption Key:\n" << dkey;
|
||||
out_file.close();
|
||||
|
||||
assert(txt_back.compare(0, text.size(), text) == 0);
|
||||
std::cout << "Passed :)\n";
|
||||
}
|
||||
|
||||
/** Main function */
|
||||
int main() {
|
||||
std::srand(std::time(nullptr));
|
||||
std::cout << "Key dictionary: (" << std::strlen(ciphers::STRKEY) << ")\n\t"
|
||||
<< ciphers::STRKEY << "\n";
|
||||
|
||||
std::string text = "This is a simple text with numb3r5 and exclamat!0n.";
|
||||
|
||||
test1(text);
|
||||
test2(text);
|
||||
|
||||
return 0;
|
||||
}
|
||||
@ -14,17 +14,17 @@ struct node {
|
||||
node *right;
|
||||
};
|
||||
|
||||
struct queue {
|
||||
struct Queue {
|
||||
node *t[100];
|
||||
int front;
|
||||
int rear;
|
||||
};
|
||||
|
||||
queue q;
|
||||
Queue queue;
|
||||
|
||||
void enqueue(node *n) { q.t[q.rear++] = n; }
|
||||
void enqueue(node *n) { queue.t[queue.rear++] = n; }
|
||||
|
||||
node *dequeue() { return (q.t[q.front++]); }
|
||||
node *dequeue() { return (queue.t[queue.front++]); }
|
||||
|
||||
void Insert(node *n, int x) {
|
||||
if (x < n->val) {
|
||||
@ -123,8 +123,8 @@ void Post(node *n) {
|
||||
}
|
||||
|
||||
int main() {
|
||||
q.front = 0;
|
||||
q.rear = 0;
|
||||
queue.front = 0;
|
||||
queue.rear = 0;
|
||||
int value;
|
||||
int ch;
|
||||
node *root = new node;
|
||||
|
||||
@ -1,18 +1,27 @@
|
||||
/* This class specifies the basic operation on a stack as a linked list */
|
||||
/**
|
||||
* @file stack.h
|
||||
* @author danghai
|
||||
* @brief This class specifies the basic operation on a stack as a linked list
|
||||
**/
|
||||
#ifndef DATA_STRUCTURES_STACK_H_
|
||||
#define DATA_STRUCTURES_STACK_H_
|
||||
|
||||
#include <cassert>
|
||||
#include <iostream>
|
||||
|
||||
/* Definition of the node */
|
||||
/** Definition of the node as a linked-list
|
||||
* \tparam Type type of data nodes of the linked list should contain
|
||||
*/
|
||||
template <class Type>
|
||||
struct node {
|
||||
Type data;
|
||||
node<Type> *next;
|
||||
Type data; ///< data at current node
|
||||
node<Type> *next; ///< pointer to the next ::node instance
|
||||
};
|
||||
|
||||
/* Definition of the stack class */
|
||||
/** Definition of the stack class
|
||||
* \tparam Type type of data nodes of the linked list in the stack should
|
||||
* contain
|
||||
*/
|
||||
template <class Type>
|
||||
class stack {
|
||||
public:
|
||||
@ -20,7 +29,7 @@ class stack {
|
||||
void display() {
|
||||
node<Type> *current = stackTop;
|
||||
std::cout << "Top --> ";
|
||||
while (current != NULL) {
|
||||
while (current != nullptr) {
|
||||
std::cout << current->data << " ";
|
||||
current = current->next;
|
||||
}
|
||||
@ -30,15 +39,45 @@ class stack {
|
||||
|
||||
/** Default constructor*/
|
||||
stack() {
|
||||
stackTop = NULL;
|
||||
stackTop = nullptr;
|
||||
size = 0;
|
||||
}
|
||||
|
||||
/** Copy constructor*/
|
||||
explicit stack(const stack<Type> &otherStack) {
|
||||
node<Type> *newNode, *current, *last;
|
||||
|
||||
/* If stack is no empty, make it empty */
|
||||
if (stackTop != nullptr) {
|
||||
stackTop = nullptr;
|
||||
}
|
||||
if (otherStack.stackTop == nullptr) {
|
||||
stackTop = nullptr;
|
||||
} else {
|
||||
current = otherStack.stackTop;
|
||||
stackTop = new node<Type>;
|
||||
stackTop->data = current->data;
|
||||
stackTop->next = nullptr;
|
||||
last = stackTop;
|
||||
current = current->next;
|
||||
/* Copy the remaining stack */
|
||||
while (current != nullptr) {
|
||||
newNode = new node<Type>;
|
||||
newNode->data = current->data;
|
||||
newNode->next = nullptr;
|
||||
last->next = newNode;
|
||||
last = newNode;
|
||||
current = current->next;
|
||||
}
|
||||
}
|
||||
size = otherStack.size;
|
||||
}
|
||||
|
||||
/** Destructor */
|
||||
~stack() {}
|
||||
|
||||
/** Determine whether the stack is empty */
|
||||
bool isEmptyStack() { return (stackTop == NULL); }
|
||||
bool isEmptyStack() { return (stackTop == nullptr); }
|
||||
|
||||
/** Add new item to the stack */
|
||||
void push(Type item) {
|
||||
@ -52,7 +91,7 @@ class stack {
|
||||
|
||||
/** Return the top element of the stack */
|
||||
Type top() {
|
||||
assert(stackTop != NULL);
|
||||
assert(stackTop != nullptr);
|
||||
return stackTop->data;
|
||||
}
|
||||
|
||||
@ -70,30 +109,30 @@ class stack {
|
||||
}
|
||||
|
||||
/** Clear stack */
|
||||
void clear() { stackTop = NULL; }
|
||||
void clear() { stackTop = nullptr; }
|
||||
|
||||
/** Overload "=" the assignment operator */
|
||||
stack<Type> &operator=(const stack<Type> &otherStack) {
|
||||
node<Type> *newNode, *current, *last;
|
||||
|
||||
/* If stack is no empty, make it empty */
|
||||
if (stackTop != NULL) {
|
||||
stackTop = NULL;
|
||||
if (stackTop != nullptr) {
|
||||
stackTop = nullptr;
|
||||
}
|
||||
if (otherStack.stackTop == NULL) {
|
||||
stackTop = NULL;
|
||||
if (otherStack.stackTop == nullptr) {
|
||||
stackTop = nullptr;
|
||||
} else {
|
||||
current = otherStack.stackTop;
|
||||
stackTop = new node<Type>;
|
||||
stackTop->data = current->data;
|
||||
stackTop->next = NULL;
|
||||
stackTop->next = nullptr;
|
||||
last = stackTop;
|
||||
current = current->next;
|
||||
/* Copy the remaining stack */
|
||||
while (current != NULL) {
|
||||
while (current != nullptr) {
|
||||
newNode = new node<Type>;
|
||||
newNode->data = current->data;
|
||||
newNode->next = NULL;
|
||||
newNode->next = nullptr;
|
||||
last->next = newNode;
|
||||
last = newNode;
|
||||
current = current->next;
|
||||
@ -105,7 +144,7 @@ class stack {
|
||||
|
||||
private:
|
||||
node<Type> *stackTop; /**< Pointer to the stack */
|
||||
int size;
|
||||
int size; ///< size of stack
|
||||
};
|
||||
|
||||
#endif // DATA_STRUCTURES_STACK_H_
|
||||
|
||||
@ -1,37 +1,38 @@
|
||||
#include <iostream>
|
||||
|
||||
int *stack;
|
||||
int top = 0, stack_size;
|
||||
int stack_idx = 0, stack_size;
|
||||
|
||||
void push(int x) {
|
||||
if (top == stack_size) {
|
||||
if (stack_idx == stack_size) {
|
||||
std::cout << "\nOverflow";
|
||||
} else {
|
||||
stack[top++] = x;
|
||||
stack[stack_idx++] = x;
|
||||
}
|
||||
}
|
||||
|
||||
void pop() {
|
||||
if (top == 0) {
|
||||
if (stack_idx == 0) {
|
||||
std::cout << "\nUnderflow";
|
||||
} else {
|
||||
std::cout << "\n" << stack[--top] << " deleted";
|
||||
std::cout << "\n" << stack[--stack_idx] << " deleted";
|
||||
}
|
||||
}
|
||||
|
||||
void show() {
|
||||
for (int i = 0; i < top; i++) {
|
||||
for (int i = 0; i < stack_idx; i++) {
|
||||
std::cout << stack[i] << "\n";
|
||||
}
|
||||
}
|
||||
|
||||
void topmost() { std::cout << "\nTopmost element: " << stack[top - 1]; }
|
||||
void topmost() { std::cout << "\nTopmost element: " << stack[stack_idx - 1]; }
|
||||
int main() {
|
||||
std::cout << "\nEnter stack_size of stack : ";
|
||||
std::cin >> stack_size;
|
||||
stack = new int[stack_size];
|
||||
int ch, x;
|
||||
do {
|
||||
std::cout << "\n0. Exit";
|
||||
std::cout << "\n1. Push";
|
||||
std::cout << "\n2. Pop";
|
||||
std::cout << "\n3. Print";
|
||||
@ -51,5 +52,7 @@ int main() {
|
||||
}
|
||||
} while (ch != 0);
|
||||
|
||||
delete[] stack;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -1,35 +1,34 @@
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
struct node {
|
||||
int val;
|
||||
node *next;
|
||||
};
|
||||
|
||||
node *top;
|
||||
node *top_var;
|
||||
|
||||
void push(int x) {
|
||||
node *n = new node;
|
||||
n->val = x;
|
||||
n->next = top;
|
||||
top = n;
|
||||
n->next = top_var;
|
||||
top_var = n;
|
||||
}
|
||||
|
||||
void pop() {
|
||||
if (top == NULL) {
|
||||
cout << "\nUnderflow";
|
||||
if (top_var == NULL) {
|
||||
std::cout << "\nUnderflow";
|
||||
} else {
|
||||
node *t = top;
|
||||
cout << "\n" << t->val << " deleted";
|
||||
top = top->next;
|
||||
node *t = top_var;
|
||||
std::cout << "\n" << t->val << " deleted";
|
||||
top_var = top_var->next;
|
||||
delete t;
|
||||
}
|
||||
}
|
||||
|
||||
void show() {
|
||||
node *t = top;
|
||||
node *t = top_var;
|
||||
while (t != NULL) {
|
||||
cout << t->val << "\n";
|
||||
std::cout << t->val << "\n";
|
||||
t = t->next;
|
||||
}
|
||||
}
|
||||
@ -37,14 +36,14 @@ void show() {
|
||||
int main() {
|
||||
int ch, x;
|
||||
do {
|
||||
cout << "\n1. Push";
|
||||
cout << "\n2. Pop";
|
||||
cout << "\n3. Print";
|
||||
cout << "\nEnter Your Choice : ";
|
||||
cin >> ch;
|
||||
std::cout << "\n1. Push";
|
||||
std::cout << "\n2. Pop";
|
||||
std::cout << "\n3. Print";
|
||||
std::cout << "\nEnter Your Choice : ";
|
||||
std::cin >> ch;
|
||||
if (ch == 1) {
|
||||
cout << "\nInsert : ";
|
||||
cin >> x;
|
||||
std::cout << "\nInsert : ";
|
||||
std::cin >> x;
|
||||
push(x);
|
||||
} else if (ch == 2) {
|
||||
pop();
|
||||
|
||||
@ -9,6 +9,8 @@
|
||||
*
|
||||
* \author [Krishna Vedala](https://github.com/kvedala)
|
||||
*/
|
||||
#include <cassert>
|
||||
#include <cmath> // for std::abs
|
||||
#include <iomanip> // for print formatting
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
@ -352,10 +354,62 @@ std::vector<float> predict_OLS_regressor(std::vector<std::vector<T>> const &X,
|
||||
return result;
|
||||
}
|
||||
|
||||
/** Self test checks */
|
||||
void ols_test() {
|
||||
int F = 3, N = 5;
|
||||
|
||||
/* test function = x^2 -5 */
|
||||
std::cout << "Test 1 (quadratic function)....";
|
||||
// create training data set with features = x, x^2, x^3
|
||||
std::vector<std::vector<float>> data1(
|
||||
{{-5, 25, -125}, {-1, 1, -1}, {0, 0, 0}, {1, 1, 1}, {6, 36, 216}});
|
||||
// create corresponding outputs
|
||||
std::vector<float> Y1({20, -4, -5, -4, 31});
|
||||
// perform regression modelling
|
||||
std::vector<float> beta1 = fit_OLS_regressor(data1, Y1);
|
||||
// create test data set with same features = x, x^2, x^3
|
||||
std::vector<std::vector<float>> test_data1(
|
||||
{{-2, 4, -8}, {2, 4, 8}, {-10, 100, -1000}, {10, 100, 1000}});
|
||||
// expected regression outputs
|
||||
std::vector<float> expected1({-1, -1, 95, 95});
|
||||
// predicted regression outputs
|
||||
std::vector<float> out1 = predict_OLS_regressor(test_data1, beta1);
|
||||
// compare predicted results are within +-0.01 limit of expected
|
||||
for (size_t rows = 0; rows < out1.size(); rows++)
|
||||
assert(std::abs(out1[rows] - expected1[rows]) < 0.01);
|
||||
std::cout << "passed\n";
|
||||
|
||||
/* test function = x^3 + x^2 - 100 */
|
||||
std::cout << "Test 2 (cubic function)....";
|
||||
// create training data set with features = x, x^2, x^3
|
||||
std::vector<std::vector<float>> data2(
|
||||
{{-5, 25, -125}, {-1, 1, -1}, {0, 0, 0}, {1, 1, 1}, {6, 36, 216}});
|
||||
// create corresponding outputs
|
||||
std::vector<float> Y2({-200, -100, -100, 98, 152});
|
||||
// perform regression modelling
|
||||
std::vector<float> beta2 = fit_OLS_regressor(data2, Y2);
|
||||
// create test data set with same features = x, x^2, x^3
|
||||
std::vector<std::vector<float>> test_data2(
|
||||
{{-2, 4, -8}, {2, 4, 8}, {-10, 100, -1000}, {10, 100, 1000}});
|
||||
// expected regression outputs
|
||||
std::vector<float> expected2({-104, -88, -1000, 1000});
|
||||
// predicted regression outputs
|
||||
std::vector<float> out2 = predict_OLS_regressor(test_data2, beta2);
|
||||
// compare predicted results are within +-0.01 limit of expected
|
||||
for (size_t rows = 0; rows < out2.size(); rows++)
|
||||
assert(std::abs(out2[rows] - expected2[rows]) < 0.01);
|
||||
std::cout << "passed\n";
|
||||
|
||||
std::cout << std::endl; // ensure test results are displayed on screen
|
||||
// (flush stdout)
|
||||
}
|
||||
|
||||
/**
|
||||
* main function
|
||||
*/
|
||||
int main() {
|
||||
ols_test();
|
||||
|
||||
size_t N, F;
|
||||
|
||||
std::cout << "Enter number of features: ";
|
||||
@ -369,7 +423,7 @@ int main() {
|
||||
std::vector<float> Y(N);
|
||||
|
||||
std::cout
|
||||
<< "Enter training data. Per sample, provide features ad one output."
|
||||
<< "Enter training data. Per sample, provide features and one output."
|
||||
<< std::endl;
|
||||
|
||||
for (size_t rows = 0; rows < N; rows++) {
|
||||
80
math/armstrong_number.cpp
Normal file
80
math/armstrong_number.cpp
Normal file
@ -0,0 +1,80 @@
|
||||
/**
|
||||
* @file
|
||||
* \brief Program to check if a number is an [Armstrong/Narcissistic
|
||||
* number](https://en.wikipedia.org/wiki/Narcissistic_number) in decimal system.
|
||||
*
|
||||
* \details
|
||||
* Armstrong number or [Narcissistic
|
||||
* number](https://en.wikipedia.org/wiki/Narcissistic_number) is a number that
|
||||
* is the sum of its own digits raised to the power of the number of digits.
|
||||
* @author iamnambiar
|
||||
*/
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
|
||||
/**
|
||||
* Function to calculate the total number of digits in the number.
|
||||
* @param num Number
|
||||
* @return Total number of digits.
|
||||
*/
|
||||
int number_of_digits(int num) {
|
||||
int total_digits = 0;
|
||||
while (num > 0) {
|
||||
num = num / 10;
|
||||
++total_digits;
|
||||
}
|
||||
return total_digits;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function to check whether the number is armstrong number or not.
|
||||
* @param num Number
|
||||
* @return `true` if the number is armstrong.
|
||||
* @return `false` if the number is not armstrong.
|
||||
*/
|
||||
bool is_armstrong(int number) {
|
||||
// If the number is less than 0, then it is not a armstrong number.
|
||||
if (number < 0) {
|
||||
return false;
|
||||
}
|
||||
int sum = 0;
|
||||
int temp = number;
|
||||
// Finding the total number of digits in the number
|
||||
int total_digits = number_of_digits(number);
|
||||
while (temp > 0) {
|
||||
int rem = temp % 10;
|
||||
// Finding each digit raised to the power total digit and add it to the
|
||||
// total sum
|
||||
sum = sum + std::pow(rem, total_digits);
|
||||
temp = temp / 10;
|
||||
}
|
||||
return number == sum;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function for testing the is_armstrong() function
|
||||
* with all the test cases.
|
||||
*/
|
||||
void test() {
|
||||
// is_armstrong(370) returns true.
|
||||
assert(is_armstrong(370) == true);
|
||||
// is_armstrong(225) returns false.
|
||||
assert(is_armstrong(225) == false);
|
||||
// is_armstrong(-23) returns false.
|
||||
assert(is_armstrong(-23) == false);
|
||||
// is_armstrong(153) returns true.
|
||||
assert(is_armstrong(153) == true);
|
||||
// is_armstrong(0) returns true.
|
||||
assert(is_armstrong(0) == true);
|
||||
// is_armstrong(12) returns false.
|
||||
assert(is_armstrong(12) == false);
|
||||
}
|
||||
|
||||
/**
|
||||
* Main Function
|
||||
*/
|
||||
int main() {
|
||||
test();
|
||||
return 0;
|
||||
}
|
||||
72
math/check_amicable_pair.cpp
Normal file
72
math/check_amicable_pair.cpp
Normal file
@ -0,0 +1,72 @@
|
||||
/**
|
||||
*
|
||||
* @file
|
||||
* \brief A C++ Program to check whether a pair of number is [amicable
|
||||
* pair](https://en.wikipedia.org/wiki/Amicable_numbers) or not.
|
||||
*
|
||||
* \details
|
||||
* Amicable Pair are two positive integers such that sum of the proper divisor
|
||||
* of each number is equal to the other number.
|
||||
* @author iamnambiar
|
||||
*/
|
||||
#include <cassert>
|
||||
#include <iostream>
|
||||
|
||||
/**
|
||||
* Function to calculate the sum of all the proper divisor
|
||||
* of an integer.
|
||||
* @param num First number.
|
||||
* @return Sum of the proper divisor of the number.
|
||||
*/
|
||||
int sum_of_divisor(int num) {
|
||||
// Variable to store the sum of all proper divisors.
|
||||
int sum = 0;
|
||||
// Below loop condition helps to reduce Time complexity by a factor of
|
||||
// square root of the number.
|
||||
for (int div = 2; div * div <= num; ++div) {
|
||||
// Check 'div' is divisor of 'num'.
|
||||
if (num % div == 0) {
|
||||
// If both divisor are same, add once to 'sum'
|
||||
if (div == (num / div)) {
|
||||
sum += div;
|
||||
} else {
|
||||
// If both divisor are not the same, add both to 'sum'.
|
||||
sum += (div + (num / div));
|
||||
}
|
||||
}
|
||||
}
|
||||
return sum + 1;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function to check whether the pair is amicable or not.
|
||||
* @param x First number.
|
||||
* @param y Second number.
|
||||
* @return `true` if the pair is amicable
|
||||
* @return `false` if the pair is not amicable
|
||||
*/
|
||||
bool are_amicable(int x, int y) {
|
||||
return (sum_of_divisor(x) == y) && (sum_of_divisor(y) == x);
|
||||
}
|
||||
|
||||
/**
|
||||
* Function for testing the is_amicable() with
|
||||
* all the test cases.
|
||||
*/
|
||||
void test() {
|
||||
// are_amicable(220, 284) returns true.
|
||||
assert(are_amicable(220, 284) == true);
|
||||
// are_amicable(6232, 6368) returns true.
|
||||
assert(are_amicable(6368, 6232) == true);
|
||||
// are_amicable(458, 232) returns false.
|
||||
assert(are_amicable(458, 232) == false);
|
||||
}
|
||||
|
||||
/**
|
||||
* Main Function
|
||||
*/
|
||||
int main() {
|
||||
test();
|
||||
std::cout << "Assertion Success." << std::endl;
|
||||
return 0;
|
||||
}
|
||||
256
math/complex_numbers.cpp
Normal file
256
math/complex_numbers.cpp
Normal file
@ -0,0 +1,256 @@
|
||||
/**
|
||||
* @author tjgurwara99
|
||||
* @file
|
||||
*
|
||||
* A basic implementation of Complex Number field as a class with operators
|
||||
* overloaded to accommodate (mathematical) field operations.
|
||||
*/
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <ctime>
|
||||
#include <iostream>
|
||||
#include <stdexcept>
|
||||
|
||||
/**
|
||||
* Class Complex to represent complex numbers as a field.
|
||||
*/
|
||||
class Complex {
|
||||
// The real value of the complex number
|
||||
double re;
|
||||
// The imaginary value of the complex number
|
||||
double im;
|
||||
|
||||
public:
|
||||
/**
|
||||
* Complex Constructor which initialises the complex number which takes two
|
||||
* arguments.
|
||||
* @param x If the third parameter is 'true' then this x is the absolute
|
||||
* value of the complex number, if the third parameter is 'false' then this
|
||||
* x is the real value of the complex number (optional).
|
||||
* @param y If the third parameter is 'true' then this y is the argument of
|
||||
* the complex number, if the third parameter is 'false' then this y is the
|
||||
* imaginary value of the complex number (optional).
|
||||
* @param is_polar 'false' by default. If we want to initialise our complex
|
||||
* number using polar form then set this to true, otherwise set it to false
|
||||
* to use initialiser which initialises real and imaginary values using the
|
||||
* first two parameters (optional).
|
||||
*/
|
||||
explicit Complex(double x = 0.f, double y = 0.f, bool is_polar = false) {
|
||||
if (!is_polar) {
|
||||
re = x;
|
||||
im = y;
|
||||
return;
|
||||
}
|
||||
|
||||
re = x * std::cos(y);
|
||||
im = x * std::sin(y);
|
||||
}
|
||||
|
||||
/**
|
||||
* Copy Constructor
|
||||
* @param other The other number to equate our number to.
|
||||
*/
|
||||
Complex(const Complex &other) : re(other.real()), im(other.imag()) {}
|
||||
|
||||
/**
|
||||
* Member function (getter) to access the class' re value.
|
||||
*/
|
||||
double real() const { return this->re; }
|
||||
|
||||
/**
|
||||
* Member function (getter) to access the class' im value.
|
||||
*/
|
||||
double imag() const { return this->im; }
|
||||
|
||||
/**
|
||||
* Member function to which gives the absolute value (modulus) of our
|
||||
* complex number
|
||||
* @return \f$ \sqrt{z \dot \bar{z}} \f$ where \f$ z \f$ is our complex
|
||||
* number.
|
||||
*/
|
||||
double abs() const {
|
||||
return std::sqrt(this->re * this->re + this->im * this->im);
|
||||
}
|
||||
|
||||
/**
|
||||
* Member function which gives the argument of our complex number.
|
||||
* @return Argument of our Complex number in radians.
|
||||
*/
|
||||
double arg() const { return std::atan2(this->im, this->re); }
|
||||
|
||||
/**
|
||||
* Operator overload to be able to add two complex numbers.
|
||||
* @param other The other number that is added to the current number.
|
||||
* @return result current number plus other number
|
||||
*/
|
||||
Complex operator+(const Complex &other) {
|
||||
Complex result(this->re + other.re, this->im + other.im);
|
||||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* Operator overload to be able to subtract two complex numbers.
|
||||
* @param other The other number being subtracted from the current number.
|
||||
* @return result current number subtract other number
|
||||
*/
|
||||
Complex operator-(const Complex &other) {
|
||||
Complex result(this->re - other.re, this->im - other.im);
|
||||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* Operator overload to be able to multiple two complex numbers.
|
||||
* @param other The other number to multiply the current number to.
|
||||
* @return result current number times other number.
|
||||
*/
|
||||
Complex operator*(const Complex &other) {
|
||||
Complex result(this->re * other.re - this->im * other.im,
|
||||
this->re * other.im + this->im * other.re);
|
||||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* Operator overload of the BITWISE NOT which gives us the conjugate of our
|
||||
* complex number. NOTE: This is overloading the BITWISE operator but its
|
||||
* not a BITWISE operation in this definition.
|
||||
* @return result The conjugate of our complex number.
|
||||
*/
|
||||
Complex operator~() const {
|
||||
Complex result(this->re, -(this->im));
|
||||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* Operator overload to be able to divide two complex numbers. This function
|
||||
* would throw an exception if the other number is zero.
|
||||
* @param other The other number we divide our number by.
|
||||
* @return result Current number divided by other number.
|
||||
*/
|
||||
Complex operator/(const Complex &other) {
|
||||
Complex result = *this * ~other;
|
||||
double denominator =
|
||||
other.real() * other.real() + other.imag() * other.imag();
|
||||
if (denominator != 0) {
|
||||
result = Complex(result.real() / denominator,
|
||||
result.imag() / denominator);
|
||||
return result;
|
||||
} else {
|
||||
throw std::invalid_argument("Undefined Value");
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Operator overload to be able to copy RHS instance of Complex to LHS
|
||||
* instance of Complex
|
||||
*/
|
||||
const Complex &operator=(const Complex &other) {
|
||||
this->re = other.real();
|
||||
this->im = other.imag();
|
||||
return *this;
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* Logical Equal overload for our Complex class.
|
||||
* @param a Left hand side of our expression
|
||||
* @param b Right hand side of our expression
|
||||
* @return 'True' If real and imaginary parts of a and b are same
|
||||
* @return 'False' Otherwise.
|
||||
*/
|
||||
bool operator==(const Complex &a, const Complex &b) {
|
||||
return a.real() == b.real() && a.imag() == b.imag();
|
||||
}
|
||||
|
||||
/**
|
||||
* Overloaded insersion operator to accommodate the printing of our complex
|
||||
* number in their standard form.
|
||||
* @param os The console stream
|
||||
* @param num The complex number.
|
||||
*/
|
||||
std::ostream &operator<<(std::ostream &os, const Complex &num) {
|
||||
os << "(" << num.real();
|
||||
if (num.imag() < 0) {
|
||||
os << " - " << -num.imag();
|
||||
} else {
|
||||
os << " + " << num.imag();
|
||||
}
|
||||
os << "i)";
|
||||
return os;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function to get random numbers to generate our complex numbers for test
|
||||
*/
|
||||
double get_rand() { return (std::rand() % 100 - 50) / 100.f; }
|
||||
|
||||
/**
|
||||
* Tests Function
|
||||
*/
|
||||
void tests() {
|
||||
std::srand(std::time(nullptr));
|
||||
double x1 = get_rand(), y1 = get_rand(), x2 = get_rand(), y2 = get_rand();
|
||||
Complex num1(x1, y1), num2(x2, y2);
|
||||
std::complex<double> cnum1(x1, y1), cnum2(x2, y2);
|
||||
Complex result;
|
||||
std::complex<double> expected;
|
||||
// Test for addition
|
||||
result = num1 + num2;
|
||||
expected = cnum1 + cnum2;
|
||||
assert(((void)"1 + 1i + 1 + 1i is equal to 2 + 2i but the addition doesn't "
|
||||
"add up \n",
|
||||
(result.real() == expected.real() &&
|
||||
result.imag() == expected.imag())));
|
||||
std::cout << "First test passes." << std::endl;
|
||||
// Test for subtraction
|
||||
result = num1 - num2;
|
||||
expected = cnum1 - cnum2;
|
||||
assert(((void)"1 + 1i - 1 - 1i is equal to 0 but the program says "
|
||||
"otherwise. \n",
|
||||
(result.real() == expected.real() &&
|
||||
result.imag() == expected.imag())));
|
||||
std::cout << "Second test passes." << std::endl;
|
||||
// Test for multiplication
|
||||
result = num1 * num2;
|
||||
expected = cnum1 * cnum2;
|
||||
assert(((void)"(1 + 1i) * (1 + 1i) is equal to 2i but the program says "
|
||||
"otherwise. \n",
|
||||
(result.real() == expected.real() &&
|
||||
result.imag() == expected.imag())));
|
||||
std::cout << "Third test passes." << std::endl;
|
||||
// Test for division
|
||||
result = num1 / num2;
|
||||
expected = cnum1 / cnum2;
|
||||
assert(((void)"(1 + 1i) / (1 + 1i) is equal to 1 but the program says "
|
||||
"otherwise.\n",
|
||||
(result.real() == expected.real() &&
|
||||
result.imag() == expected.imag())));
|
||||
std::cout << "Fourth test passes." << std::endl;
|
||||
// Test for conjugates
|
||||
result = ~num1;
|
||||
expected = std::conj(cnum1);
|
||||
assert(((void)"(1 + 1i) has a conjugate which is equal to (1 - 1i) but the "
|
||||
"program says otherwise.\n",
|
||||
(result.real() == expected.real() &&
|
||||
result.imag() == expected.imag())));
|
||||
std::cout << "Fifth test passes.\n";
|
||||
// Test for Argument of our complex number
|
||||
assert(((void)"(1 + 1i) has argument PI / 4 but the program differs from "
|
||||
"the std::complex result.\n",
|
||||
(num1.arg() == std::arg(cnum1))));
|
||||
std::cout << "Sixth test passes.\n";
|
||||
// Test for absolute value of our complex number
|
||||
assert(((void)"(1 + 1i) has absolute value sqrt(2) but the program differs "
|
||||
"from the std::complex result. \n",
|
||||
(num1.abs() == std::abs(cnum1))));
|
||||
std::cout << "Seventh test passes.\n";
|
||||
}
|
||||
|
||||
/**
|
||||
* Main function
|
||||
*/
|
||||
int main() {
|
||||
tests();
|
||||
return 0;
|
||||
}
|
||||
@ -1,8 +1,9 @@
|
||||
/**
|
||||
* @file
|
||||
* @brief Compute double factorial: \f$n!!\f$
|
||||
* @brief Compute [double
|
||||
* factorial](https://en.wikipedia.org/wiki/Double_factorial): \f$n!!\f$
|
||||
*
|
||||
* Double factorial of a non-negative integer n, is defined as the product of
|
||||
* Double factorial of a non-negative integer `n`, is defined as the product of
|
||||
* all the integers from 1 to n that have the same parity (odd or even) as n.
|
||||
* <br/>It is also called as semifactorial of a number and is denoted by
|
||||
* \f$n!!\f$
|
||||
@ -32,10 +33,38 @@ uint64_t double_factorial_recursive(uint64_t n) {
|
||||
return n * double_factorial_recursive(n - 2);
|
||||
}
|
||||
|
||||
/// main function
|
||||
int main() {
|
||||
uint64_t n;
|
||||
std::cin >> n;
|
||||
assert(n >= 0);
|
||||
std::cout << double_factorial_iterative(n);
|
||||
/** Wrapper to run tests using both recursive and iterative implementations.
|
||||
* The checks are only valid in debug builds due to the use of `assert()`
|
||||
* statements.
|
||||
* \param [in] n number to check double factorial for
|
||||
* \param [in] expected expected result
|
||||
*/
|
||||
void test(uint64_t n, uint64_t expected) {
|
||||
assert(double_factorial_iterative(n) == expected);
|
||||
assert(double_factorial_recursive(n) == expected);
|
||||
}
|
||||
|
||||
/**
|
||||
* Test implementations
|
||||
*/
|
||||
void tests() {
|
||||
std::cout << "Test 1:\t n=5\t...";
|
||||
test(5, 15);
|
||||
std::cout << "passed\n";
|
||||
|
||||
std::cout << "Test 2:\t n=15\t...";
|
||||
test(15, 2027025);
|
||||
std::cout << "passed\n";
|
||||
|
||||
std::cout << "Test 3:\t n=0\t...";
|
||||
test(0, 1);
|
||||
std::cout << "passed\n";
|
||||
}
|
||||
|
||||
/**
|
||||
* Main function
|
||||
*/
|
||||
int main() {
|
||||
tests();
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -19,28 +19,32 @@
|
||||
#include <cstdio>
|
||||
#include <iostream>
|
||||
|
||||
/** maximum number that can be computed - The result after 93 cannot be stored
|
||||
* in a `uint64_t` data type. */
|
||||
const uint64_t MAX = 93;
|
||||
/**
|
||||
* maximum number that can be computed - The result after 93 cannot be stored
|
||||
* in a `uint64_t` data type.
|
||||
*/
|
||||
|
||||
/** Array of computed fibonacci numbers */
|
||||
uint64_t f[MAX] = {0};
|
||||
#define MAX 93
|
||||
|
||||
/** Algorithm */
|
||||
uint64_t fib(uint64_t n) {
|
||||
if (n == 0)
|
||||
static uint64_t f1 = 1,
|
||||
f2 = 1; // using static keyword will retain the values of
|
||||
// f1 and f2 for the next function call.
|
||||
|
||||
if (n <= 2)
|
||||
return f2;
|
||||
if (n >= 93) {
|
||||
std::cerr
|
||||
<< "Cannot compute for n>93 due to limit of 64-bit integers\n";
|
||||
return 0;
|
||||
if (n == 1 || n == 2)
|
||||
return (f[n] = 1);
|
||||
}
|
||||
|
||||
if (f[n])
|
||||
return f[n];
|
||||
uint64_t temp = f2; // we do not need temp to be static
|
||||
f2 += f1;
|
||||
f1 = temp;
|
||||
|
||||
uint64_t k = (n % 2 != 0) ? (n + 1) / 2 : n / 2;
|
||||
|
||||
f[n] = (n % 2 != 0) ? (fib(k) * fib(k) + fib(k - 1) * fib(k - 1))
|
||||
: (2 * fib(k - 1) + fib(k)) * fib(k);
|
||||
return f[n];
|
||||
return f2;
|
||||
}
|
||||
|
||||
/** Main function */
|
||||
|
||||
@ -35,7 +35,7 @@ class stats_computer1 {
|
||||
n++;
|
||||
T tmp = x - K;
|
||||
Ex += tmp;
|
||||
Ex2 += tmp * tmp;
|
||||
Ex2 += static_cast<double>(tmp) * tmp;
|
||||
}
|
||||
|
||||
/** return sample mean computed till last sample */
|
||||
|
||||
@ -4,76 +4,18 @@
|
||||
* square matrix
|
||||
* \author [Krishna Vedala](https://github.com/kvedala)
|
||||
*/
|
||||
#include <cassert>
|
||||
#include <ctime>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
/** Perform LU decomposition on matrix
|
||||
* \param[in] A matrix to decompose
|
||||
* \param[out] L output L matrix
|
||||
* \param[out] U output U matrix
|
||||
* \returns 0 if no errors
|
||||
* \returns negative if error occurred
|
||||
*/
|
||||
int lu_decomposition(const std::vector<std::vector<double>> &A,
|
||||
std::vector<std::vector<double>> *L,
|
||||
std::vector<std::vector<double>> *U) {
|
||||
int row, col, j;
|
||||
int mat_size = A.size();
|
||||
|
||||
if (mat_size != A[0].size()) {
|
||||
// check matrix is a square matrix
|
||||
std::cerr << "Not a square matrix!\n";
|
||||
return -1;
|
||||
}
|
||||
|
||||
// regularize each row
|
||||
for (row = 0; row < mat_size; row++) {
|
||||
// Upper triangular matrix
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
for (col = row; col < mat_size; col++) {
|
||||
// Summation of L[i,j] * U[j,k]
|
||||
double lu_sum = 0.;
|
||||
for (j = 0; j < row; j++) lu_sum += L[0][row][j] * U[0][j][col];
|
||||
|
||||
// Evaluate U[i,k]
|
||||
U[0][row][col] = A[row][col] - lu_sum;
|
||||
}
|
||||
|
||||
// Lower triangular matrix
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
for (col = row; col < mat_size; col++) {
|
||||
if (row == col) {
|
||||
L[0][row][col] = 1.;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Summation of L[i,j] * U[j,k]
|
||||
double lu_sum = 0.;
|
||||
for (j = 0; j < row; j++) lu_sum += L[0][col][j] * U[0][j][row];
|
||||
|
||||
// Evaluate U[i,k]
|
||||
L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
#include "./lu_decomposition.h"
|
||||
|
||||
/**
|
||||
* operator to print a matrix
|
||||
*/
|
||||
template <typename T>
|
||||
std::ostream &operator<<(std::ostream &out,
|
||||
std::vector<std::vector<T>> const &v) {
|
||||
std::ostream &operator<<(std::ostream &out, matrix<T> const &v) {
|
||||
const int width = 10;
|
||||
const char separator = ' ';
|
||||
|
||||
@ -87,26 +29,21 @@ std::ostream &operator<<(std::ostream &out,
|
||||
return out;
|
||||
}
|
||||
|
||||
/** Main function */
|
||||
int main(int argc, char **argv) {
|
||||
/**
|
||||
* Test LU decomposition
|
||||
* \todo better ways to self-check a matrix output?
|
||||
*/
|
||||
void test1() {
|
||||
int mat_size = 3; // default matrix size
|
||||
const int range = 50;
|
||||
const int range2 = range >> 1;
|
||||
|
||||
if (argc == 2)
|
||||
mat_size = atoi(argv[1]);
|
||||
|
||||
std::srand(std::time(NULL)); // random number initializer
|
||||
|
||||
/* Create a square matrix with random values */
|
||||
std::vector<std::vector<double>> A(mat_size);
|
||||
std::vector<std::vector<double>> L(mat_size); // output
|
||||
std::vector<std::vector<double>> U(mat_size); // output
|
||||
matrix<double> A(mat_size, std::valarray<double>(mat_size));
|
||||
matrix<double> L(mat_size, std::valarray<double>(mat_size)); // output
|
||||
matrix<double> U(mat_size, std::valarray<double>(mat_size)); // output
|
||||
for (int i = 0; i < mat_size; i++) {
|
||||
// calloc so that all valeus are '0' by default
|
||||
A[i] = std::vector<double>(mat_size);
|
||||
L[i] = std::vector<double>(mat_size);
|
||||
U[i] = std::vector<double>(mat_size);
|
||||
for (int j = 0; j < mat_size; j++)
|
||||
/* create random values in the limits [-range2, range-1] */
|
||||
A[i][j] = static_cast<double>(std::rand() % range - range2);
|
||||
@ -121,6 +58,33 @@ int main(int argc, char **argv) {
|
||||
std::cout << "A = \n" << A << "\n";
|
||||
std::cout << "L = \n" << L << "\n";
|
||||
std::cout << "U = \n" << U << "\n";
|
||||
}
|
||||
|
||||
/**
|
||||
* Test determinant computation using LU decomposition
|
||||
*/
|
||||
void test2() {
|
||||
std::cout << "Determinant test 1...";
|
||||
matrix<int> A1({{1, 2, 3}, {4, 9, 6}, {7, 8, 9}});
|
||||
assert(determinant_lu(A1) == -48);
|
||||
std::cout << "passed\n";
|
||||
|
||||
std::cout << "Determinant test 2...";
|
||||
matrix<int> A2({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}});
|
||||
assert(determinant_lu(A2) == 0);
|
||||
std::cout << "passed\n";
|
||||
|
||||
std::cout << "Determinant test 3...";
|
||||
matrix<float> A3({{1.2, 2.3, 3.4}, {4.5, 5.6, 6.7}, {7.8, 8.9, 9.0}});
|
||||
assert(determinant_lu(A3) == 3.63);
|
||||
std::cout << "passed\n";
|
||||
}
|
||||
|
||||
/** Main function */
|
||||
int main(int argc, char **argv) {
|
||||
std::srand(std::time(NULL)); // random number initializer
|
||||
|
||||
test1();
|
||||
test2();
|
||||
return 0;
|
||||
}
|
||||
|
||||
102
numerical_methods/lu_decomposition.h
Normal file
102
numerical_methods/lu_decomposition.h
Normal file
@ -0,0 +1,102 @@
|
||||
/**
|
||||
* @file lu_decomposition.h
|
||||
* @author [Krishna Vedala](https://github.com/kvedala)
|
||||
* @brief Functions associated with [LU
|
||||
* Decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
|
||||
* of a square matrix.
|
||||
*/
|
||||
#pragma once
|
||||
|
||||
#include <iostream>
|
||||
#include <valarray>
|
||||
#include <vector>
|
||||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
|
||||
/** Define matrix type as a `std::vector` of `std::valarray` */
|
||||
template <typename T>
|
||||
using matrix = std::vector<std::valarray<T>>;
|
||||
|
||||
/** Perform LU decomposition on matrix
|
||||
* \param[in] A matrix to decompose
|
||||
* \param[out] L output L matrix
|
||||
* \param[out] U output U matrix
|
||||
* \returns 0 if no errors
|
||||
* \returns negative if error occurred
|
||||
*/
|
||||
template <typename T>
|
||||
int lu_decomposition(const matrix<T> &A, matrix<double> *L, matrix<double> *U) {
|
||||
int row, col, j;
|
||||
int mat_size = A.size();
|
||||
|
||||
if (mat_size != A[0].size()) {
|
||||
// check matrix is a square matrix
|
||||
std::cerr << "Not a square matrix!\n";
|
||||
return -1;
|
||||
}
|
||||
|
||||
// regularize each row
|
||||
for (row = 0; row < mat_size; row++) {
|
||||
// Upper triangular matrix
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
for (col = row; col < mat_size; col++) {
|
||||
// Summation of L[i,j] * U[j,k]
|
||||
double lu_sum = 0.;
|
||||
for (j = 0; j < row; j++) {
|
||||
lu_sum += L[0][row][j] * U[0][j][col];
|
||||
}
|
||||
|
||||
// Evaluate U[i,k]
|
||||
U[0][row][col] = A[row][col] - lu_sum;
|
||||
}
|
||||
|
||||
// Lower triangular matrix
|
||||
#ifdef _OPENMP
|
||||
#pragma omp for
|
||||
#endif
|
||||
for (col = row; col < mat_size; col++) {
|
||||
if (row == col) {
|
||||
L[0][row][col] = 1.;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Summation of L[i,j] * U[j,k]
|
||||
double lu_sum = 0.;
|
||||
for (j = 0; j < row; j++) {
|
||||
lu_sum += L[0][col][j] * U[0][j][row];
|
||||
}
|
||||
|
||||
// Evaluate U[i,k]
|
||||
L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute determinant of an NxN square matrix using LU decomposition.
|
||||
* Using LU decomposition, the determinant is given by the product of diagonal
|
||||
* elements of matrices L and U.
|
||||
*
|
||||
* @tparam T datatype of input matrix - int, unsigned int, double, etc
|
||||
* @param A input square matrix
|
||||
* @return determinant of matrix A
|
||||
*/
|
||||
template <typename T>
|
||||
double determinant_lu(const matrix<T> &A) {
|
||||
matrix<double> L(A.size(), std::valarray<double>(A.size()));
|
||||
matrix<double> U(A.size(), std::valarray<double>(A.size()));
|
||||
|
||||
if (lu_decomposition(A, &L, &U) < 0)
|
||||
return 0;
|
||||
|
||||
double result = 1.f;
|
||||
for (size_t i = 0; i < A.size(); i++) {
|
||||
result *= L[i][i] * U[i][i];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
@ -39,18 +39,15 @@ using std::vector;
|
||||
#define pb push_back
|
||||
#define MOD 1000000007
|
||||
|
||||
/** returns absolute value */
|
||||
inline ll ab(ll x) { return x > 0LL ? x : -x; }
|
||||
|
||||
/** global variable k
|
||||
/** global variable mat_size
|
||||
* @todo @stepfencurryxiao add documetnation
|
||||
*/
|
||||
ll k;
|
||||
ll mat_size;
|
||||
|
||||
/** global vector variables
|
||||
/** global vector variables used in the ::ans function.
|
||||
* @todo @stepfencurryxiao add documetnation
|
||||
*/
|
||||
vector<ll> a, b, c;
|
||||
vector<ll> fib_b, fib_c;
|
||||
|
||||
/** To multiply 2 matrices
|
||||
* \param [in] A matrix 1 of size (m\f$\times\f$n)
|
||||
@ -59,10 +56,10 @@ vector<ll> a, b, c;
|
||||
*/
|
||||
vector<vector<ll>> multiply(const vector<vector<ll>> &A,
|
||||
const vector<vector<ll>> &B) {
|
||||
vector<vector<ll>> C(k + 1, vector<ll>(k + 1));
|
||||
for (ll i = 1; i <= k; i++) {
|
||||
for (ll j = 1; j <= k; j++) {
|
||||
for (ll z = 1; z <= k; z++) {
|
||||
vector<vector<ll>> C(mat_size + 1, vector<ll>(mat_size + 1));
|
||||
for (ll i = 1; i <= mat_size; i++) {
|
||||
for (ll j = 1; j <= mat_size; j++) {
|
||||
for (ll z = 1; z <= mat_size; z++) {
|
||||
C[i][j] = (C[i][j] + (A[i][z] * B[z][j]) % MOD) % MOD;
|
||||
}
|
||||
}
|
||||
@ -94,24 +91,24 @@ vector<vector<ll>> power(const vector<vector<ll>> &A, ll p) {
|
||||
ll ans(ll n) {
|
||||
if (n == 0)
|
||||
return 0;
|
||||
if (n <= k)
|
||||
return b[n - 1];
|
||||
if (n <= mat_size)
|
||||
return fib_b[n - 1];
|
||||
// F1
|
||||
vector<ll> F1(k + 1);
|
||||
for (ll i = 1; i <= k; i++) F1[i] = b[i - 1];
|
||||
vector<ll> F1(mat_size + 1);
|
||||
for (ll i = 1; i <= mat_size; i++) F1[i] = fib_b[i - 1];
|
||||
|
||||
// Transpose matrix
|
||||
vector<vector<ll>> T(k + 1, vector<ll>(k + 1));
|
||||
for (ll i = 1; i <= k; i++) {
|
||||
for (ll j = 1; j <= k; j++) {
|
||||
if (i < k) {
|
||||
vector<vector<ll>> T(mat_size + 1, vector<ll>(mat_size + 1));
|
||||
for (ll i = 1; i <= mat_size; i++) {
|
||||
for (ll j = 1; j <= mat_size; j++) {
|
||||
if (i < mat_size) {
|
||||
if (j == i + 1)
|
||||
T[i][j] = 1;
|
||||
else
|
||||
T[i][j] = 0;
|
||||
continue;
|
||||
}
|
||||
T[i][j] = c[k - j];
|
||||
T[i][j] = fib_c[mat_size - j];
|
||||
}
|
||||
}
|
||||
// T^n-1
|
||||
@ -119,7 +116,7 @@ ll ans(ll n) {
|
||||
|
||||
// T*F1
|
||||
ll res = 0;
|
||||
for (ll i = 1; i <= k; i++) {
|
||||
for (ll i = 1; i <= mat_size; i++) {
|
||||
res = (res + (T[1][i] * F1[i]) % MOD) % MOD;
|
||||
}
|
||||
return res;
|
||||
@ -133,19 +130,19 @@ int main() {
|
||||
cin >> t;
|
||||
ll i, j, x;
|
||||
while (t--) {
|
||||
cin >> k;
|
||||
for (i = 0; i < k; i++) {
|
||||
cin >> mat_size;
|
||||
for (i = 0; i < mat_size; i++) {
|
||||
cin >> x;
|
||||
b.pb(x);
|
||||
fib_b.pb(x);
|
||||
}
|
||||
for (i = 0; i < k; i++) {
|
||||
for (i = 0; i < mat_size; i++) {
|
||||
cin >> x;
|
||||
c.pb(x);
|
||||
fib_c.pb(x);
|
||||
}
|
||||
cin >> x;
|
||||
cout << ans(x) << endl;
|
||||
b.clear();
|
||||
c.clear();
|
||||
fib_b.clear();
|
||||
fib_c.clear();
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -20,13 +20,13 @@
|
||||
char stack[MAX];
|
||||
|
||||
//! pointer to track stack index
|
||||
int top = -1;
|
||||
int stack_idx = -1;
|
||||
|
||||
//! push byte to stack variable
|
||||
void push(char ch) { stack[++top] = ch; }
|
||||
void push(char ch) { stack[++stack_idx] = ch; }
|
||||
|
||||
//! pop a byte out of stack variable
|
||||
char pop() { return stack[top--]; }
|
||||
char pop() { return stack[stack_idx--]; }
|
||||
|
||||
//! @}-------------- end stack -----------
|
||||
|
||||
@ -56,7 +56,7 @@ int main() {
|
||||
while (valid == 1 && i < exp.length()) {
|
||||
if (exp[i] == '(' || exp[i] == '{' || exp[i] == '[' || exp[i] == '<') {
|
||||
push(exp[i]);
|
||||
} else if (top >= 0 && stack[top] == opening(exp[i])) {
|
||||
} else if (stack_idx >= 0 && stack[stack_idx] == opening(exp[i])) {
|
||||
pop();
|
||||
} else {
|
||||
valid = 0;
|
||||
@ -65,7 +65,7 @@ int main() {
|
||||
}
|
||||
|
||||
// makes sure the stack is empty after processsing (above)
|
||||
if (valid == 1 && top == -1) {
|
||||
if (valid == 1 && stack_idx == -1) {
|
||||
std::cout << "\nCorrect Expression";
|
||||
} else {
|
||||
std::cout << "\nWrong Expression";
|
||||
|
||||
@ -14,7 +14,7 @@ void beadSort(int *a, int len) {
|
||||
|
||||
// allocating memory
|
||||
unsigned char *beads = new unsigned char[max * len];
|
||||
memset(beads, 0, max * len);
|
||||
memset(beads, 0, static_cast<size_t>(max) * len);
|
||||
|
||||
// mark the beads
|
||||
for (int i = 0; i < len; i++)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user