lizixuan/Algo
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add algo files

Browse Source
This commit is contained in:
lzx
2025-07-04 17:47:19 +08:00
parent 155a338571
commit e7a707ea39
122 changed files with 18845 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

192
heavyHash/DiagonalMatrix.h Normal file
View File

@@ -0,0 +1,192 @@
#ifndef _SINGULAR_DIAGONAL_MATRIX_H
#define _SINGULAR_DIAGONAL_MATRIX_H
#include "singular.h"
//#define L(M,N) (M < N ? M : N)
#define L(M,N) (M*N)
#if 1
typedef struct class_DiagonalMatrix DiagonalMatrix_t;
struct class_DiagonalMatrix {
double *pBlock;
double (*operator)(struct class_DiagonalMatrix *p, int i, int j);
void (*release)(struct class_DiagonalMatrix *p);
};
#else
#include <algorithm>
#include <cassert>
#include <cstring>
//namespace singular {
/**
* Diagonal matrix.
*/
template < int M, int N >
class DiagonalMatrix {
public:
enum {
/** Number of diagonal elements. */
L = M < N ? M : N
};
private:
/**
* Memory block for the diagonal elements.
* The ith row and ith column is given by `elements[i]`.
*/
double* pBlock;
public:
/** Initializes a diagonal matrix filled with 0. */
DiagonalMatrix() {
this->pBlock = new double[L];
std::fill(this->pBlock, this->pBlock + L, 0.0);
}
/**
* Initializes a diagonal matrix with given diagonal values.
*
* The diagonal matrix will look like,
* \f[
* \begin{bmatrix}
* \text{values[0]} & & \\
* & \ddots & \\
* & & \text{values[min(M, N)-1]}
* \end{bmatrix}
* \f]
*
* The behavior is undefined if `values` has less than `min(M, N)`
* elements.
*
* @param values
* Diagonal values of the matrix.
*/
explicit DiagonalMatrix(const double values[]) {
this->pBlock = new double[L];
memcpy(this->pBlock, values, sizeof(double) * L);
}
/**
* Steals the memory block from a given diagonal matrix.
*
* @param[in,out] copyee
* Diagonal matrix from which the memory block is to be stolen.
* No loger valid after this call.
*/
#if SINGULAR_RVALUE_REFERENCE_SUPPORTED
DiagonalMatrix(DiagonalMatrix&& copyee) : pBlock(copyee.pBlock) {
copyee.pBlock = nullptr;
}
#else
DiagonalMatrix(const DiagonalMatrix& copyee) : pBlock(copyee.pBlock) {
const_cast< DiagonalMatrix& >(copyee).pBlock = nullptr;
}
#endif
/** Releases the memory block of this diagonal matrix. */
~DiagonalMatrix() {
this->release();
}
/**
* Steals the memory block from a given diagonal matrix.
*
* @param[in,out] copyee
* Diagonal matrix from which the memory block is to be stolen.
* No longer valid after this call.
* @return
* Reference to this diagonal matrix.
*/
#if SINGULAR_RVALUE_REFERENCE_SUPPORTED
DiagonalMatrix& operator =(DiagonalMatrix&& copyee) {
#else
DiagonalMatrix& operator =(const DiagonalMatrix& copyee) {
#endif
this->release();
this->pBlock = copyee.pBlock;
#if SINGULAR_RVALUE_REFERENCE_SUPPORTED
copyee.pBlock = nullptr;
#else
const_cast< DiagonalMatrix& >(copyee).pBlock = nullptr;
#endif
return *this;
}
/**
* Returns a clone of this matrix.
*
* @return
* Clone of this matrix.
*/
inline DiagonalMatrix clone() const {
return DiagonalMatrix(this->pBlock);
}
/**
* Returns the element at a given row and column.
*
* The behavior is undefined,
* - if `i < 0` or `i >= M`,
* - or if `j < 0` or `j >= N`
*
* @param i
* Index of the row to be obtained.
* @param j
* Index of the column to be obtained.
* @return
* Element at the ith row and jth column.
* 0 if `i != j`.
*/
double operator ()(int i, int j) const {
assert(i >= 0 && i < M);
assert(j >= 0 && j < N);
if (i == j) {
return this->pBlock[i];
} else {
return 0.0;
}
}
/**
* Transposes this matrix.
*
* @return
* Transposed matrix.
*/
DiagonalMatrix< N, M > transpose() const {
return DiagonalMatrix< N, M >(this->pBlock);
}
private:
#if SINGULAR_FUNCTION_DELETION_SUPPORTED
/** Copy constructor is not allowed. */
DiagonalMatrix(const DiagonalMatrix& copyee) = delete;
/** Copy assignment is not allowed. */
DiagonalMatrix& operator =(const DiagonalMatrix& copyee) = delete;
#elif SINGULAR_RVALUE_REFERENCE_SUPPORTED
/** Copy constructor is not allowed. */
DiagonalMatrix(const DiagonalMatrix& copyee) {}
/** Copy assignment is not allowed. */
DiagonalMatrix& operator =(const DiagonalMatrix& copyee) {
return *this;
}
#endif
/**
* Releases the memory block of this matrix.
* Has no effect if the memory block has already been released.
*/
inline void release() {
delete[] this->pBlock;
this->pBlock = nullptr;
}
};
//}
#endif
#endif

14
heavyHash/Makefile Normal file
View File

@@ -0,0 +1,14 @@
SRCS = heavyhash.c obtc.c sha3.c
OBJS = $(SRCS:.c=.o)
CC = gcc
CCFLAGS = -Wall
libkas.a:$(OBJS)
ar -rv libkas.a $(OBJS)
%.o:%.c
$(CC) $(CCFLAGS) -c $< -o $@
clean:
rm -rf *.o *.a

25
heavyHash/Matrix.h Normal file
View File

@@ -0,0 +1,25 @@
#ifndef _SINGULAR_MATRIX_H
#define _SINGULAR_MATRIX_H
#include "singular.h"
#include "Vector.h"
//#include <algorithm>
//#include <cstring>
//#include <iostream>
typedef struct class_Matrix Matrix_t;
struct class_Matrix {
double* pBlock;
Matrix_t (*clone)(struct class_Matrix *p);
void (*filledwith)(struct class_Matrix *p,const double values[]);
double (*operator)(struct class_Matrix *p, int i, int j);
Vector_t (*row)(struct class_Matrix *p, int i);
Vector_t (*column)(struct class_Matrix *p, int j);
void (*release)(struct class_Matrix *p);
};
#endif

19
heavyHash/Reflector.h Normal file
View File

@@ -0,0 +1,19 @@
#ifndef _SINGULAR_REFLECTOR_H
#define _SINGULAR_REFLECTOR_H
#include "Matrix.h"
#include "singular.h"
typedef struct class_Reflector Reflector_t;
struct class_Reflector {
Vector_t u;
double gamma;
size_t L;
double* ptr;
};
#endif

17
heavyHash/Rotator.h Normal file
View File

@@ -0,0 +1,17 @@
#ifndef _SINGULAR_ROTATOR_H
#define _SINGULAR_ROTATOR_H
#include "Matrix.h"
#include "singular.h"
typedef struct class_Rotator Rotator_t;
struct class_Rotator {
double elements[4];
double (*operator)(struct class_Rotator *p, int i, int j);
void (*applyFromLeftTo)(struct class_Rotator *p, Matrix_t rhs, int k);
void (*applyFromRightTo)(struct class_Rotator *p, Matrix_t rhs, int k);
};
#endif

45
heavyHash/Svd.h Normal file
View File

@@ -0,0 +1,45 @@
#ifndef _SINGULAR_SVD_H
#define _SINGULAR_SVD_H
#include "DiagonalMatrix.h"
#include "Matrix.h"
#include "Reflector.h"
#include "Rotator.h"
//#include "singular.h"
//#include <algorithm>
//#include <cassert>
//#include <tuple>
typedef struct Svd Svd_t;
struct Svd {
//USV decomposeUSV(const Matrix< M, N >& m)
bool (*isFullRank)(Svd_t *p, DiagonalMatrix_t singularValues, const int size);
};
typedef struct class_BidiagonalMatrix BidiagonalMatrix_t;
struct class_BidiagonalMatrix {
double* pBlock;
double (*operator)(struct class_BidiagonalMatrix *p, int i, int j);
double (*applyFirstRotatorFromRight)(struct class_BidiagonalMatrix *p, Rotator_t *r);
double (*applyRotatorFromRight)(struct class_BidiagonalMatrix *p, Rotator_t *r, int n, double bulge);
double (*applyRotatorFromLeft)(struct class_BidiagonalMatrix *p, Rotator_t *r, int n, double bulge);
BidiagonalMatrix_t (*bidiagonalize)(struct class_BidiagonalMatrix *p, Matrix_t m);
void (*doFrancis)(struct class_BidiagonalMatrix *m,int n);
double (*calculateShift)(struct class_BidiagonalMatrix *m, int n);
void (*releases)(struct class_BidiagonalMatrix *p);
};
void BidiagonalMatrix_doFrancis(BidiagonalMatrix_t *m, int n);
double BidiagonalMatrix_calculateShift(BidiagonalMatrix_t *m, int n);
double BidiagonalMatrix_applyRotatorFromLeft(BidiagonalMatrix_t *ptr, Rotator_t *r, int n, double bulge);
double BidiagonalMatrix_applyRotatorFromRight(BidiagonalMatrix_t *ptr, Rotator_t *r, int n, double bulge);
double BidiagonalMatrix_applyFirstRotatorFromRight(BidiagonalMatrix_t *p, Rotator_t *r);
double BidiagonalMatrix_operator(BidiagonalMatrix_t *p, int i, int j);
void BidiagonalMatrix_release(BidiagonalMatrix_t *p);
void BidiagonalMatrix_init(BidiagonalMatrix_t *p, Matrix_t *m);
void BidiagonalMatrix_def(BidiagonalMatrix_t *p);
BidiagonalMatrix_t BidiagonalMatrix_bidiagonalize(BidiagonalMatrix_t *p, Matrix_t m);
#endif

22
heavyHash/Vector.h Normal file
View File

@@ -0,0 +1,22 @@
#ifndef _SINGULAR_VECTOR_H
#define _SINGULAR_VECTOR_H
#include <stddef.h>
#include "singular.h"
typedef struct class_Vector Vector_t;
struct class_Vector {
double* pBlock;
size_t len;
ptrdiff_t delta;
double* ptr;
void (*move)(struct class_Vector *p, ptrdiff_t delta);
double (*operator)(struct class_Vector *p, size_t idx);
Vector_t (*slice)(struct class_Vector *p, size_t start);
};
#endif

150
heavyHash/heavyhash.c Normal file
View File

@@ -0,0 +1,150 @@
#include "sha3.h"
#include "obtc.h"
void CSHA3_256_Write(CSHA3_256 *p, const unsigned char* data, size_t len) {
sha3_update(&p->context, data, len);
//return *this;
}
void CSHA3_256_Finalize(CSHA3_256 *p, unsigned char hash[OUTPUT_SIZE]) {
sha3_final(hash, &p->context);
}
/*void CSHA3_256_Reset(Obtc_t *Obtc, CSHA3_256 *p) {
sha3_init(Obtc,&p->context, OUTPUT_SIZE);
//return *this;
}*/
void CSHA3_256_init(Obtc_t *Obtc, CSHA3_256 *p) {
sha3_init(Obtc, &p->context, OUTPUT_SIZE);
p->Write = CSHA3_256_Write;
p->Finalize = CSHA3_256_Finalize;
//p->Reset = CSHA3_256_Reset;
}
void CSHA3_256_CSHA3_256(Obtc_t *Obtc,CSHA3_256 *p) {
sha3_init(Obtc,&p->context, OUTPUT_SIZE);
}
void CHeavyHash_Write(CHeavyHash *p, const unsigned char* data, size_t len) {
p->hasher.Write(&p->hasher,data, len);
//sha3_update(&CSHA3_256_p.context, data, len);
//CSHA3_256_Write(&CSHA3_256_p, data, OUTPUT_SIZE);
}
void CHeavyHash_Finalize(Obtc_t *Obtc, CHeavyHash *p, unsigned char hash[OUTPUT_SIZE]) {
uint256 hash_first;
uint8_t a[32];
p->hasher.Finalize(&p->hasher,&Obtc->g_hash_first.bb.data[0]);
memcpy(a,&Obtc->g_hash_first.bb.data[0],32);
uint256 product = MultiplyUsing4bitPrecision(p->matrix, Obtc->g_hash_first);
uint256 hash_xored;
for (size_t i = 0; i < OUTPUT_SIZE; ++i) {
//hash_xored.begin()[i] = hash_first.begin()[i] ^ product.begin()[i];
hash_xored.bb.data[i] = Obtc->g_hash_first.bb.data[i] ^ product.bb.data[i];
}
uint8_t temp[200]={
0x16,0x19,0x32,0x7d,0x10,0xb9,0xda,0x35,0x54,0x9a,0xe0,0x31,0x2f,0x9f,0xc6,0x15,0x92,0xbb,0x39,0x9d,
0xb5,0x29,0x0c,0x0a,0x47,0xc3,0x9f,0x67,0x51,0x12,0xc2,0x2e,0xc7,0x76,0xc5,0x04,0x84,0x81,0xb9,0x57,
0xb9,0x92,0xf2,0xd3,0x7b,0x34,0xca,0x58,0xea,0x8f,0xdb,0x80,0xba,0xc4,0x6d,0x39,0x7e,0x8f,0x1d,0xb1,
0x77,0x65,0xcc,0x07,0x87,0xe9,0x61,0xb0,0x36,0xbc,0x94,0x16,0x77,0x4c,0x86,0x83,0x54,0x34,0xf2,0xb0,
0x4e,0xf7,0x4b,0x3a,0x99,0xcd,0xb0,0x44,0x2e,0xc6,0x5b,0xd3,0x56,0x24,0x93,0xe4,0x6c,0x6b,0x7d,0x01,
0xa7,0x69,0xcc,0x3d,0xd3,0x1f,0x4c,0xc3,0x54,0xc1,0x8c,0x3f,0xf4,0x31,0xc0,0x5d,0xd0,0xa9,0xa2,0x26,
0xa0,0xbc,0xaa,0x9f,0x79,0x2a,0x3d,0x0c,0x80,0x39,0xf9,0xa6,0x0d,0xcf,0x6a,0x48,0x5e,0x21,0x90,0x40,
0x25,0x0f,0xc4,0x62,0xc1,0x00,0xff,0x2a,0x93,0x89,0x35,0xba,0x72,0xc7,0xd8,0x2e,0x14,0xf3,0x40,0x69,
0xe7,0x20,0xe0,0xdf,0x44,0xee,0xce,0xde,0x11,0xa7,0x5f,0x4c,0x80,0x05,0x64,0x98,0x7a,0x14,0xff,0x48,
0x16,0xc7,0xf8,0xee,0x79,0x62,0x9b,0x0e,0x2f,0x9f,0x42,0x16,0x3a,0xd7,0x4c,0x52,0xb2,0x24,0x85,0x09,
};
for(int i = 0 ;i< 200 ;i++)Obtc->const_data[i] = temp[i];
// CSHA3_256().Write(hash_xored.begin(), OUTPUT_SIZE).Finalize(hash);
CSHA3_256_CSHA3_256(Obtc, &p->hasher);
CSHA3_256_Write(&p->hasher, &hash_xored.bb.data[0], OUTPUT_SIZE);
CSHA3_256_Finalize(&p->hasher, hash) ;
}
void CHeavyHash_Reset(CHeavyHash *p, uint64_t matrix_[64*64]) {
for (int i = 0; i < 64*64; ++i)
p->matrix[i] = matrix_[i];
}
void CHeavyHash_init(Obtc_t *Obtc, CHeavyHash *p, uint64_t matrix_[64*64]){
p->Write = CHeavyHash_Write;
p->Finalize = CHeavyHash_Finalize;
p->Reset = CHeavyHash_Reset;
p->hasher.Write = CSHA3_256_Write;
p->hasher.Finalize = CSHA3_256_Finalize;
//p->hasher.Reset = CSHA3_256_Reset;
sha3_init(Obtc, &p->hasher.context, OUTPUT_SIZE);
for (int i = 0; i < 64*64; ++i)
p->matrix[i] = matrix_[i];
}
void MultiplyMatrices(uint64_t matrix[64*64], uint64_t vector[64], uint64_t product[64]){
for (int i = 0; i < 64; ++i) {
for (int j = 0; j < 64; ++j) {
product[i] += matrix[64*i + j]*vector[j];
}
}
}
uint256 MultiplyUsing4bitPrecision(uint64_t matrix[64*64], const uint256 hash) {
// conversion to matrix with 4 bit values
uint64_t vector[64] = {0};
ConvertTo4BitPrecisionVector(hash, vector);
// perform matrix multiplication
uint64_t product[64] = {0};
MultiplyMatrices(matrix, vector, product);
for (int i = 0; i < 64; ++i) {
product[i] >>= 10;
}
return Convert4bitVectorToUint(product);
}
void ConvertTo4BitPrecisionVector(uint256 bit_sequence, uint64_t vector[64]) {
int index = 0;
int i;
for (i = 0; i < WIDTH; i++) {
vector[index] = bit_sequence.bb.data[i] >> 4;
vector[index+1] = bit_sequence.bb.data[i] & 0xF;
index += 2;
}
}
uint256 Convert4bitVectorToUint(const uint64_t x[64]) {
uint256 bit_sequence;
int index = 0;
int i;
for (i = 0; i < WIDTH; i++) {
bit_sequence.bb.data[i] = ( x[index] << 4) | x[index+1];
index += 2;
}
return bit_sequence;
}

98
heavyHash/heavyhash.h Normal file
View File

@@ -0,0 +1,98 @@
#ifndef OPOW_CRYPTO_HEAVYHASH_H
#define OPOW_CRYPTO_HEAVYHASH_H
#include <stdint.h>
#include <stdlib.h>
#include "sha3.h"
//#include <memory>
//#include "obtc.h"
#define OUTPUT_SIZE 32
typedef struct class_CSHA3_256 CSHA3_256;
struct class_CSHA3_256
{
sha3_ctx_t context;
// static const size_t OUTPUT_SIZE = 32;
//CSHA3_256& Write(const unsigned char* data, size_t len);
void (*Write)(struct class_CSHA3_256 *p, const unsigned char* data, size_t len);
void (*Finalize)(struct class_CSHA3_256 *p, unsigned char hash[OUTPUT_SIZE]);
//CSHA3_256& Reset();
};
typedef struct class_CHeavyHash CHeavyHash;
struct class_CHeavyHash
{
uint64_t matrix[64*64];
CSHA3_256 hasher;
//static const size_t OUTPUT_SIZE = 32;
//explicit CHeavyHash(uint64_t matrix_[64*64]);
//CHeavyHash& Reset(uint64_t matrix_[64*64]);
//CHeavyHash& Write(const unsigned char* data, size_t len);
//void Finalize(unsigned char hash[OUTPUT_SIZE]);
void (*Reset)(struct class_CHeavyHash *p, uint64_t matrix_[64*64]);
void (*Write)(struct class_CHeavyHash *p, const unsigned char* data, size_t len);
void (*Finalize)(struct class_CHeavyHash *p, unsigned char hash[OUTPUT_SIZE]);
};
#if 0
/** A hasher class for SHA3-256. */
class CSHA3_256
{
private:
sha3_ctx_t context;
public:
static const size_t OUTPUT_SIZE = 32;
CSHA3_256();
CSHA3_256& Write(const unsigned char* data, size_t len);
void Finalize(unsigned char hash[OUTPUT_SIZE]);
CSHA3_256& Reset();
};
class CHeavyHash
{
private:
uint64_t matrix[64*64];
CSHA3_256 hasher;
public:
static const size_t OUTPUT_SIZE = 32;
explicit CHeavyHash(uint64_t matrix_[64*64]);
CHeavyHash& Reset(uint64_t matrix_[64*64]);
CHeavyHash& Write(const unsigned char* data, size_t len);
void Finalize(unsigned char hash[OUTPUT_SIZE]);
};
#endif
uint256 MultiplyUsing4bitPrecision(uint64_t matrix[64*64], const uint256 hash);
void ConvertTo4BitPrecisionVector(uint256 bit_sequence, uint64_t vector[64]);
uint256 Convert4bitVectorToUint(const uint64_t x[64]);
//zzj add
/*extern void CSHA3_256_init(struct Obtc_opt *Obtc, CSHA3_256 *p);
void CSHA3_256_CSHA3_256(struct Obtc_opt *Obtc, CSHA3_256 *p);
void CSHA3_256_Write(CSHA3_256 *p, const unsigned char* data, size_t len);
void CSHA3_256_Finalize(CSHA3_256 *p, unsigned char hash[OUTPUT_SIZE]);
//
void CHeavyHash_init(struct Obtc_opt *Obtc, CHeavyHash *p, uint64_t matrix_[64*64]);
void CHeavyHash_Write(CHeavyHash *p, const unsigned char* data, size_t len);
void CHeavyHash_Finalize(struct Obtc_opt *Obtc, CHeavyHash *p, unsigned char hash[OUTPUT_SIZE]);
*/
#endif // OPOW_CRYPTO_HEAVYHASH_H

BIN
heavyHash/heavyhash.o Normal file
View File

Binary file not shown.

BIN
heavyHash/libkas.a Normal file
View File

Binary file not shown.

907
heavyHash/obtc.c Normal file
View File

@@ -0,0 +1,907 @@
//! heavyhash extracted from optical bitcoin
//! 2022 barrystyle
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <search.h>//qsort
#include<time.h>
#include "obtc.h"
#define M 64
#define N 64
bool Is4BitPrecision(const uint64_t matrix[64*64])
{
for (int i = 0; i < 64; ++i) {
for (int j = 0; j < 64; ++j) {
if (matrix[ i*64 + j] > 0xF)
return false;
}
}
return true;
}
double DiagonalMatrix_operator(DiagonalMatrix_t *p, int i, int j)
{
assert(i >= 0 && i < 64);
assert(j >= 0 && j < 64);
if (i == j) {
return p->pBlock[i];
} else {
return 0.0;
}
}
void DiagonalMatrix_release(DiagonalMatrix_t *p)
{
if (p->pBlock != NULL){
free(p->pBlock);
p->pBlock = NULL;
}
}
void DiagonalMatrix_init(DiagonalMatrix_t *p, const double values[])
{
p->pBlock = (double *)malloc(sizeof(double)*M);
//memset(pBlock, 0.0, sizeof(double)*L(64,64));
memcpy(p->pBlock, values, sizeof(double) * M);
p->operator = DiagonalMatrix_operator;
p->release = DiagonalMatrix_release;
}
void DiagonalMatrix_DiagonalMatrix(DiagonalMatrix_t *p)
{
p->operator = DiagonalMatrix_operator;
p->release = DiagonalMatrix_release;
}
//-----------------------------vector-------------------------------//
void vector_move(Vector_t *p, ptrdiff_t delta) {
p->ptr += delta;
}
Vector_t vector_slice(Vector_t v, size_t start) {
//assert(start >= 0 && start <= p->len);
Vector_t v_tmp;
v_tmp.pBlock = v.pBlock + start * v.delta;
v_tmp.len = v.len - start;
v_tmp.delta = v.delta;
return v_tmp;
}
double Vector_column_operator(Vector_t *p, size_t idx){
return p->pBlock[idx * p->delta];
}
double Vector_row_operator(Vector_t *p, size_t idx){
return p->pBlock[idx * p->delta];
}
void Vector_sync(Matrix_t *p, size_t idx, Vector_t vec, int offset){
for(int i = 0; i < vec.len; i++){
p->pBlock[idx+(offset+i)*N] = vec.pBlock[i];
}
}
void Vector_row_sync(Matrix_t *p, size_t idx, Vector_t vec, int offset){
for(int i = 0; i < vec.len; i++){
p->pBlock[offset+idx*N+i] = vec.pBlock[i];
}
}
//-----------------------------Martrix-------------------------------//
Matrix_t Matrix_clone(Matrix_t *p)
{
Matrix_t m;
m.pBlock = (double *)malloc(sizeof(double)*L(64,64));
memcpy(m.pBlock, p->pBlock, sizeof(double)*L(64,64));
return m;
}
void Matrix_filledwith(Matrix_t *p, const double values[])
{
//p->pBlock = (double *)malloc(sizeof(double)*L(64,64));
//memset(pBlock, 0.0, sizeof(double)*L(64,64));
memcpy(p->pBlock, values, sizeof(double) * L(64,64));
}
double Matrix_operator(Matrix_t *p, int i, int j)
{
assert(i >= 0 && i < N);
assert(j >= 0 && j < N);
return p->pBlock[i*N+j];
}
Vector_t Matrix_row(Matrix_t *p, int i)
{
Vector_t vec_tmp;
vec_tmp.len = N;
vec_tmp.delta = 1;
vec_tmp.pBlock = p->pBlock + i*N;
//return Vector< const double >(this->pBlock + i * N, N, 1);
return vec_tmp;
}
Vector_t Matrix_column(Matrix_t *p, int j)
{
Vector_t vec_tmp;
vec_tmp.len = M;
vec_tmp.delta = N;
vec_tmp.pBlock = p->pBlock + j;
return vec_tmp;
//return Vector< double >(this->pBlock + j, M, N);
}
void Matrix_release(Matrix_t *p)
{
if (p->pBlock != NULL){
free(p->pBlock);
p->pBlock = NULL;
}
}
void Matrix_init(Matrix_t *p)
{
p->pBlock = (double *)malloc(sizeof(double)*L(64,64));
memset(p->pBlock, 0.0, sizeof(double)*L(64,64));
//memcpy(p->pBlock, values, sizeof(double) * L(64,64));
}
void Matrix_def(Matrix_t *p)
{
//p->clone = Matrix_clone;
p->filledwith = Matrix_filledwith;
p->operator = Matrix_operator;
p->row = Matrix_row;
p->column = Matrix_column;
p->release = Matrix_release;
}
//-----------------------------Rotator-------------------------------//
double max(double a, double b)
{
return a > b ? a : b;
}
double Rotator_operator(Rotator_t *p, int i, int j){
assert(0 <= i && i < 2);
assert(0 <= j && j < 2);
return p->elements[i * 2 + j];
}
void Rotator_init(Rotator_t *p, double x1, double x2)
{
// normalizes by the maximum magnitude
// to avoid harmful underflow and overflow
double mx = max(fabs(x1), fabs(x2));
x1 /= mx;
x2 /= mx;
double norm = sqrt(x1 * x1 + x2 * x2);
double cs = x1 / norm;
double sn = x2 / norm;
p->elements[0] = cs;
p->elements[1] = -sn;
p->elements[2] = sn;
p->elements[3] = cs;
p->operator = Rotator_operator;
}
//-----------------------------Reflector-------------------------------//
void Reflector_transform(Reflector_t *p, double u0, size_t len){
int i;
for (i = 0; i < len; i++){
p->u.pBlock[i] = p->u.pBlock[i] /u0;
}
}
void Reflector_transform_left(Reflector_t *src1, Vector_t src2, Vector_t dst, double gUM, size_t len){
int i;
for (i = 0; i < len; i++){
dst.pBlock[i] = src2.pBlock[i] - src1->u.pBlock[i] * gUM;
}
}
void Reflector_transform_right(Reflector_t *src1, Vector_t src2, Vector_t dst, double gMU, size_t len){
int i;
for (i = 0; i < len; i++){
dst.pBlock[i] = src2.pBlock[i] - gMU * src1->u.pBlock[i];
}
}
void Reflector_init(Reflector_t *p, Vector_t v) {
//assert(v.size() > 0 && v.size() <= L);
//const size_t N = v.size();
//const size_t p->L = sizeof(v)/sizeof(double);
p->L = v.len;
p->u.pBlock = (double *)malloc(sizeof(double)*v.len);
memcpy(p->u.pBlock, v.pBlock, sizeof(double)*v.len);
// normalizes elements by the maximum amplitude
// to avoid harmful underflow and overflow
double mx = 0.0;
for (size_t i = 0; i < p->L; ++i) {
mx = max(fabs(p->u.pBlock[i]), mx);
}
if (mx > 0.0) {
// calculates the normalized norm
double tau = 0.0;
for (size_t i = 0; i < p->L; ++i) {
double x = p->u.pBlock[i] / mx;
p->u.pBlock[i] = x;
tau += x * x;
}
tau = sqrt(tau);
// tau's sign should be the same as the first element in `u`
if (p->u.pBlock[0] < 0.0) {
tau = -tau;
}
double u0 = p->u.pBlock[0] + tau;
p->u.pBlock[0] = u0;
Reflector_transform(p, u0, p->L);
p->gamma = u0 / tau;
} else {
// v is a zero vector
p->gamma = 0.0;
memset(p->u.pBlock, 0.0, p->L);
}
}
void Reflector_release(Reflector_t *p){
if (p->u.pBlock != NULL){
free(p->u.pBlock);
p->u.pBlock = NULL;
}
}
double inner_product(double *a,double *b,int n){
int i;
double sum = 0.0;
for(i = 0; i < n; i++)
{
sum += (*(a+i))*(*(b+i));
}
return sum;
}
Matrix_t Reflector_applyFromLeftTo(Reflector_t *p, Matrix_t m){
// H * m = m - gamma * u * u^T * m
Matrix_t m2 = Matrix_clone(&m);//m->clone(m);
Vector_t vec_m;
Vector_t vec_m2;
int offset = N - p->L;
for (int i = 0; i < N; ++i) {
// caches gamma * u^T * m
vec_m = Matrix_column(&m, i);
Vector_t srcColumn = vector_slice(vec_m, offset);
double v_src_column[srcColumn.len];
for(size_t i = 0; i < srcColumn.len; i++){
v_src_column[i] = Vector_column_operator(&srcColumn, i);
}
srcColumn.pBlock = v_src_column;
double gUM = inner_product(p->u.pBlock, srcColumn.pBlock, p->L);
//Vector< const double > srcColumn = m->column(m, i).slice(offset);
gUM *= p->gamma;
// H * m = m - u * gUM
vec_m2 = Matrix_column(&m2, i);
Vector_t dstColumn = vector_slice(vec_m2, offset);
double v_dstcolumn[dstColumn.len];
for(size_t i = 0; i < dstColumn.len; i++){
v_dstcolumn[i] = Vector_column_operator(&dstColumn, i);
}
dstColumn.pBlock = v_dstcolumn;
Reflector_transform_left(p, srcColumn, dstColumn, gUM, p->L);
Vector_sync(&m2, i, dstColumn, offset);
}
Matrix_release(&m);
return m2;
}
Matrix_t Reflector_applyFromRightTo(Reflector_t *p, Matrix_t m){
// m * H = m - m * gamma * u * u^T
Matrix_t m2 = Matrix_clone(&m);
Vector_t vec_m;
Vector_t vec_m2;
int offset = 64 - p->L;
for (int i = 0; i < M; ++i) {
// caches gamma * m * u
vec_m = Matrix_row(&m, i);
Vector_t srcRow = vector_slice(vec_m, offset);
double v_src_row[srcRow.len];
for(size_t j = 0; j< srcRow.len; j++){
v_src_row[j] = Vector_row_operator(&srcRow, j);
}
srcRow.pBlock = v_src_row;
double gMU = inner_product(p->u.pBlock, srcRow.pBlock, p->L);
gMU *= p->gamma;
// m * H = m - gMU * u^T
vec_m2 = Matrix_row(&m2, i);
Vector_t dstRow = vector_slice(vec_m2, offset);
double v_dstrow[dstRow.len];
for(size_t j = 0; j < dstRow.len; j++){
v_dstrow[j] = Vector_row_operator(&dstRow, j);
}
dstRow.pBlock = v_dstrow;
Reflector_transform_right(p ,srcRow, dstRow, gMU, p->L);
Vector_row_sync(&m2, i, dstRow, offset);
}
Matrix_release(&m);
return m2;
}
//-----------------------------Svd-------------------------------//
int cmp_double(const void* e1, const void* e2)
{
if ((*(double*)e2 - *(double*)e1) > 0.00000)
return 1;
else if ((*(double*)e2 - *(double*)e1) == 0.000000)
return 0;
else
return -1;
}
DiagonalMatrix_t Svd_decomposeUSV(BidiagonalMatrix_t *p, Matrix_t *m) {
const int MAX_ITERATIONS = N * 10;
// allocates matrices
Matrix_t m1 = Matrix_clone(m);
Matrix_def(&m1);
// bidiagonalizes a given matrix
BidiagonalMatrix_t m2 = p->bidiagonalize(p, m1);
// repeats Francis iteration
int iteration = 0;
int n = N;
while (n >= 2) {
// processes the n-1 x n-1 submatrix
// if the current n x n submatrix has converged
double bn = m2.operator(&m2, n - 1, n - 1);
if (bn == 0.0 || fabs(m2.operator(&m2, n - 2, n - 1) / bn) < 1.0e-15) {
--n;
} else {
// aborts if too many iterations
++iteration;
if (iteration > MAX_ITERATIONS) {
break;
}
m2.doFrancis(&m2, n);
}
}
// copies the diagonal elements
// and makes all singular values positive
double ss[N];
for (int i = 0; i < N; ++i) {
if (m2.operator(&m2, i, i) < 0) {
ss[i] = -m2.operator(&m2, i, i);
// inverts the sign of the right singular vector
//Vector< double > vi = v.column(i);
//std::transform(
// vi.begin(), vi.end(), vi.begin(),
// [](double x) {
// return -x;
// });
} else {
ss[i] = m2.operator(&m2, i, i);
}
}
// sorts singular values in descending order if necessary
int shuffle[M]; // M >= N
bool sortNeeded = false;
for (int i = 0; i < M; ++i) {
shuffle[i] = i;
sortNeeded = sortNeeded || (i < N - 1 && ss[i] < ss[i + 1]);
}
m1.release(&m1);
BidiagonalMatrix_release(p);
DiagonalMatrix_t dm;
if (sortNeeded) {
// shuffles the N (<= M) singular values
qsort(ss, N,sizeof(double), cmp_double);
double ss2[M];
memcpy(ss2, ss, M*sizeof(double));
DiagonalMatrix_init(&dm, ss2);
return dm;
} else {
DiagonalMatrix_init(&dm, ss);
return dm;
}
}
bool Svd_isFullRank(DiagonalMatrix_t *p, const int size) {
const double round_off = 1.000009e-12;
for (int i = 0; i < size; ++i) {
if (fabs( p->operator(p, i, i) ) < round_off){
p->release(p);
return false;
}
}
p->release(p);
return true;
}
//-----------------------------BidiagonalMatrix_t-------------------------------//
BidiagonalMatrix_t BidiagonalMatrix_bidiagonalize(BidiagonalMatrix_t *p, Matrix_t m)
{
assert(M >= N);
Vector_t vec_m;
Vector_t vec_m2;
for (int i = 0; i < N; ++i) {
Reflector_t rU;
vec_m = Matrix_column(&m, i);
Vector_t column_slice = vector_slice(vec_m, i);
// applies a householder transform to the column vector i
double v_column[column_slice.len];
for(size_t i = 0; i < column_slice.len; i++){
v_column[i] = Vector_column_operator(&column_slice, i);
}
column_slice.pBlock = v_column;
Reflector_init(&rU, column_slice);
m = Reflector_applyFromLeftTo(&rU, m);
Reflector_release(&rU);
//u = rU.applyFromRightTo(u); // U1^T*U0^T = U0*U1
if (i < N - 1) {
// applies a householder transform to the row vector i + 1
//Reflector< N > rV(m.row(i).slice(i + 1));
Reflector_t rV;
vec_m2 = Matrix_row(&m, i);
Vector_t row_slice = vector_slice(vec_m2, i+1);
double v_row[row_slice.len];
for(size_t i = 0; i < row_slice.len; i++){
v_row[i] = Vector_row_operator(&row_slice, i);
}
row_slice.pBlock = v_row;
Reflector_init(&rV, row_slice);
m = Reflector_applyFromRightTo(&rV, m);
//m = rV.applyFromRightTo(m);
//v = rV.applyFromRightTo(v);
Reflector_release(&rV);
}
}
BidiagonalMatrix_init(p, &m);
return *p;
}
void BidiagonalMatrix_release(BidiagonalMatrix_t *p)
{
if (p->pBlock != NULL){
free(p->pBlock);
p->pBlock = NULL;
}
}
double BidiagonalMatrix_operator(BidiagonalMatrix_t *p, int i, int j)
{
assert(i >= 0 && i < M);
assert(j >= 0 && j < N);
if (i == j) {
return p->pBlock[2 * i];
} else if (i + 1 == j) {
return p->pBlock[2 * i + 1];
} else {
return 0.0;
}
}
double BidiagonalMatrix_applyFirstRotatorFromRight(BidiagonalMatrix_t *p, Rotator_t *r)
{
double b1 = p->pBlock[0];
double g1 = p->pBlock[1];
double b2 = p->pBlock[2];
double r11 = Rotator_operator(r, 0, 0);//r->operator(r, 0, 0);
double r12 = Rotator_operator(r, 0, 1);//r->operator(r, 0, 1);
double r21 = Rotator_operator(r, 1, 0);//r->operator(r, 1, 0);
double r22 = Rotator_operator(r, 1, 1);//r->operator(r, 1, 1);
//Rotator_operator
p->pBlock[0] = b1 * r11 + g1 * r21;
p->pBlock[1] = b1 * r12 + g1 * r22;
p->pBlock[2] = b2 * r22;
return b2 * r21;
}
double BidiagonalMatrix_applyRotatorFromRight(BidiagonalMatrix_t *ptr, Rotator_t *r, int n, double bulge)
{
double* p = ptr->pBlock + n * 2;
double g0 = p[-1];
double b1 = p[0];
double g1 = p[1];
double b2 = p[2];
double r11 = r->operator(r, 0, 0);
double r12 = r->operator(r, 0, 1);
double r21 = r->operator(r, 1, 0);
double r22 = r->operator(r, 1, 1);
p[-1] = g0 * r11 + bulge * r21;
p[0] = b1 * r11 + g1 * r21;
p[1] = b1 * r12 + g1 * r22;
p[2] = b2 * r22;
return b2 * r21;
}
double BidiagonalMatrix_applyRotatorFromLeft(BidiagonalMatrix_t *ptr, Rotator_t *r, int n, double bulge)
{
double* p = ptr->pBlock + n * 2;
double b1 = p[0];
double g1 = p[1];
double b2 = p[2];
double r11 = r->operator(r, 0, 0);
double r12 = r->operator(r, 0, 1);
double r21 = r->operator(r, 1, 0);
double r22 = r->operator(r, 1, 1);
p[0] = r11 * b1 + r21 * bulge;
p[1] = r11 * g1 + r21 * b2;
p[2] = r12 * g1 + r22 * b2;
double newBulge;
if (n < N - 2) {
double g2 = p[3];
newBulge = r21 * g2;
p[3] = r22 * g2;
} else {
newBulge = 0.0;
}
return newBulge;
}
double BidiagonalMatrix_calculateShift(BidiagonalMatrix_t *m, int n)
{
assert(M >= N);
assert(n >= 2);
double b1 = m->operator(m, n - 2, n - 2);
double b2 = m->operator(m, n - 1, n - 1);
double g1 = m->operator(m, n - 2, n - 1);
// solves lambda^4 - d*lambda^2 + e = 0
// where
// d = b1^2 + b2^2 + g1^2
// e = b1^2 * b2^2
// chooses lambda (rho) closest to b2
double rho;
double d = b1 * b1 + b2 * b2 + g1 * g1;
double e = b1 * b1 * b2 * b2;
// lambda^2 = (d +- sqrt(d^2 - 4e)) / 2
// so, f = d^2 - 4e must be positive
double f = d * d - 4 * e;
if (f >= 0) {
f = sqrt(f);
// lambda = +-sqrt(d +- f) (d >= 0, f >= 0)
// if d > f, both d+f and d-f have real square roots
// otherwise considers only d+f
if (d > f) {
// lets l1 > l2
double l1 = sqrt((d + f) * 0.5);
double l2 = sqrt((d - f) * 0.5);
// if b2 >= 0, chooses a positive shift
// otherwise chooses a negative shift
if (b2 >= 0) {
if (fabs(b2 - l1) < fabs(b2 - l2)) {
rho = l1;
} else {
rho = l2;
}
} else {
if (fabs(b2 + l1) < fabs(b2 + l2)) {
rho = -l1;
} else {
rho = -l2;
}
}
} else {
double l1 = sqrt((d + f) * 0.5);
if (fabs(b2 - l1) <= fabs(b2 + l1)) {
rho = l1;
} else {
rho = -l1;
}
}
} else {
// no solution. chooses b2 as the shift
rho = b2;
}
return rho;
}
void BidiagonalMatrix_doFrancis(BidiagonalMatrix_t *m, int n)
{
assert(M >= N);
assert(n >= 2);
// calculates the shift
double rho = m->calculateShift(m, n);
// applies the first right rotator
double b1 = m->operator(m, 0, 0);
double g1 = m->operator(m, 0, 1);
double mx = max(fabs(rho), max(fabs(b1), fabs(g1)));
rho /= mx;
b1 /= mx;
g1 /= mx;
//Rotator_t r0(b1 * b1 - rho * rho, b1 * g1);
Rotator_t r0;
Rotator_init(&r0, b1 * b1 - rho * rho, b1 * g1);
double bulge = m->applyFirstRotatorFromRight(m, &r0);
//v = r0.applyFromRightTo(&r0, v, 0);
// applies the first left rotator
Rotator_t r1;
Rotator_init(&r1, m->operator(m, 0, 0), bulge);
//Rotator_t r1(m(0, 0), bulge);
bulge = m->applyRotatorFromLeft(m, &r1, 0, bulge);
//u = r1.applyFromRightTo(&r1, u, 0); // U1^T*U0^T = U0*U1
for (int i = 1; i + 1 < n; ++i) {
// calculates (i+1)-th right rotator
//Rotator rV(m(i - 1, i), bulge);
Rotator_t rV;
Rotator_init(&rV, m->operator(m, i - 1, i), bulge);
bulge = m->applyRotatorFromRight(m, &rV, i, bulge);
//v = rV.applyFromRightTo(&rV, v, i);
// calculates (i+1)-th left rotator
//Rotator rU(m(i, i), bulge);
Rotator_t rU;
Rotator_init(&rU, m->operator(m, i, i), bulge);
bulge = m->applyRotatorFromLeft(m, &rU, i, bulge);
//u = rU.applyFromRightTo(rU, u, i); // U1^T*U0^T = U0*U1
}
}
void BidiagonalMatrix_def(BidiagonalMatrix_t *p)
{
p->applyFirstRotatorFromRight = BidiagonalMatrix_applyFirstRotatorFromRight;
p->applyRotatorFromLeft = BidiagonalMatrix_applyRotatorFromLeft;
p->applyRotatorFromRight = BidiagonalMatrix_applyRotatorFromRight;
p->bidiagonalize = BidiagonalMatrix_bidiagonalize;
p->calculateShift = BidiagonalMatrix_calculateShift;
p->doFrancis = BidiagonalMatrix_doFrancis;
p->operator = BidiagonalMatrix_operator;
p->releases = BidiagonalMatrix_release;
}
void BidiagonalMatrix_init(BidiagonalMatrix_t *p, Matrix_t *m)
{
assert(M >= N);
int len;
len = 2 * N - 1;
p->pBlock = (double *)malloc(sizeof(double)*len);
memset(p->pBlock, 0.0,sizeof(double)*len);
for (int i = 0; i < N; ++i) {
p->pBlock[i * 2] = Matrix_operator(m, i, i);//m->operator(m, i, i);
if (i < N - 1) {
p->pBlock[i * 2 + 1] = Matrix_operator(m, i, i + 1);//m->operator(m, i, i + 1);
}
}
}
bool IsFullRank(const uint64_t matrix_[64*64])
{
double matrix__ [64*64];
// Matrix<64, 64> matrix;
for (int i = 0; i < 64; ++i) {
for (int j = 0; j < 64; ++j) {
matrix__[64*i + j] = (double) matrix_[64*i + j];
}
}
DiagonalMatrix_t dm;
Matrix_t mt;
BidiagonalMatrix_t bt;
DiagonalMatrix_init(&dm, matrix__);
//matrix.fill(matrix__);
Matrix_init(&mt);
Matrix_def(&mt);
mt.filledwith(&mt, matrix__);
BidiagonalMatrix_def(&bt);
DiagonalMatrix_t usv = Svd_decomposeUSV(&bt, &mt);
DiagonalMatrix_t singularValues = usv;
mt.release(&mt);
dm.release(&dm);
//DiagonalMatrix_release(&dm);
return Svd_isFullRank(&usv,64);
}
uint64_t GetUint64_t(uint8_t *data, int pos)
{
const uint8_t* ptr = data + pos * 8;
return ((uint64_t)ptr[0]) | \
((uint64_t)ptr[1]) << 8 | \
((uint64_t)ptr[2]) << 16 | \
((uint64_t)ptr[3]) << 24 | \
((uint64_t)ptr[4]) << 32 | \
((uint64_t)ptr[5]) << 40 | \
((uint64_t)ptr[6]) << 48 | \
((uint64_t)ptr[7]) << 56;
}
void XoShiRo256PlusPlus_init(Obtc_t *Obtc, uint64_t *s, uint256 seed) {
for (int i = 0; i < 4; ++i) {
//p->s[i] = seed.GetUint64(i);
s[i] = GetUint64_t(Obtc->data_r,i);
}
}
uint64_t RotateLeft64(const uint64_t x, int k) {
return (x << k) | (x >> (64 - k));
}
uint64_t XoShiRo256PlusPlus_operator(uint64_t *s){
const uint64_t result = RotateLeft64(s[0] + s[3], 23) + s[0];
const uint64_t t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = RotateLeft64(s[3], 45);
return result;
}
void GenerateHeavyHashMatrix_t(Obtc_t *Obtc, uint256 matrix_seed, uint64_t matrix[64*64])
{
XoShiRo256PlusPlus_init(Obtc, Obtc->ss, matrix_seed);
do {
for (int i = 0; i < 64; ++i) {
for (int j = 0; j < 64; j += 16) {
uint64_t value = XoShiRo256PlusPlus_operator(Obtc->ss);//generator();
for (int shift = 0; shift < 16; ++shift) {
matrix[64*i + j + shift] = (value >> (4 * shift)) & 0xF;
}
}
}
//} while (!Is4BitPrecision(matrix) || !IsFullRank(matrix));
}while(!Is4BitPrecision(matrix));
}
void serialize_heavyhash(Obtc_t *Obtc, uint64_t matrix[64*64], const char* in, char* out, int len)
{
uint8_t temp[200]={
0x02,0xb9,0x7c,0x78,0x6f,0x82,0x43,0x83,0x5d,0x11,0x29,0xcf,0x82,0xaf,0xa5,0xbc,0xb1,0xfc,0xce,0x9c,
0xe7,0x8b,0x52,0x72,0x48,0xb0,0x94,0x27,0xa8,0x74,0x2e,0xdb,0x89,0xca,0x4e,0x84,0x9b,0xce,0xcf,0x4a,
0xd1,0x02,0x57,0x41,0x05,0x09,0x5f,0x8d,0xba,0x1d,0xe5,0xe4,0x45,0x16,0x68,0xe4,0xc1,0xa2,0x02,0x1d,
0x56,0x3b,0xb1,0x42,0x8f,0x06,0xdd,0x1c,0x7a,0x2f,0x85,0x1a,0x34,0x85,0x54,0x90,0x64,0xa3,0x6a,0x46,
0xb2,0x1a,0x60,0x1f,0x85,0xb4,0xb2,0x23,0xe6,0xc8,0x5d,0x8f,0x82,0xe9,0xda,0x89,0xec,0x70,0xf1,0xa4,
0x25,0xb1,0x37,0x15,0x44,0xe3,0x67,0x87,0x5b,0x29,0x91,0x52,0x0f,0x96,0x07,0x05,0x40,0xf1,0x4a,0x0e,
0x2e,0x65,0x1c,0x3c,0x43,0x28,0x5f,0xf0,0xf8,0xeb,0xf1,0x33,0x88,0x66,0x31,0x40,0x77,0x6b,0xf6,0x0c,
0x78,0x9b,0xc2,0x9c,0x18,0x3a,0x98,0x1e,0xad,0x41,0x5b,0x10,0x4a,0xef,0x61,0xd6,0x29,0xdc,0xe2,0x46,
0x7b,0x2f,0xaf,0xca,0x87,0x5e,0x2d,0x65,0x1b,0xa5,0xa4,0xa3,0xf5,0x98,0x69,0xa0,0x1e,0x5f,0x2e,0x72,
0x0e,0xfb,0x44,0xd2,0x29,0xbf,0x88,0x55,0xb7,0x02,0x7e,0x3c,0x11,0x3c,0xff,0x0d,0xa1,0xf6,0xd8,0x3d
};
for(int i = 0 ;i< 200 ;i++)Obtc->const_data[i] = temp[i];
CHeavyHash_init(Obtc, &Obtc->CHeavyHash_p, matrix);
CHeavyHash_Write(&Obtc->CHeavyHash_p, (const unsigned char*)in, len);
CHeavyHash_Finalize(Obtc, &Obtc->CHeavyHash_p, (unsigned char*)out);
}
void opticalbtc_hash(const char* in, char* out, int len)
{
uint8_t *ptr = (uint8_t*) in;
uint256 seed, hashprev;
uint64_t matrix[64*64];
Obtc_t Obtc;
CSHA3_256_init(&Obtc, &Obtc.CSHA3_256_p);
memcpy(Obtc.data_r,ptr, 32);
GenerateHeavyHashMatrix_t(&Obtc, seed, matrix);
serialize_heavyhash(&Obtc, matrix, in, out, len);
}

52
heavyHash/obtc.h Normal file
View File

@@ -0,0 +1,52 @@
#ifndef OBTC_H
#define OBTC_H
#include "uint256.h"
#include "xoshiro256pp.h"
#include "Svd.h"
#include "DiagonalMatrix.h"
#include "Matrix.h"
#include "Rotator.h"
#include "heavyhash.h"
typedef struct Obtc_opt Obtc_t;
struct Obtc_opt{
uint8_t data_r[32];
uint64_t ss[4];
uint8_t const_data[200];
CSHA3_256 CSHA3_256_p;
CHeavyHash CHeavyHash_p;
uint256 g_hash_first;
XoShiRo256PlusPlus_t *xo;
DiagonalMatrix_t g_DiagonalMatrix;
};
//struct Obtc_opt;
bool Is4BitPrecision(const uint64_t matrix[64*64]);
bool IsFullRank(const uint64_t matrix_[64*64]);
void GenerateHeavyHashMatrix(uint256 matrix_seed, uint64_t matrix[64*64]);
void serialize_heavyhash(Obtc_t *Obtc, uint64_t matrix[64*64], const char* in, char* out, int len);
void opticalbtc_hash(const char* in, char* out, int len);
extern void CSHA3_256_init(Obtc_t *Obtc, CSHA3_256 *p);
extern void CSHA3_256_CSHA3_256(Obtc_t *Obtc, CSHA3_256 *p);
extern void CSHA3_256_Write(CSHA3_256 *p, const unsigned char* data, size_t len);
extern void CSHA3_256_Finalize(CSHA3_256 *p, unsigned char hash[OUTPUT_SIZE]);
//extern void CSHA3_256_Reset(Obtc_t *Obtc, CSHA3_256 *p);
extern void CHeavyHash_init(Obtc_t *Obtc, CHeavyHash *p, uint64_t matrix_[64*64]);
extern void CHeavyHash_Write(CHeavyHash *p, const unsigned char* data, size_t len);
extern void CHeavyHash_Finalize(Obtc_t *Obtc, CHeavyHash *p, unsigned char hash[OUTPUT_SIZE]);
extern int sha3_init(Obtc_t *Obtc,sha3_ctx_t *c, int mdlen); // mdlen = hash output in bytes
#endif // OBTC_H

BIN
heavyHash/obtc.o Normal file
View File

Binary file not shown.

199
heavyHash/sha3.c Normal file
View File

@@ -0,0 +1,199 @@
// sha3.c
// 19-Nov-11 Markku-Juhani O. Saarinen <mjos@iki.fi>
// Revised 07-Aug-15 to match with official release of FIPS PUB 202 "SHA3"
// Revised 03-Sep-15 for portability + OpenSSL - style API
#include <stdio.h>
#include "sha3.h"
#include "obtc.h"
// update the state with given number of rounds
void sha3_keccakf(uint64_t st[25])
{
// constants
const uint64_t keccakf_rndc[24] = {
0x0000000000000001, 0x0000000000008082, 0x800000000000808a,
0x8000000080008000, 0x000000000000808b, 0x0000000080000001,
0x8000000080008081, 0x8000000000008009, 0x000000000000008a,
0x0000000000000088, 0x0000000080008009, 0x000000008000000a,
0x000000008000808b, 0x800000000000008b, 0x8000000000008089,
0x8000000000008003, 0x8000000000008002, 0x8000000000000080,
0x000000000000800a, 0x800000008000000a, 0x8000000080008081,
0x8000000000008080, 0x0000000080000001, 0x8000000080008008
};
const int keccakf_rotc[24] = {
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 2, 14,
27, 41, 56, 8, 25, 43, 62, 18, 39, 61, 20, 44
};
const int keccakf_piln[24] = {
10, 7, 11, 17, 18, 3, 5, 16, 8, 21, 24, 4,
15, 23, 19, 13, 12, 2, 20, 14, 22, 9, 6, 1
};
// variables
int i, j, r;
uint64_t t, bc[5];
#if __BYTE_ORDER__ != __ORDER_LITTLE_ENDIAN__
uint8_t *v;
// endianess conversion. this is redundant on little-endian targets
for (i = 0; i < 25; i++) {
v = (uint8_t *) &st[i];
st[i] = ((uint64_t) v[0]) | (((uint64_t) v[1]) << 8) |
(((uint64_t) v[2]) << 16) | (((uint64_t) v[3]) << 24) |
(((uint64_t) v[4]) << 32) | (((uint64_t) v[5]) << 40) |
(((uint64_t) v[6]) << 48) | (((uint64_t) v[7]) << 56);
}
#endif
// actual iteration
for (r = 0; r < KECCAKF_ROUNDS; r++) {
// Theta
for (i = 0; i < 5; i++)
bc[i] = st[i] ^ st[i + 5] ^ st[i + 10] ^ st[i + 15] ^ st[i + 20];
for (i = 0; i < 5; i++) {
t = bc[(i + 4) % 5] ^ ROTL64(bc[(i + 1) % 5], 1);
for (j = 0; j < 25; j += 5)
st[j + i] ^= t;
}
// Rho Pi
t = st[1];
for (i = 0; i < 24; i++) {
j = keccakf_piln[i];
bc[0] = st[j];
st[j] = ROTL64(t, keccakf_rotc[i]);
t = bc[0];
}
// Chi
for (j = 0; j < 25; j += 5) {
for (i = 0; i < 5; i++)
bc[i] = st[j + i];
for (i = 0; i < 5; i++)
st[j + i] ^= (~bc[(i + 1) % 5]) & bc[(i + 2) % 5];
}
// Iota
st[0] ^= keccakf_rndc[r];
}
#if __BYTE_ORDER__ != __ORDER_LITTLE_ENDIAN__
// endianess conversion. this is redundant on little-endian targets
for (i = 0; i < 25; i++) {
v = (uint8_t *) &st[i];
t = st[i];
v[0] = t & 0xFF;
v[1] = (t >> 8) & 0xFF;
v[2] = (t >> 16) & 0xFF;
v[3] = (t >> 24) & 0xFF;
v[4] = (t >> 32) & 0xFF;
v[5] = (t >> 40) & 0xFF;
v[6] = (t >> 48) & 0xFF;
v[7] = (t >> 56) & 0xFF;
}
#endif
}
// Initialize the context for SHA3
int sha3_init(Obtc_t *Obtc, sha3_ctx_t *c, int mdlen)
{
int i;
for (i = 0; i < 200; i++){
c->st.b[i] = Obtc->const_data[199-i];
}
c->mdlen = mdlen;
c->rsiz = 200 - 2 * mdlen;
c->pt = 0;
return 1;
}
// update state with more data
int sha3_update(sha3_ctx_t *c, const void *data, size_t len)
{
size_t i;
int j;
j = c->pt;
for (i = 0; i < len; i++) {
c->st.b[j++] ^= ((const uint8_t *) data)[i];
if (j >= c->rsiz) {
sha3_keccakf(c->st.q);
j = 0;
}
}
c->pt = j;
return 1;
}
// finalize and output a hash
int sha3_final(void *md, sha3_ctx_t *c)
{
int i;
// c->st.b[c->pt] ^= 0x06;
c->st.b[c->pt] ^= 0x04;
c->st.b[c->rsiz - 1] ^= 0x80;
sha3_keccakf(c->st.q);
for (i = 0; i < c->mdlen; i++) {
((uint8_t *) md)[i] = c->st.b[i];
}
return 1;
}
// compute a SHA-3 hash (md) of given byte length from "in"
/*void *sha3(const void *in, size_t inlen, void *md, int mdlen)
{
sha3_ctx_t sha3;
sha3_init(&sha3, mdlen);
sha3_update(&sha3, in, inlen);
sha3_final(md, &sha3);
return md;
}*/
// SHAKE128 and SHAKE256 extensible-output functionality
void shake_xof(sha3_ctx_t *c)
{
c->st.b[c->pt] ^= 0x1F;
c->st.b[c->rsiz - 1] ^= 0x80;
sha3_keccakf(c->st.q);
c->pt = 0;
}
void shake_out(sha3_ctx_t *c, void *out, size_t len)
{
size_t i;
int j;
j = c->pt;
for (i = 0; i < len; i++) {
if (j >= c->rsiz) {
sha3_keccakf(c->st.q);
j = 0;
}
((uint8_t *) out)[i] = c->st.b[j++];
}
c->pt = j;
}

51
heavyHash/sha3.h Normal file
View File

@@ -0,0 +1,51 @@
// sha3.h
// 19-Nov-11 Markku-Juhani O. Saarinen <mjos@iki.fi>
#ifndef SHA3_H
#define SHA3_H
#include <stddef.h>
#include <stdint.h>
#ifndef KECCAKF_ROUNDS
#define KECCAKF_ROUNDS 24
#endif
#ifndef ROTL64
#define ROTL64(x, y) (((x) << (y)) | ((x) >> (64 - (y))))
#endif
// state context
typedef struct {
union { // state:
uint8_t b[200]; // 8-bit bytes
uint64_t q[25]; // 64-bit words
} st;
int pt, rsiz, mdlen; // these don't overflow
} sha3_ctx_t;
// Compression function.
void sha3_keccakf(uint64_t st[25]);
// OpenSSL - like interfece
int sha3_update(sha3_ctx_t *c, const void *data, size_t len);
int sha3_final(void *md, sha3_ctx_t *c); // digest goes to md
// compute a sha3 hash (md) of given byte length from "in"
void *sha3(const void *in, size_t inlen, void *md, int mdlen);
// SHAKE128 and SHAKE256 extensible-output functions
//#define shake128_init(c) sha3_init(c, 16)
//#define shake256_init(c) sha3_init(c, 32)
//#define shake_update sha3_update
void shake_xof(sha3_ctx_t *c);
void shake_out(sha3_ctx_t *c, void *out, size_t len);
#endif

BIN
heavyHash/sha3.o Normal file
View File

Binary file not shown.

42
heavyHash/singular.h Normal file
View File

@@ -0,0 +1,42 @@
#ifndef _SINGULAR_SINGULAR_H
#define _SINGULAR_SINGULAR_H
/** The version of the singular library. */
#define SINGULAR_VERSION "@PROJECT_VERSION@"
/**
* Whether rvalue references are supported.
*
* Visual Studio 2010 and lower do not have rvalue references so far.
*/
#if defined(_MSC_VER) && _MSC_VER < 1700
#define SINGULAR_RVALUE_REFERENCE_SUPPORTED 0
#else
#define SINGULAR_RVALUE_REFERENCE_SUPPORTED 1
#endif
/**
* Whether function deletions are supported.
*
* Visual Studio 2012 and lower do not like "delete" stuff so far.
*/
#if defined(_MSC_VER) && _MSC_VER < 1800
#define SINGULAR_FUNCTION_DELETION_SUPPORTED 0
#else
#define SINGULAR_FUNCTION_DELETION_SUPPORTED 1
#endif
/**
* Whether template friend operator overalodings are supported.
*
* Visual Studio 2012 and lower do not like overloading a template firend
* operators.
* Neither does GCC.
*/
#if (defined(_MSC_VER) && _MSC_VER < 1800) || (defined(__GNUC__) && !defined(__clang__))
#define SINGULAR_TEMPLATE_FRIEND_OPERATOR_OVERLOADING_SUPPORTED 0
#else
#define SINGULAR_TEMPLATE_FRIEND_OPERATOR_OVERLOADING_SUPPORTED 1
#endif
#endif

183
heavyHash/test.c Normal file
View File

@@ -0,0 +1,183 @@
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include "obtc.h"
#include "singular.h"
#include<time.h>
//uint8_t const_data[200];
static const int hex2bin_tbl[256] = {
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, -1, -1, -1, -1, -1, -1,
-1, 10, 11, 12, 13, 14, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, 10, 11, 12, 13, 14, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
};
bool hex2bin(unsigned char *p, const char *hexstr, size_t len)
{
int nibble1, nibble2;
unsigned char idx;
bool ret = false;
while (*hexstr && len) {
if ((!hexstr[1])) {
printf("hex2bin str truncated");
return ret;
}
idx = *hexstr++;
nibble1 = hex2bin_tbl[idx];
idx = *hexstr++;
nibble2 = hex2bin_tbl[idx];
if (((nibble1 < 0) || (nibble2 < 0))) {
printf("hex2bin scan failed");
return ret;
}
*p++ = (((unsigned char)nibble1) << 4) | ((unsigned char)nibble2);
--len;
}
if ((len == 0 && *hexstr == 0))
ret = true;
return ret;
}
int main(int argc, char **argv)
{
uint8_t genesis_block[80];
uint8_t hash[32];
uint8_t last_prehash[32];
uint8_t last_prehash2[32];
uint8_t prehash_tab[32];
uint8_t nonce_tab[8];
char *prehash_str = "d76ffb1d8e31ec04579b0452b52bde7dbd088e912ab1b11ba924ff309ab44a43";//argv[1];
char *nonce_str = "80aa59a7901f2502";//argv[2];
//char *last_prehash_str = argv[3];
//char *last_prehash_str2 = argv[4];
hex2bin(prehash_tab, prehash_str, strlen(prehash_str)/2);
hex2bin(nonce_tab, nonce_str, strlen(nonce_str)/2);
//hex2bin(last_prehash, last_prehash_str, strlen(last_prehash_str)/2);
//hex2bin(last_prehash2, last_prehash_str2, strlen(last_prehash_str2)/2);
/*for (uint8_t i = 0; i<32;i++){
printf("0x%x, ",prehash_tab[i]);
}
printf("\n");
for (uint8_t i = 0; i<8;i++){
printf("0x%x, ",nonce_tab[i]);
}
printf("\n");*/
//uint8_t prehash[32] = {0x81,0x55,0x3a,0x69,0x5a,0x05,0x88,0x99,0x8c,0x41,0x37,0x92,0xe7,0x4c,0xe8,0xb8,0xf8,0xa0,0x96,0xd6,0x4b,0x3e,0xe4,0x73,0x87,0x37,0x24,0x34,0x48,0x5c,0x0b,0x6f};
//uint8_t utime[8] = {0x00,0x00,0x01,0x84,0x8c,0xa8,0x7c,0x49};
uint8_t pad[32] = {0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00};
//uint8_t nonce[8] = {0x2f,0x84,0x00,0x00,0x0e,0xba,0x16,0x7c};
#if 0
//uint8_t prehash[32] = {0xa4,0x8f,0xae,0x69,0xeb,0x28,0xc7,0xe0,0x14,0x11,0x4f,0x01,0xae,0x60,0xc8,0xc3,0x82,0x73,0xc4,0x60,0x66,0xcf,0x95,0xd6,0x77,0x1a,0x55,0xd6,0x16,0xd7,0xa1,0x9a};//大端
//uint8_t utime[8] = {0x00,0x00,0x01,0x87,0x22,0x1e,0xad,0x44};
//uint8_t nonce[8] = {0x8e,0xd4,0x00,0x10,0x6b,0xe7,0xe4,0x00};
//uint8_t nonce[8] = {0x8e,0xd4,0x00,0x12,0x27,0xc6,0x90,0xa0};
//uint8_t nonce[8] = {0x8e,0xd4,0x00,0x32,0x0b,0x6b,0xd6,0xd1};
//3f 9a aa c6 32 af 1a 4e 0e 1f ea 8a f8 e3 d5 32 b7 5a a4 71 b2 e4 ef fe a5 bd cc fa 3b dd b6 61
uint8_t prehash[32] = {0x3f,0x9a,0xaa,0xc6,0x32,0xaf,0x1a,0x4e,0x0e,0x1f,0xea,0x8a,0xf8,0xe3,0xd5,0x32,0xb7,0x5a,0xa4,0x71,0xb2,0xe4,0xef,0xfe,0xa5,0xbd,0xcc,0xfa,0x3b,0xdd,0xb6,0x61};//大端
uint8_t utime[8] = {0x00,0x00,0x01,0x87,0x21,0xeb,0x73,0x79};
uint8_t nonce[8] = {0xa3,0xdd,0x02,0x10,0x1a,0x87,0xb4,0x70};
#else
/*443e01000000ffff00000000
e0af2a3ba173157d3f70c94aad742fdf16d9930fdfc9d6301e869bcef04ced6c
e0af2a3ba173157d3f70c94aad742fdf16d9930fdfc9d6301e869bcef04ced6c
dbee84288701000000000000901f25020000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
[2023-03-28 22:00:46.549] 00 cc 01 11 70 83 85 16 90 1f 25 02
kas_pow_hash: in:e0af2a3ba173157d3f70c94aad742fdf16d9930fdfc9d6301e869bcef04ced6cdbee842887010000000000000000000000000000000000000000000000000000000000000000000070838516901f2502
kas_pow_hash: out:dae78f5008d3b66f
01a740ce33c812ba
772b3f5763da7bc6
da24cb6c00000000*/
uint8_t prehash[32] = {0xe0,0xaf,0x2a,0x3b,0xa1,0x73,0x15,0x7d,0x3f,0x70,0xc9,0x4a,0xad,0x74,0x2f,0xdf,0x16,0xd9,0x93,0x0f,0xdf,0xc9,0xd6,0x30,0x1e,0x86,0x9b,0xce,0xf0,0x4c,0xed,0x6c};
//uint8_t utime[8] = {0x00,0x00,0x01,0x87,0x28,0x84,0xee,0xdb};
uint8_t nonce[8] = {0x02,0x25,0x1f,0x90,0x16,0x85,0x83,0x70};
#endif
/*for (int i = 0; i < 32; ++i) genesis_block[i] = prehash[i];
for (int i = 0; i < 8; ++i) genesis_block[i+32] = utime[7-i];
for (int i = 0; i < 32; ++i) genesis_block[i+40] = pad[31-i];
for (int i = 0; i < 8; ++i) genesis_block[i+72] = nonce[7-i];*/
//uint8_t utime[8] = {0x00,0x00,0x01,0x87,0x21,0xeb,0x73,0x79};
//dbee8428870100000
uint8_t utime[8] = {0x00,0x00,0x01,0x87,0x28,0x84,0xee,0xdb};
for (int i = 0; i < 32; ++i) genesis_block[i] = prehash_tab[i];
for (int i = 0; i < 8; ++i) genesis_block[i+32] = utime[7-i];
for (int i = 0; i < 32; ++i) genesis_block[i+40] = pad[31-i];
for (int i = 0; i < 8; ++i) genesis_block[i+72] = nonce_tab[i];
clock_t start, finish;
double Total_time;
uint32_t cnt = 0;;
//while(1)
{
start = clock();
opticalbtc_hash((const char*)&genesis_block, (char*)&hash, sizeof(genesis_block));
finish = clock();
Total_time = (double)(finish-start) / CLOCKS_PER_SEC;
printf( "\n cnt = %d, opticalbtc_hash run times %f seconds\n", cnt++, Total_time);
for (int i=31; i>-1; i--) {
printf("%02hhx", hash[i]);
}
printf("\n");
}
//if (hash[31] != 0 || hash[30] != 0){
// for (int i = 0; i < 32; ++i) genesis_block[i] = last_prehash[i];
// opticalbtc_hash((const char*)&genesis_block, (char*)&hash, sizeof(genesis_block));
//}
//if (hash[31] != 0 || hash[30] != 0){
// for (int i = 0; i < 32; ++i) genesis_block[i] = last_prehash2[i];
// opticalbtc_hash((const char*)&genesis_block, (char*)&hash, sizeof(genesis_block));
//}
if (hash[31] != 0 && hash[30] != 0){
printf("reject\n");
}
return 0;
}
//g++ -std=c++11 *.cpp

44
heavyHash/uint256.h Normal file
View File

@@ -0,0 +1,44 @@
// Copyright (c) 2009-2010 Satoshi Nakamoto
// Copyright (c) 2009-2016 The Bitcoin Core developers
// Distributed under the MIT software license, see the accompanying
// file COPYING or http://www.opensource.org/licenses/mit-license.php.
#ifndef BITCOIN_UINT256_H
#define BITCOIN_UINT256_H
#include <assert.h>
//#include <cstring>
//#include <stdexcept>
#include <stdint.h>
#include <string.h>
//#include <vector>
#include <stdbool.h>
/** 256-bit opaque blob.
* @note This type is called uint256 for historical reasons only. It is an
* opaque blob of 256 bits and has no integer operations. Use arith_uint256 if
* those are required.
*/
#define UPPER_P(x) x->elements[0]
#define LOWER_P(x) x->elements[1]
#define UPPER(x) x.elements[0]
#define LOWER(x) x.elements[1]
#define WIDTH 32
typedef struct class_base_blob base_blob_t;
struct class_base_blob{
uint8_t data[WIDTH];
};
typedef struct uint128_t { uint64_t elements[2]; } uint128_t;
typedef struct uint256_t {
uint128_t elements[2];
base_blob_t bb;
} uint256;
#endif // BITCOIN_UINT256_H

15
heavyHash/xoshiro256pp.h Normal file
View File

@@ -0,0 +1,15 @@
#ifndef OPOW_CRYPTO_XOSHIRO256PP_H
#define OPOW_CRYPTO_XOSHIRO256PP_H
#include <stdint.h>
#include "uint256.h"
typedef struct class_XoShiRo256PlusPlus XoShiRo256PlusPlus_t;
struct class_XoShiRo256PlusPlus{
uint64_t s[4];
};
#endif //OPOW_CRYPTO_XOSHIRO256PP_H