398 lines
11 KiB
C
398 lines
11 KiB
C
/*******************************************************************
|
|
*
|
|
* M4RI: Linear Algebra over GF(2)
|
|
*
|
|
* Copyright (C) 2011 Martin Albrecht <martinralbrecht@googlemail.com>
|
|
*
|
|
* Distributed under the terms of the GNU General Public License (GPL)
|
|
* version 2 or higher.
|
|
*
|
|
* This code is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
* General Public License for more details.
|
|
*
|
|
* The full text of the GPL is available at:
|
|
*
|
|
* http://www.gnu.org/licenses/
|
|
*
|
|
********************************************************************/
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include "m4ri_config.h"
|
|
|
|
#if __M4RI_HAVE_LIBPNG
|
|
#include <png.h>
|
|
#endif //__M4RI_HAVE_LIBPNG
|
|
|
|
|
|
#include "io.h"
|
|
#include "echelonform.h"
|
|
|
|
void mzd_info(const mzd_t *A, int do_rank) {
|
|
printf("nrows: %6zu, ncols: %6zu, density: %6.5f, hash: 0x%016zx",(size_t)A->nrows,(size_t)A->ncols,mzd_density(A,1),mzd_hash(A));
|
|
if(do_rank) {
|
|
mzd_t *AA = mzd_copy(NULL, A);
|
|
printf(", rank: %6zu\n",(size_t)mzd_echelonize(AA,0));
|
|
mzd_free(AA);
|
|
} else {
|
|
printf("\n");
|
|
}
|
|
}
|
|
|
|
#define SAFECHAR (m4ri_radix + m4ri_radix / 4 + 1)
|
|
|
|
void mzd_print(mzd_t const *M) {
|
|
char temp[SAFECHAR];
|
|
for (rci_t i = 0; i < M->nrows; ++i) {
|
|
printf("[");
|
|
word *row = M->rows[i];
|
|
for (wi_t j = 0; j < M->width - 1; ++j) {
|
|
m4ri_word_to_str(temp, row[j], 1);
|
|
printf("%s|", temp);
|
|
}
|
|
row = row + M->width - 1;
|
|
int const wide = (M->ncols % m4ri_radix) ? M->ncols % m4ri_radix : m4ri_radix;
|
|
for (int j = 0; j < wide; ++j) {
|
|
if(j != 0 && (j % 4) == 0)
|
|
printf(":");
|
|
if (__M4RI_GET_BIT(*row, j))
|
|
printf("1");
|
|
else
|
|
printf(" ");
|
|
}
|
|
printf("]\n");
|
|
}
|
|
}
|
|
|
|
void mzd_print_row(mzd_t const *M, const rci_t i) {
|
|
char temp[SAFECHAR];
|
|
printf("[");
|
|
word *row = M->rows[i];
|
|
for (wi_t j = 0; j < M->width - 1; ++j) {
|
|
m4ri_word_to_str(temp, row[j], 1);
|
|
printf("%s|", temp);
|
|
}
|
|
row = row + M->width - 1;
|
|
int const wide = (M->ncols % m4ri_radix) ? M->ncols % m4ri_radix : m4ri_radix;
|
|
for (int j = 0; j < wide; ++j) {
|
|
if(j != 0 && (j % 4) == 0)
|
|
printf(":");
|
|
if (__M4RI_GET_BIT(*row, j))
|
|
printf("1");
|
|
else
|
|
printf(" ");
|
|
}
|
|
printf("]\n");
|
|
}
|
|
|
|
|
|
#if __M4RI_HAVE_LIBPNG
|
|
#define PNGSIGSIZE 8
|
|
|
|
mzd_t * mzd_from_png(const char *fn, int verbose) {
|
|
int retval = 0;
|
|
mzd_t *A = NULL;
|
|
png_byte pngsig[PNGSIGSIZE];
|
|
|
|
FILE *fh = fopen(fn,"rb");
|
|
|
|
if (!fh) {
|
|
if (verbose)
|
|
printf("Could not open file '%s' for reading\n",fn);
|
|
return NULL;
|
|
};
|
|
|
|
if (fread((char*)pngsig, PNGSIGSIZE, 1, fh) != 1) {
|
|
if (verbose)
|
|
printf("Could not read file '%s'\n",fn);
|
|
retval = 1;
|
|
goto from_png_close_fh;
|
|
}
|
|
|
|
if (png_sig_cmp(pngsig, 0, PNGSIGSIZE) != 0) {
|
|
if (verbose)
|
|
printf("'%s' is not a PNG file.\n",fn);
|
|
retval = 2;
|
|
goto from_png_close_fh;
|
|
}
|
|
|
|
png_structp png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
|
|
|
|
if (!png_ptr) {
|
|
if (verbose)
|
|
printf("failed to initialise PNG read struct.\n");
|
|
retval = 3;
|
|
goto from_png_close_fh;
|
|
}
|
|
png_set_user_limits(png_ptr, 0x7fffffffL, 0x7fffffffL);
|
|
|
|
png_infop info_ptr = png_create_info_struct(png_ptr);
|
|
|
|
if (!info_ptr) {
|
|
if (verbose)
|
|
printf("failed to initialise PNG info struct\n");
|
|
retval = 3;
|
|
goto from_png_destroy_read_struct;
|
|
}
|
|
|
|
png_init_io(png_ptr, fh);
|
|
png_set_sig_bytes(png_ptr, PNGSIGSIZE);
|
|
png_read_info(png_ptr, info_ptr);
|
|
|
|
const png_uint_32 m = png_get_image_height(png_ptr, info_ptr);
|
|
const png_uint_32 n = png_get_image_width(png_ptr, info_ptr);
|
|
const png_uint_32 bit_depth = png_get_bit_depth(png_ptr, info_ptr);
|
|
const png_uint_32 channels = png_get_channels(png_ptr, info_ptr);
|
|
const png_uint_32 color_type = png_get_color_type(png_ptr, info_ptr);
|
|
const png_uint_32 compression_type = png_get_compression_type(png_ptr, info_ptr);
|
|
const png_uint_32 interlace_type = png_get_interlace_type(png_ptr, info_ptr);
|
|
|
|
if (interlace_type != PNG_INTERLACE_NONE) {
|
|
if (verbose)
|
|
printf("interlaced images not supported\n");
|
|
goto from_png_destroy_read_struct;
|
|
};
|
|
|
|
if (verbose)
|
|
printf("reading %u x %u matrix (bit depth: %u, channels: %u, color type: %u, compression type: %u)\n",(unsigned int)m,
|
|
(unsigned int)n, (unsigned int)bit_depth, (unsigned int)channels, (unsigned int)color_type, (unsigned int)compression_type);
|
|
|
|
if(color_type != 0 && color_type != 3) {
|
|
if (verbose)
|
|
printf("only graycscale and palette colors are supported.\n");
|
|
goto from_png_destroy_read_struct;
|
|
}
|
|
|
|
A = mzd_init(m, n);
|
|
const word bitmask_end = A->high_bitmask;
|
|
png_bytep row = m4ri_mm_calloc(sizeof(char),n/8+1);
|
|
|
|
word tmp;
|
|
wi_t j;
|
|
|
|
png_set_packswap(png_ptr);
|
|
//png_set_invert_mono(png_ptr);
|
|
|
|
for(rci_t i=0; i<m; i++) {
|
|
png_read_row(png_ptr, row, NULL);
|
|
for(j=0; j<A->width-1; j++) {
|
|
tmp = ((word)row[8*j+7])<<56 | ((word)row[8*j+6])<<48 \
|
|
| ((word)row[8*j+5])<<40 | ((word)row[8*j+4])<<32 \
|
|
| ((word)row[8*j+3])<<24 | ((word)row[8*j+2])<<16 \
|
|
| ((word)row[8*j+1])<< 8 | ((word)row[8*j+0])<< 0;
|
|
A->rows[i][j] = ~tmp;
|
|
}
|
|
tmp = 0;
|
|
switch((n/8 + ((n%8) ? 1 : 0))%8) {
|
|
case 0: tmp |= ((word)row[8*j+7])<<56;
|
|
case 7: tmp |= ((word)row[8*j+6])<<48;
|
|
case 6: tmp |= ((word)row[8*j+5])<<40;
|
|
case 5: tmp |= ((word)row[8*j+4])<<32;
|
|
case 4: tmp |= ((word)row[8*j+3])<<24;
|
|
case 3: tmp |= ((word)row[8*j+2])<<16;
|
|
case 2: tmp |= ((word)row[8*j+1])<< 8;
|
|
case 1: tmp |= ((word)row[8*j+0])<< 0;
|
|
};
|
|
A->rows[i][j] |= (~tmp & bitmask_end);
|
|
}
|
|
|
|
m4ri_mm_free(row);
|
|
png_read_end(png_ptr, NULL);
|
|
|
|
from_png_destroy_read_struct:
|
|
png_destroy_read_struct(&png_ptr, &info_ptr,(png_infopp)0);
|
|
|
|
from_png_close_fh:
|
|
fclose(fh);
|
|
|
|
if (retval != 0 && A) {
|
|
mzd_free(A);
|
|
return NULL;
|
|
} else {
|
|
return A;
|
|
}
|
|
}
|
|
|
|
int mzd_to_png(const mzd_t *A, const char *fn, int compression_level, const char *comment, int verbose) {
|
|
FILE *fh = fopen(fn, "wb");
|
|
|
|
if (!fh) {
|
|
if(verbose)
|
|
printf("Could not open file '%s' for writing\n",fn);
|
|
return 1;
|
|
}
|
|
|
|
png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
|
|
|
|
if (!png_ptr) {
|
|
if(verbose)
|
|
printf("failed to initialise PNG write struct.\n");
|
|
fclose(fh);
|
|
return 3;
|
|
}
|
|
png_set_user_limits(png_ptr, 0x7fffffffL, 0x7fffffffL);
|
|
|
|
png_infop info_ptr = png_create_info_struct(png_ptr);
|
|
|
|
if (!info_ptr) {
|
|
if (verbose)
|
|
printf("failed to initialise PNG info struct\n");
|
|
png_destroy_write_struct(&png_ptr, &info_ptr);
|
|
fclose(fh);
|
|
return 3;
|
|
}
|
|
|
|
if (setjmp(png_jmpbuf(png_ptr))) {
|
|
if (verbose)
|
|
printf("error writing PNG file\n");
|
|
png_destroy_write_struct(&png_ptr, &info_ptr);
|
|
fclose(fh);
|
|
return 1;
|
|
}
|
|
|
|
png_init_io(png_ptr, fh);
|
|
png_set_compression_level(png_ptr, compression_level);
|
|
|
|
png_set_IHDR(png_ptr, info_ptr, A->ncols, A->nrows, 1, \
|
|
PNG_COLOR_TYPE_GRAY, \
|
|
PNG_INTERLACE_NONE, \
|
|
PNG_COMPRESSION_TYPE_DEFAULT, \
|
|
PNG_FILTER_TYPE_DEFAULT);
|
|
|
|
png_text txt_ptr[3];
|
|
|
|
char pdate[21];
|
|
time_t ptime=time(NULL);
|
|
struct tm *ltime=localtime(&ptime);
|
|
sprintf(pdate,"%04d/%02d/%02d %02d:%02d:%02d",ltime->tm_year+1900,ltime->tm_mon+1,ltime->tm_mday,ltime->tm_hour,ltime->tm_min,ltime->tm_sec);
|
|
|
|
txt_ptr[0].key="Software";
|
|
txt_ptr[0].text="M4RI";
|
|
txt_ptr[0].compression=PNG_TEXT_COMPRESSION_NONE;
|
|
txt_ptr[1].key="Date";
|
|
txt_ptr[1].text=pdate;
|
|
txt_ptr[1].compression=PNG_TEXT_COMPRESSION_NONE;
|
|
txt_ptr[2].key="Comment";
|
|
txt_ptr[2].text=(char*)comment;
|
|
txt_ptr[2].compression=PNG_TEXT_COMPRESSION_NONE;
|
|
|
|
png_set_text(png_ptr, info_ptr, txt_ptr, 3);
|
|
|
|
png_write_info(png_ptr, info_ptr);
|
|
|
|
png_set_packswap(png_ptr);
|
|
png_set_invert_mono(png_ptr);
|
|
|
|
png_bytep row = m4ri_mm_calloc(sizeof(char),A->ncols/8+8);
|
|
|
|
wi_t j=0;
|
|
word tmp = 0;
|
|
for(rci_t i=0; i<A->nrows; i++) {
|
|
word *rowptr = A->rows[i];
|
|
for(j=0; j<A->width-1; j++) {
|
|
tmp = rowptr[j];
|
|
row[8*j+0] = (png_byte)((tmp>> 0) & 0xff);
|
|
row[8*j+1] = (png_byte)((tmp>> 8) & 0xff);
|
|
row[8*j+2] = (png_byte)((tmp>>16) & 0xff);
|
|
row[8*j+3] = (png_byte)((tmp>>24) & 0xff);
|
|
row[8*j+4] = (png_byte)((tmp>>32) & 0xff);
|
|
row[8*j+5] = (png_byte)((tmp>>40) & 0xff);
|
|
row[8*j+6] = (png_byte)((tmp>>48) & 0xff);
|
|
row[8*j+7] = (png_byte)((tmp>>56) & 0xff);
|
|
}
|
|
tmp = rowptr[j];
|
|
switch( (A->ncols/8 + ((A->ncols%8) ? 1 : 0)) %8 ) {
|
|
case 0: row[8*j+7] = (png_byte)((tmp>>56) & 0xff);
|
|
case 7: row[8*j+6] = (png_byte)((tmp>>48) & 0xff);
|
|
case 6: row[8*j+5] = (png_byte)((tmp>>40) & 0xff);
|
|
case 5: row[8*j+4] = (png_byte)((tmp>>32) & 0xff);
|
|
case 4: row[8*j+3] = (png_byte)((tmp>>24) & 0xff);
|
|
case 3: row[8*j+2] = (png_byte)((tmp>>16) & 0xff);
|
|
case 2: row[8*j+1] = (png_byte)((tmp>> 8) & 0xff);
|
|
case 1: row[8*j+0] = (png_byte)((tmp>> 0) & 0xff);
|
|
};
|
|
png_write_row(png_ptr, row);
|
|
}
|
|
m4ri_mm_free(row);
|
|
|
|
png_write_end(png_ptr, info_ptr);
|
|
png_destroy_write_struct(&png_ptr, &info_ptr);
|
|
fclose(fh);
|
|
return 0;
|
|
}
|
|
|
|
#endif //__M4RI_HAVE_LIBPNG
|
|
|
|
mzd_t *mzd_from_jcf(const char *fn, int verbose) {
|
|
int retval = 0;
|
|
mzd_t *A = NULL;
|
|
FILE *fh = fopen(fn,"r");
|
|
|
|
rci_t m,n;
|
|
long p = 0;
|
|
long nonzero = 0;
|
|
|
|
if (!fh) {
|
|
if (verbose)
|
|
printf("Could not open file '%s' for reading\n",fn);
|
|
return NULL;
|
|
}
|
|
|
|
if (fscanf(fh, "%d %d %ld\n%ld\n\n",&m,&n,&p,&nonzero) != 4) {
|
|
if (verbose)
|
|
printf("File '%s' does not seem to be in JCF format.",fn);
|
|
retval = 1;
|
|
goto from_jcf_close_fh;
|
|
}
|
|
|
|
if(p != 2) {
|
|
if (verbose)
|
|
printf("Expected p==2 but found p==%ld\n",p);
|
|
retval = 1;
|
|
goto from_jcf_close_fh;
|
|
}
|
|
|
|
if (verbose)
|
|
printf("reading %lu x %lu matrix with at most %ld non-zero entries (density at most: %6.5f)\n",
|
|
(unsigned long)m, (unsigned long)n, (unsigned long)nonzero, ((double)nonzero)/((double)m*n));
|
|
|
|
A = mzd_init(m,n);
|
|
|
|
long i = -1;
|
|
long j = 0;
|
|
|
|
while(fscanf(fh,"%ld\n",&j) == 1) {
|
|
if (j<0) {
|
|
i++, j = -j;
|
|
}
|
|
if (((j-1) >= n) || (i>= m))
|
|
m4ri_die("trying to write to (%ld,%ld) in %ld x %ld matrix\n", i, j-1, m, n);
|
|
mzd_write_bit(A, i, j-1, 1);
|
|
};
|
|
|
|
from_jcf_close_fh:
|
|
fclose(fh);
|
|
|
|
if(retval != 0 && A) {
|
|
mzd_free(A);
|
|
return NULL;
|
|
} else {
|
|
return A;
|
|
}
|
|
}
|
|
|
|
mzd_t *mzd_from_str(rci_t m, rci_t n, const char *str) {
|
|
int idx = 0;
|
|
mzd_t *A = mzd_init(m, n);
|
|
for(rci_t i=0; i<A->nrows; i++) {
|
|
for(rci_t j=0; j<A->ncols; j++) {
|
|
mzd_write_bit(A, i, j, str[idx++] == '1');
|
|
}
|
|
}
|
|
return A;
|
|
}
|