- 🧑💻 I'm a master's student.
- 🌱 I’m currently learning lattice theory!
- 📫 How to reach me: Twitter - @satoshin_NCD
C | C++ | C# | Cython | HSP |
---|---|---|---|---|
maxima | processing | Python | Risa/Asir | SageMath |
wenyan | ||||
DXライブラリ | Eigen | IPython | matplotlib | Numpy |
---|---|---|---|---|
pandas | plotly | scipy | sympy | |
clang | GCC |
---|---|
CSS | HTML | Markdown | LaTeX |
---|---|---|---|
Android | Ubuntu | Windows |
---|---|---|
Unity |
---|
Git | GitHub | Sourcetree |
---|---|---|
sublime text | VScode |
---|---|
Visual Studio |
---|
AviUtl | ibis | Inkscape | MuseScore | NEUTRINO |
---|---|---|---|---|
SoundEngine | UTAU | VOICEVOX | COEIROINK | Domino |
C
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* inner product */
double dot_dbl_dbl(double *x, double *y, const int n){
double s = 0.0;
for(int i = 0; i < n; ++i) s += x[i] * y[i];
return s;
}
double dot_int_dbl(int *x, double *y, const int n){
double s = 0.0;
for(int i = 0; i < n; ++i) s += y[i] * x[i];
return s;
}
/* Gram-Schmidt's method */
void GSO(int **b, double *B, double **mu, const int n, const int m){
int i, j, k;
double t, s, **GSOb;
GSOb = (double **)malloc(n * sizeof(double *));
for(i = 0; i < n; ++i) GSOb[i] = (double *)malloc(m * sizeof(double));
for(i = 0; i < n; ++i){
mu[i][i] = 1.0;
for(j = 0; j < m; ++j) GSOb[i][j] = b[i][j];
for(j = 0; j < i; ++j){
mu[i][j] = dot_int_dbl(b[i], GSOb[j], m) / dot_dbl_dbl(GSOb[j], GSOb[j], m);
for(k = 0; k < m; ++k) GSOb[i][k] -= mu[i][j] * GSOb[j][k];
}
B[i] = dot_dbl_dbl(GSOb[i], GSOb[i], m);
}
}
/* size-reduction */
void SizeReduce(int **b, double **mu, const int i, const int j, const int m){
int k;
if(mu[i][j] > 0.5 || mu[i][j] < -0.5){
const int q = round(mu[i][j]);
for(k = 0; k < m; ++k) b[i][k] -= q * b[j][k];
for(k = 0; k <= j; ++k) mu[i][k] -= mu[j][k] * q;
}
}
/* LLL-reduction */
void LLLReduce(int **b, const double d, const int n, const int m){
int j, i, h;
double **mu, *B, nu, BB, C, t;
mu = (double **)malloc(n * sizeof(double *));
B = (double *)malloc(n * sizeof(double));
for(i = 0; i < n; ++i) mu[i] = (double *)malloc(n * sizeof(double));
GSO(b, B, mu, n, m);
int tmp;
for(int k = 1; k < n;){
h = k - 1;
for(j = h; j > -1; --j) SizeReduce(b, mu, k, j, m);
if(k > 0 && B[k] < (d - mu[k][h] * mu[k][h]) * B[h]){
for(i = 0; i < m; ++i){tmp = b[h][i]; b[h][i] = b[k][i]; b[k][i] = tmp;}
nu = mu[k][k - 1]; BB = B[k] + nu * nu * B[k - 1]; C = 1.0 / BB;
mu[k][k - 1] = nu * B[k - 1] * C; B[k] *= B[k - 1] * C; B[k - 1] = BB;
for(i = 0; i <= k - 2; ++i){
t = mu[k - 1][i]; mu[k - 1][i] = mu[k][i]; mu[k][i] = t;
}
for(i = k + 1; i < n; ++i){
t = mu[i][k]; mu[i][k] = mu[i][k - 1] - nu * t;
mu[i][k - 1] = t + mu[k][k - 1] * mu[i][k];
}
k = h;
}else ++k;
}
}
C++
#include <iostream>
#include <vector>
#include <tuple>
/* inner product */
double dot(const std::vector<int> x, const std::vector<double> y){
double z = 0.0;
const int n = x.size();
for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
return z;
}
double dot(const std::vector<double> x, const std::vector<double> y){
double z = 0.0;
const int n = x.size();
for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
return z;
}
double dot(const std::vector<int> x, const std::vector<int> y){
double z = 0.0;
const int n = x.size();
for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
return z;
}
/* Gram-Schmidt's method */
std::tuple<std::vector<double>, std::vector<std::vector<double>>> Gram_Schmidt_squared(const std::vector<std::vector<int>> b){
const int n = b.size(), m = b.at(0).size(); int i, j, k;
std::vector<double> B(n);
std::vector<std::vector<double>> GSOb(n, std::vector<double>(m)), mu(n, std::vector<double>(n));
for(i = 0; i < n; ++i){
mu.at(i).at(i)= 1.0;
for(j = 0; j < m; ++j) GSOb.at(i).at(j) = b.at(i).at(j);
for(j = 0; j < i; ++j){
mu.at(i).at(j) = dot(b.at(i), GSOb.at(j)) / dot(GSOb.at(j), GSOb.at(j));
for(k = 0; k < m; ++k) GSOb.at(i).at(k) -= mu.at(i).at(j) * GSOb.at(j).at(k);
}
B.at(i) = dot(GSOb.at(i), GSOb.at(i));
}
return std::forward_as_tuple(B, mu);
}
/* size-reduction */
void SizeReduce(std::vector<std::vector<int>>& b, std::vector<std::vector<double>>& mu, const int i, const int j){
int q;
const int m = b.at(0).size();
if(mu.at(i).at(j) > 0.5 || mu.at(i).at(j) < -0.5){
q = round(mu.at(i).at(j));
for(int k = 0; k < m; ++k) b.at(i).at(k) -= q * b.at(j).at(k);
for(int k = 0; k <= j; ++k) mu.at(i).at(k) -= mu.at(j).at(k) * q;
}
}
/* LLL-reduction */
void LLLReduce(std::vector<std::vector<int>>& b, const float d = 0.99){
const int n = b.size(), m = b.at(0).size(); int j, i, h;
double t, nu, BB, C;
std::vector<std::vector<double>> mu;
std::vector<double> B; std::tie(B, mu) = Gram_Schmidt_squared(b);
int tmp;
for(int k = 1; k < n;){
h = k - 1;
for(j = h; j > -1; --j) SizeReduce(b, mu, k, j);
//Checks if the lattice basis matrix b satisfies Lovasz condition.
if(k > 0 && B.at(k) < (d - mu.at(k).at(h) * mu.at(k).at(h)) * B.at(h)){
for(i = 0; i < m; ++i){tmp = b.at(h).at(i); b.at(h).at(i) = b.at(k).at(i); b.at(k).at(i) = tmp;}
nu = mu.at(k).at(h); BB = B.at(k) + nu * nu * B.at(h); C = 1.0 / BB;
mu.at(k).at(h) = nu * B.at(h) * C; B[k] *= B.at(h) * C; B.at(h) = BB;
for(i = 0; i <= k - 2; ++i){
t = mu.at(h).at(i); mu.at(h).at(i) = mu.at(k).at(i); mu.at(k).at(i) = t;
}
for(i = k + 1; i < n; ++i){
t = mu.at(i).at(k); mu.at(i).at(k) = mu.at(i).at(h) - nu * t;
mu.at(i).at(h) = t + mu.at(k).at(h) * mu.at(i).at(k);
}
--k;
}else ++k;
}
}
C++ with Eigen library
#include <iostream>
#include <eigen3/Eigen/Dense>
/* Gram-Schmidt's method */
void GSO(const Eigen::MatrixXi b, Eigen::VectorXd& B, Eigen::MatrixXd& mu, const int n, const int m){
int j;
Eigen::MatrixXd GSOb(n, m);
for(int i = 0; i < n; ++i){
mu.coeffRef(i, i) = 1.0;
GSOb.row(i) = b.row(i).cast<double>();
for(j = 0; j < i; ++j){
mu.coeffRef(i, j) = b.row(i).cast<double>().dot(GSOb.row(j)) / GSOb.row(j).dot(GSOb.row(j));
GSOb.row(i) -= mu.coeff(i, j) * GSOb.row(j);
}
B.coeffRef(i) = GSOb.row(i).dot(GSOb.row(i));
}
}
/* size-reduction */
void SizeReduce(Eigen::MatrixXi& b, Eigen::MatrixXd& mu, const int i, const int j, const int m){
if(mu.coeff(i, j) > 0.5 || mu.coeff(i, j) < -0.5){
const int q = round(mu.coeff(i, j));
b.row(i) -= q * b.row(j);
mu.row(i).head(j + 1) -= (double)q * mu.row(j).head(j + 1);
}
}
/* LLL-reduction */
void LLLReduce(Eigen::MatrixXi& b, const long double d, const int n, const int m){
double nu, BB, C, t;
Eigen::VectorXd B(n), logB(n);
Eigen::MatrixXd mu(n, n);
GSO(b, B, mu, n, m);
for(int k = 1, j, i, k1; k < n;){
k1 = k - 1;
for(j = k1; j > -1; --j) SizeReduce(b, mu, k, j, m);
if(k > 0 && B.coeff(k) < (d - mu.coeff(k, k1) * mu.coeff(k, k1)) * B.coeff(k1)){
b.row(k).swap(b.row(k1));
nu = mu.coeff(k, k1); BB = B.coeff(k) + nu * nu * B.coeff(k1); C = 1.0 / BB;
mu.coeffRef(k, k1) = nu * B.coeff(k1) * C;
B.coeffRef(k) *= B.coeff(k1) * C; B.coeffRef(k1) = BB;
mu.row(k1).head(k - 1).swap(mu.row(k).head(k - 1));
for(i = k + 1; i < n; ++i){
t = mu.coeff(i, k); mu.coeffRef(i, k) = mu.coeff(i, k1) - nu * t;
mu.coeffRef(i, k1) = t + mu.coeff(k, k1) * mu.coeff(i, k);
}
k = k1;
}else ++k;
}
}