ADD file via upload
This commit is contained in:
parent
e04c13cf0d
commit
5ae15858dd
|
|
@ -0,0 +1,304 @@
|
|||
#ifndef MATRIX_H
|
||||
#define MATRIX_H
|
||||
#include "conf.h"
|
||||
#include "polynomial.h"
|
||||
#include "polynomial_store.h"
|
||||
|
||||
struct Matrix {
|
||||
Polynomial* R[reductors];
|
||||
Sparse_Polynomial* SR[reductors];
|
||||
int R_leader[reductors];
|
||||
bool R_sparse[reductors]; //0稠密1稀疏
|
||||
Polynomial* E[MAX_POLYS];
|
||||
int end_elim[MAX_POLYS];
|
||||
int reductor_rows;
|
||||
int elim_poly_rows;
|
||||
int reductor_block; //当前消元块首块：算法以64位为单位进行消元
|
||||
int reductors_not_null;
|
||||
int reduce_record[MAX_POLYS][reductors];
|
||||
int reduce_record_len[MAX_POLYS];
|
||||
int reduce_record_sparse[MAX_POLYS][reductors];
|
||||
int reduce_record_sparse_len[MAX_POLYS];
|
||||
|
||||
inline Matrix()
|
||||
{
|
||||
clear();
|
||||
}
|
||||
inline int add_reductor(Store& ss, Temp_Polynomial_Store& ps_temp, int begin, int end, int mode = 1) //添加消元子
|
||||
{
|
||||
ps_temp.clear();
|
||||
memset(R, 0, sizeof(R));
|
||||
memset(SR, 0, sizeof(SR));
|
||||
memset(R_sparse, 0, sizeof(R_sparse));
|
||||
reductors_not_null = 0;
|
||||
//int now_row = 0; //添加消元子到now_row行
|
||||
if (mode == 1)
|
||||
{
|
||||
parallel_for(begin, end, [&](int i) {
|
||||
int now_row = i - begin;
|
||||
if (ss[i].type == 0)
|
||||
{
|
||||
R[now_row] = ss[i].poly;
|
||||
SR[now_row] = NULL;
|
||||
}
|
||||
else if (ss[i].type == 1)
|
||||
{
|
||||
int set_type = ps_temp.set(ss[i].mul_poly, ss[i].mul_mono, now_row);
|
||||
if (set_type == 1)
|
||||
{
|
||||
R[now_row] = ps_temp.get_dense(now_row);
|
||||
SR[now_row] = NULL;
|
||||
}
|
||||
else
|
||||
{
|
||||
R[now_row] = NULL;
|
||||
SR[now_row] = ps_temp.get_sparse(now_row);
|
||||
R_sparse[now_row] = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
R[now_row] = NULL;
|
||||
SR[now_row] = NULL;
|
||||
}
|
||||
R_leader[now_row] = i;
|
||||
});
|
||||
}
|
||||
if (mode == 2) //得到既约Groebner基，默认约化得到最少的多项式
|
||||
{
|
||||
parallel_for(begin, end, [&](int i) {
|
||||
int now_row = i - begin;
|
||||
if (ss[i].type == 1)
|
||||
{
|
||||
int set_type = ps_temp.set(ss[i].mul_poly, ss[i].mul_mono, now_row);
|
||||
if (set_type == 1)
|
||||
{
|
||||
R[now_row] = ps_temp.get_dense(now_row);
|
||||
SR[now_row] = NULL;
|
||||
}
|
||||
else
|
||||
{
|
||||
R[now_row] = NULL;
|
||||
SR[now_row] = ps_temp.get_sparse(now_row);
|
||||
R_sparse[now_row] = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
R[now_row] = NULL;
|
||||
SR[now_row] = NULL;
|
||||
}
|
||||
R_leader[now_row] = i;
|
||||
});
|
||||
}
|
||||
reductor_rows = end - begin;
|
||||
reductor_block = begin / 64;
|
||||
for (int i = 0; i < reductor_rows; i++)
|
||||
{
|
||||
if (R[i] != NULL || SR[i] != NULL)
|
||||
{
|
||||
reductors_not_null++;
|
||||
}
|
||||
}
|
||||
return end - begin;
|
||||
}
|
||||
inline int add_elim_poly(Polynomial_Store& ps, int begin = 0, int end = M, int degree = D) //添加被消元多项式
|
||||
{
|
||||
memset(E, 0, sizeof(E));
|
||||
int row = 0;
|
||||
for (int i = begin; i < end; i++)
|
||||
{
|
||||
if (ps[i].degree > degree)
|
||||
continue;
|
||||
if (!ps[i].iszero()) //不导入零多项式
|
||||
{
|
||||
E[row] = &(ps[i]);
|
||||
row++;
|
||||
}
|
||||
}
|
||||
elim_poly_rows = row;
|
||||
return row;
|
||||
}
|
||||
inline void clear()
|
||||
{
|
||||
memset(end_elim, 0, sizeof(end_elim));
|
||||
memset(R, 0, sizeof(R));
|
||||
memset(R_leader, 0, sizeof(R_leader));
|
||||
memset(E, 0, sizeof(E));
|
||||
memset(SR, 0, sizeof(SR));
|
||||
memset(R_sparse, 0, sizeof(R_sparse));
|
||||
reductor_rows = 0;
|
||||
elim_poly_rows = 0;
|
||||
reductor_block = 0;
|
||||
}
|
||||
inline Polynomial*& get_elim_poly(int n)
|
||||
{
|
||||
return E[n];
|
||||
}
|
||||
inline void Gauss_elimination()
|
||||
{
|
||||
/*
|
||||
memset(reduce_record, 0, sizeof(reduce_record));
|
||||
memset(reduce_record_len, 0, sizeof(reduce_record_len));
|
||||
memset(reduce_record_sparse, 0, sizeof(reduce_record_sparse));
|
||||
memset(reduce_record_sparse_len, 0, sizeof(reduce_record_sparse_len));
|
||||
*/
|
||||
//只将不初始化会引起错误的内存初始化，没有初始化的内存并不会被直接读取
|
||||
memset(reduce_record_len, 0, elim_poly_rows * sizeof(int));
|
||||
memset(reduce_record_sparse_len, 0, elim_poly_rows * sizeof(int));
|
||||
if (reductors_not_null == reductors) //消元子满编：本部分不会产生新首项
|
||||
{
|
||||
parallel_for(0, elim_poly_rows, [&](int i) {
|
||||
uint64* temp_block = E[i]->mono.data + reductor_block;
|
||||
bool zero_block = true;
|
||||
if (end_elim[i])
|
||||
return; //跳出本次循环，相当于for循环中的continue
|
||||
for (int j = 0; j < reductors / 64; j++)
|
||||
{
|
||||
if (temp_block[j] != 0)
|
||||
{
|
||||
zero_block = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (zero_block)
|
||||
return;
|
||||
for (int j = reductor_rows - 1; j >= 0; j--)
|
||||
{
|
||||
if (temp_block[j / 64] & (1ull << (j % 64)))
|
||||
{
|
||||
if (R_sparse[j]) //该位置一定存在一个稀疏多项式
|
||||
{
|
||||
xor_sparse(*E[i], *SR[j]);
|
||||
}
|
||||
else
|
||||
{
|
||||
E[i]->xor_only(*R[j], 0, reductor_block + reductors / 64);
|
||||
}
|
||||
}
|
||||
}
|
||||
if (E[i]->mono.data_len > reductor_block)
|
||||
E[i]->mono.data_len = reductor_block;
|
||||
});
|
||||
}
|
||||
else
|
||||
{
|
||||
for (int i = 0; i < elim_poly_rows; i++)
|
||||
{
|
||||
int& temp_count = reduce_record_len[i];
|
||||
int& sparse_count = reduce_record_sparse_len[i];
|
||||
uint64 temp_block[reductors / 64] = { 0 };
|
||||
bool zero_block = true;
|
||||
for (int j = 0; j < reductors / 64; j++)
|
||||
temp_block[j] = E[i]->mono.data[reductor_block + j];
|
||||
if (end_elim[i])
|
||||
continue;
|
||||
for (int j = 0; j < reductors / 64; j++)
|
||||
{
|
||||
if (temp_block[j] != 0)
|
||||
{
|
||||
zero_block = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (int j = reductor_rows - 1; j >= 0; j--)
|
||||
{
|
||||
if (temp_block[j / 64] & (1ull << (j % 64)))
|
||||
{
|
||||
if (R_sparse[j]) //该位置一定存在一个稀疏多项式
|
||||
{
|
||||
for (int k = 0; k < reductors / 64; k++)
|
||||
{
|
||||
temp_block[k] ^= SR[j]->leader_block[k];
|
||||
}
|
||||
reduce_record_sparse[i][sparse_count] = j;
|
||||
sparse_count++;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (R[j] == NULL)
|
||||
{
|
||||
//该行会消元得出新首项，将其直接计算出
|
||||
for (int k = 0; k < reduce_record_len[i]; k++)
|
||||
{
|
||||
E[i]->xor_only(*R[reduce_record[i][k]], 0, reductor_block + reductors / 64);
|
||||
}
|
||||
for (int k = 0; k < reduce_record_sparse_len[i]; k++)
|
||||
{
|
||||
xor_sparse(*E[i], *SR[reduce_record_sparse[i][k]]);
|
||||
}
|
||||
E[i]->set_leader_degree(reductor_block * 64 + reductors - 1);
|
||||
R[j] = E[i];
|
||||
reduce_record_len[i] = 0;
|
||||
reduce_record_sparse_len[i] = 0;
|
||||
end_elim[i] = true;
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
for (int k = 0; k < reductors / 64; k++)
|
||||
temp_block[k] ^= R[j]->mono.data[reductor_block + k];
|
||||
//temp_block ^= R[j]->mono.data[reductor_block];
|
||||
reduce_record[i][temp_count] = j;
|
||||
temp_count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (!end_elim[i])
|
||||
{
|
||||
if (E[i]->mono.data_len > reductor_block)
|
||||
E[i]->mono.data_len = reductor_block;
|
||||
}
|
||||
}
|
||||
|
||||
parallel_for(0, 96, [&](int index) {
|
||||
for (int i = 0; i < elim_poly_rows; i++)
|
||||
{
|
||||
for (int j = 0; j < reduce_record_len[i]; j++)
|
||||
{
|
||||
int blocks = reductor_block + reductors / 64;
|
||||
E[i]->xor_only(*R[reduce_record[i][j]], index * blocks / 96, (index + 1) * blocks / 96);
|
||||
}
|
||||
}
|
||||
});
|
||||
|
||||
parallel_for(0, 240, [&](int index) {
|
||||
for (int i = index; i < elim_poly_rows; i += 240)
|
||||
{
|
||||
for (int j = 0; j < reduce_record_sparse_len[i]; j++)
|
||||
{
|
||||
xor_sparse(*E[i], *SR[reduce_record_sparse[i][j]]);
|
||||
}
|
||||
}
|
||||
});
|
||||
}
|
||||
if (reductor_block == 0)
|
||||
{
|
||||
for (int i = 0; i < elim_poly_rows; i++)
|
||||
{
|
||||
if (!end_elim[i]) //说明该行已被消到0
|
||||
{
|
||||
E[i]->leader_order = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
inline int get_max_leader()
|
||||
{
|
||||
int max_leader = -1;
|
||||
for (int i = 0; i < elim_poly_rows; i++)
|
||||
{
|
||||
if (E[i] != 0)
|
||||
{
|
||||
if (max_leader < E[i]->leader_order)
|
||||
{
|
||||
max_leader = E[i]->leader_order;
|
||||
}
|
||||
}
|
||||
}
|
||||
return max_leader;
|
||||
}
|
||||
};
|
||||
#endif
|
||||
Loading…
Reference in New Issue