diff --git a/.github/workflows/msolve.yml b/.github/workflows/msolve.yml index e3804686..cfc381fb 100644 --- a/.github/workflows/msolve.yml +++ b/.github/workflows/msolve.yml @@ -9,15 +9,30 @@ on: jobs: build: runs-on: ${{ matrix.os }} - timeout-minutes: 10 + timeout-minutes: 20 strategy: fail-fast: false matrix: - os: - - ubuntu-latest - - macos-latest + include: + - { os: ubuntu-latest, shell: bash } + - { os: macos-latest, shell: bash } + - { os: windows-latest, shell: msys2 } + defaults: + run: + shell: ${{ matrix.shell }} {0} steps: + - name: Disable CRLF conversion on Windows + if: runner.os == 'Windows' + shell: pwsh + run: | + git config --global core.autocrlf false + git config --global core.eol lf - uses: actions/checkout@v3 + - name: Set up MSYS2 + if: runner.os == 'Windows' + uses: msys2/setup-msys2@v2 + with: + update: true - name: "Install dependencies" run: | if [ "$RUNNER_OS" == "Linux" ]; then @@ -25,6 +40,8 @@ jobs: sudo apt install libgmp-dev libflint-dev libmpfr-dev libntl-dev elif [ "$RUNNER_OS" == "macOS" ]; then brew install autoconf automake libtool gmp flint mpfr ntl + elif [ "$RUNNER_OS" == "Windows" ]; then + pacman -S --noconfirm autotools mingw-w64-x86_64-toolchain mingw-w64-x86_64-gmp mingw-w64-x86_64-flint mingw-w64-x86_64-mpfr mingw-w64-x86_64-ntl else echo "$RUNNER_OS not supported" exit 1 @@ -33,7 +50,7 @@ jobs: run: ./autogen.sh - name: configure run: | - if [ "$RUNNER_OS" == "Linux" ]; then + if [ "$RUNNER_OS" = "Linux" ] || [ "$RUNNER_OS" = "Windows" ]; then ./configure elif [ "$RUNNER_OS" == "macOS" ]; then ./configure LDFLAGS="-L$(brew --prefix)/lib/" CFLAGS="-g -O2 -I$(brew --prefix)/include/" @@ -47,7 +64,7 @@ jobs: run: make check - name: make distcheck run: | - if [ "$RUNNER_OS" == "Linux" ]; then + if [ "$RUNNER_OS" = "Linux" ] || [ "$RUNNER_OS" = "Windows" ]; then make distcheck elif [ "$RUNNER_OS" == "macOS" ]; then make distcheck LDFLAGS="-L$(brew --prefix)/lib/" CFLAGS="-g -O2 -I$(brew --prefix)/include/" diff --git a/src/fglm/Makefile.am b/src/fglm/Makefile.am index dc3d956f..42173c6f 100644 --- a/src/fglm/Makefile.am +++ b/src/fglm/Makefile.am @@ -6,6 +6,7 @@ libfglm_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) EXTRA_DIST = fglm.h \ libfglm.h \ + aligned_alloc.h \ berlekamp_massey.c \ data_fglm.c \ fglm_core.c \ diff --git a/src/fglm/aligned_alloc.h b/src/fglm/aligned_alloc.h new file mode 100644 index 00000000..1553d05c --- /dev/null +++ b/src/fglm/aligned_alloc.h @@ -0,0 +1,50 @@ +/* This file is part of msolve. + * + * msolve is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * msolve 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. + * + * You should have received a copy of the GNU General Public License + * along with msolve. If not, see + * + * Authors: + * Jérémy Berthomieu + * Christian Eder + * Mohab Safey El Din */ + +#ifndef ALIGNED_ALLOC_HEADER_H +#define ALIGNED_ALLOC_HEADER_H + +#ifdef _WIN32 + +#include +#include + +static inline int posix_memalign(void **__memptr, size_t __alignment, size_t __size) +{ + void *p = _aligned_malloc(__size, __alignment); + if (!p) + { + return ENOMEM; + } + *__memptr = p; + return 0; +} +#endif + +static inline void posix_memalign_free(void *__p) +{ +#ifdef _WIN32 + _aligned_free(__p); +#else + free(__p); +#endif +} + +#endif /* ALIGNED_ALLOC_HEADER_H */ diff --git a/src/fglm/data_fglm.c b/src/fglm/data_fglm.c index 368c85e5..8720b963 100644 --- a/src/fglm/data_fglm.c +++ b/src/fglm/data_fglm.c @@ -26,14 +26,15 @@ #include #include +#include "aligned_alloc.h" static inline void free_sp_mat_fglm(sp_matfglm_t *mat){ if(mat!=NULL){ - free(mat->dense_mat); - free(mat->triv_idx); - free(mat->triv_pos); - free(mat->dense_idx); - free(mat->dst); + posix_memalign_free(mat->dense_mat); + posix_memalign_free(mat->triv_idx); + posix_memalign_free(mat->triv_pos); + posix_memalign_free(mat->dense_idx); + posix_memalign_free(mat->dst); free(mat); } } @@ -80,10 +81,10 @@ static inline fglm_data_t *allocate_fglm_data(szmat_t nrows, szmat_t ncols, szma static inline void free_fglm_data(fglm_data_t *data){ - free(data->vecinit); - free(data->res); - free(data->vecmult); - free(data->vvec); + posix_memalign_free(data->vecinit); + posix_memalign_free(data->res); + posix_memalign_free(data->vecmult); + posix_memalign_free(data->vvec); free(data->pts); free(data); } diff --git a/src/fglm/fglm_core.c b/src/fglm/fglm_core.c index 87f61455..b761d75d 100644 --- a/src/fglm/fglm_core.c +++ b/src/fglm/fglm_core.c @@ -70,6 +70,8 @@ double omp_get_wtime(void) { return realtime();} #include "../upolmat/nmod_poly_mat_pmbasis.c" #endif +#include "aligned_alloc.h" + void print_fglm_data( FILE *file, const md_t * const st, @@ -775,9 +777,9 @@ static void generate_matrix_sequence(sp_matfglm_t *matxn, fglm_data_t *data, RED_32, RED_64); } - free(Rmat); - free(res); - free(tres); + posix_memalign_free(Rmat); + posix_memalign_free(res); + posix_memalign_free(tres); } @@ -1544,7 +1546,7 @@ param_t *nmod_fglm_compute_trace_data(sp_matfglm_t *matrix, mod_t prime, #endif #if DEBUGFGLM >= 1 - FILE *fmat = fopen("/tmp/matrix.fglm", "w"); + FILE *fmat = fopen("/tmp/matrix.fglm", "wb"); display_fglm_matrix(fmat, matrix); fclose(fmat); #endif @@ -1749,7 +1751,7 @@ int nmod_fglm_compute_apply_trace_data(sp_matfglm_t *matrix, #endif #if DEBUGFGLM >= 1 - FILE *fmat = fopen("/tmp/matrix.fglm", "w"); + FILE *fmat = fopen("/tmp/matrix.fglm", "wb"); display_fglm_matrix(fmat, matrix); fclose(fmat); #endif @@ -2098,7 +2100,7 @@ param_t *nmod_fglm_guess_colon(sp_matfglmcol_t *matrix, #endif #if DEBUGFGLM >= 1 - FILE *fmat = fopen("/tmp/matrix.fglm", "w"); + FILE *fmat = fopen("/tmp/matrix.fglm", "wb"); display_fglm_colon_matrix(fmat, matrix); fclose(fmat); #endif diff --git a/src/msolve/Makefile.am b/src/msolve/Makefile.am index bfe5f6d9..cb5d0f7f 100644 --- a/src/msolve/Makefile.am +++ b/src/msolve/Makefile.am @@ -8,6 +8,7 @@ libmsolve_la_LIBADD = ../usolve/libusolve.la ../fglm/libfglm.la ../neogb/libneo EXTRA_DIST = msolve-data.h \ msolve.h \ + getdelim.h \ msolve-data.c \ duplicate.c \ hilbert.c \ diff --git a/src/msolve/duplicate.c b/src/msolve/duplicate.c index f1db24f0..5cf3742e 100644 --- a/src/msolve/duplicate.c +++ b/src/msolve/duplicate.c @@ -394,7 +394,7 @@ static inline void duplicate_data_mthread_gbtrace(int nthreads, trace_t **btrace){ - const len_t len = num_gb[0] * (st->nvars); + const len_t len = num_gb[0] * (st->nvars - st->nev); for(int i = 0; i < nthreads; i++){ leadmons_current[i] = (int32_t *)calloc(len, sizeof(int32_t)); diff --git a/src/msolve/getdelim.h b/src/msolve/getdelim.h new file mode 100644 index 00000000..70033b73 --- /dev/null +++ b/src/msolve/getdelim.h @@ -0,0 +1,131 @@ +/* getdelim.h --- Implementation of replacement getdelim/getline function. + Copyright (C) 1994, 1996-1998, 2001, 2003, 2005-2025 Free Software + Foundation, Inc. + + This file is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as + published by the Free Software Foundation; either version 2.1 of the + License, or (at your option) any later version. + + This file 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 Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program. If not, see . */ + +/* Ported from glibc by Simon Josefsson. */ + +#ifndef GETDELIM_HEADER_H +#define GETDELIM_HEADER_H + +#ifdef _WIN32 + +#include + +#include +#include +#include +#include + +static inline void +alloc_failed (void) +{ +#if defined _WIN32 && ! defined __CYGWIN__ + /* Avoid errno problem without using the realloc module; see: + https://lists.gnu.org/r/bug-gnulib/2016-08/msg00025.html */ + errno = ENOMEM; +#endif +} + +/* Read up to (and including) a DELIMITER from FP into *LINEPTR (and + NUL-terminate it). *LINEPTR is a pointer returned from malloc (or + NULL), pointing to *N characters of space. It is realloc'ed as + necessary. Returns the number of characters read (not including + the null terminator), or -1 on error or EOF. */ + +static inline ssize_t +getdelim (char **lineptr, size_t *n, int delimiter, FILE *fp) +{ + ssize_t result; + size_t cur_len = 0; + + if (lineptr == NULL || n == NULL || fp == NULL) + { + errno = EINVAL; + return -1; + } + + if (*lineptr == NULL || *n == 0) + { + char *new_lineptr; + *n = 120; + new_lineptr = (char *) realloc (*lineptr, *n); + if (new_lineptr == NULL) + { + alloc_failed (); + return -1; + } + *lineptr = new_lineptr; + } + + for (;;) + { + int i; + + i = getc (fp); + if (i == EOF) + { + result = -1; + break; + } + + /* Make enough space for len+1 (for final NUL) bytes. */ + if (cur_len + 1 >= *n) + { + size_t needed_max = + SSIZE_MAX < SIZE_MAX ? (size_t) SSIZE_MAX + 1 : SIZE_MAX; + size_t needed = 2 * *n + 1; /* Be generous. */ + char *new_lineptr; + + if (needed_max < needed) + needed = needed_max; + if (cur_len + 1 >= needed) + { + errno = EOVERFLOW; + return -1; + } + + new_lineptr = (char *) realloc (*lineptr, needed); + if (new_lineptr == NULL) + { + alloc_failed (); + return -1; + } + + *lineptr = new_lineptr; + *n = needed; + } + + (*lineptr)[cur_len] = i; + cur_len++; + + if (i == delimiter) + break; + } + (*lineptr)[cur_len] = '\0'; + result = cur_len ? cur_len : result; + + return result; +} + +static inline ssize_t +getline (char **lineptr, size_t *n, FILE *stream) +{ + return getdelim (lineptr, n, '\n', stream); +} + +#endif + +#endif /* GETDELIM_HEADER_H */ diff --git a/src/msolve/hilbert.c b/src/msolve/hilbert.c index af9b4174..8fb518dc 100644 --- a/src/msolve/hilbert.c +++ b/src/msolve/hilbert.c @@ -1535,10 +1535,10 @@ static inline sp_matfglm_t * build_matrixn(int32_t *lmb, long dquot, int32_t bld count++; if(len_xn < count && i < dquot){ fprintf(stderr, "One should not arrive here (build_matrix)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); free(matrix); free(len_gb_xn); @@ -1552,11 +1552,11 @@ static inline sp_matfglm_t * build_matrixn(int32_t *lmb, long dquot, int32_t bld fprintf(stderr, "Multiplication by "); display_monomial_full(stderr, nv, NULL, 0, exp); fprintf(stderr, " gets outside the staircase\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); matrix = NULL; @@ -1917,12 +1917,12 @@ build_matrixn_colon(int32_t *lmb, long dquot, int32_t bld, free(lens); free(exps); free(cfs); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->zero_idx); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->zero_idx); + posix_memalign_free(matrix->dst); free(matrix); free(evi); free(len_gb_xn); @@ -2375,12 +2375,12 @@ build_matrixn_colon_no_zero(int32_t *lmb, long dquot, int32_t bld, free(lens); free(exps); free(cfs); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); - /* free(matrix->zero_idx); */ + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); + /* posix_memalign_free(matrix->zero_idx); */ free(matrix); free(evi); free(len_gb_xn); @@ -2593,11 +2593,11 @@ static inline sp_matfglm_t * build_matrixn_trace(int32_t **bdiv_xn, count++; if(len_xn < count && i < dquot){ fprintf(stderr, "One should not arrive here (build_matrix)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -2611,11 +2611,11 @@ static inline sp_matfglm_t * build_matrixn_trace(int32_t **bdiv_xn, fprintf(stderr, "Multiplication by "); display_monomial_full(stderr, nv, NULL, 0, exp); fprintf(stderr, " gets outside the staircase\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -2790,11 +2790,11 @@ static inline sp_matfglm_t * build_matrixn_from_bs(int32_t *lmb, long dquot, count++; if(len_xn < count && i < dquot){ fprintf(stderr, "One should not arrive here (build_matrix with trace)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -2810,11 +2810,11 @@ static inline sp_matfglm_t * build_matrixn_from_bs(int32_t *lmb, long dquot, display_monomial_full(stderr, nv, NULL, 0, exp); #endif fprintf(stderr, " gets outside the staircase\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -2985,11 +2985,11 @@ static inline void build_matrixn_from_bs_trace_application(sp_matfglm_t *matrix, count++; if(len_xn < count && i < dquot){ fprintf(stderr, "One should not arrive here (build_matrix with trace)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -3003,11 +3003,11 @@ static inline void build_matrixn_from_bs_trace_application(sp_matfglm_t *matrix, fprintf(stderr, "Multiplication by "); display_monomial_full(stderr, nv, NULL, 0, exp); fprintf(stderr, " gets outside the staircase\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free(len_gb_xn); @@ -3159,11 +3159,11 @@ static inline void build_matrixn_unstable_from_bs_trace_application(sp_matfglm_t count++; if(len_xn < count && i < dquot){ fprintf(stderr, "One should not arrive here (build_matrix with trace)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free_basis_without_hash_table(&tbr); @@ -3188,11 +3188,11 @@ static inline void build_matrixn_unstable_from_bs_trace_application(sp_matfglm_t count_nf++; if (count_not_lm < count_nf && i < dquot) { fprintf(stderr, "One should not arrive here (build_matrix with trace)\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free_basis_without_hash_table(&tbr); @@ -3212,11 +3212,11 @@ static inline void build_matrixn_unstable_from_bs_trace_application(sp_matfglm_t fprintf(stderr, "Multiplication by "); display_monomial_full(stderr, nv, NULL, 0, exp); fprintf(stderr, " gets outside the staircase\n"); - free(matrix->dense_mat); - free(matrix->dense_idx); - free(matrix->triv_idx); - free(matrix->triv_pos); - free(matrix->dst); + posix_memalign_free(matrix->dense_mat); + posix_memalign_free(matrix->dense_idx); + posix_memalign_free(matrix->triv_idx); + posix_memalign_free(matrix->triv_pos); + posix_memalign_free(matrix->dst); free(matrix); free_basis_without_hash_table(&tbr); diff --git a/src/msolve/iofiles.c b/src/msolve/iofiles.c index 8cf823dd..0c0c1975 100644 --- a/src/msolve/iofiles.c +++ b/src/msolve/iofiles.c @@ -18,6 +18,8 @@ * Christian Eder * Mohab Safey El Din */ +#include "getdelim.h" + static inline void store_exponent(const char *term, data_gens_ff_t *gens, int32_t pos) { len_t i, j, k; diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c index e25c2eb2..7f1a0527 100644 --- a/src/msolve/lifting-gb.c +++ b/src/msolve/lifting-gb.c @@ -1833,7 +1833,7 @@ void print_msolve_gbtrace_qq(data_gens_ff_t *gens, FILE *ofile; if (flags->files->out_file != NULL) { - ofile = fopen(flags->files->out_file, "w+"); + ofile = fopen(flags->files->out_file, "wb+"); } else { ofile = stdout; } @@ -1870,7 +1870,7 @@ void print_msolve_gbtrace_qq(data_gens_ff_t *gens, if(flags->print_gb > 1){ if(flags->files->out_file != NULL){ - FILE *ofile = fopen(flags->files->out_file, "a+"); + FILE *ofile = fopen(flags->files->out_file, "ab+"); display_gbmodpoly_cf_qq(ofile, (*modgbsp), gens); fclose(ofile); } @@ -1880,7 +1880,7 @@ void print_msolve_gbtrace_qq(data_gens_ff_t *gens, } if(flags->print_gb == 1){ if(flags->files->out_file != NULL){ - FILE *ofile = fopen(flags->files->out_file, "a+"); + FILE *ofile = fopen(flags->files->out_file, "ab+"); display_lm_gbmodpoly_cf_qq(ofile, (*modgbsp), gens); fclose(ofile); } diff --git a/src/msolve/main.c b/src/msolve/main.c index 87607fee..49eae86c 100644 --- a/src/msolve/main.c +++ b/src/msolve/main.c @@ -405,7 +405,7 @@ int main(int argc, char **argv){ /* clear out_file if given */ if(files->out_file != NULL){ - FILE *ofile = fopen(files->out_file, "w"); + FILE *ofile = fopen(files->out_file, "wb"); if(ofile == NULL){ fprintf(stderr, "Cannot open output file\n"); exit(1); diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index 4f2623c9..d95f03be 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -1550,7 +1550,7 @@ static inline void print_groebner_basis(files_gb *files, if (md->print_gb) { int32_t gfc = md->gfc; md->gfc = fc; - print_ff_basis_data(files->out_file, "a", bs, bs->ht, md, gens, + print_ff_basis_data(files->out_file, "ab", bs, bs->ht, md, gens, md->print_gb); md->gfc = gfc; } @@ -3787,7 +3787,7 @@ int real_msolve_qq(mpz_param_t *mpz_paramp, param_t **nmod_param, int *dim_ptr, void display_arrays_of_real_roots(files_gb *files, int32_t len, real_point_t **lreal_pts, long *lnbr) { if (files->out_file != NULL) { - FILE *ofile = fopen(files->out_file, "a+"); + FILE *ofile = fopen(files->out_file, "ab+"); fprintf(ofile, "["); for (int i = 0; i < len - 1; i++) { display_real_points(ofile, lreal_pts[i], lnbr[i]); @@ -3814,7 +3814,7 @@ void display_output(int b, int dim, int dquot, files_gb *files, real_point_t **real_pts_ptr, int info_level) { if (dquot == 0) { if (files->out_file != NULL) { - FILE *ofile = fopen(files->out_file, "a+"); + FILE *ofile = fopen(files->out_file, "ab+"); fprintf(ofile, "[-1]:\n"); fclose(ofile); } else { @@ -3826,7 +3826,7 @@ void display_output(int b, int dim, int dquot, files_gb *files, if (dim == 0 && dquot >= 0) { (*mpz_paramp)->nvars = gens->nvars; if (files->out_file != NULL) { - FILE *ofile = fopen(files->out_file, "a+"); + FILE *ofile = fopen(files->out_file, "ab+"); fprintf(ofile, "[0, "); if (get_param >= 1 || gens->field_char) { mpz_param_out_str_maple(ofile, gens, dquot, *mpz_paramp, param); @@ -3858,7 +3858,7 @@ void display_output(int b, int dim, int dquot, files_gb *files, fprintf(stdout, "The ideal has positive dimension\n"); } if (files->out_file != NULL) { - FILE *ofile2 = fopen(files->out_file, "a+"); + FILE *ofile2 = fopen(files->out_file, "ab+"); // 1 because dim is >0 fprintf(ofile2, "[1, %d, -1, []]:\n", gens->nvars); fclose(ofile2); @@ -4097,7 +4097,7 @@ int core_msolve( return -1; } if (print_gb) { - print_ff_basis_data(files->out_file, "a", bs, bht, + print_ff_basis_data(files->out_file, "ab", bs, bht, st, gens, print_gb); } } @@ -4574,7 +4574,7 @@ int core_msolve( } /* print all reduced elements in tbr, first normal_form ones * are the input elements */ - print_ff_nf_data(files->out_file, "a", 0, normal_form, tbr, bht, st, gens, 2); + print_ff_nf_data(files->out_file, "ab", 0, normal_form, tbr, bht, st, gens, 2); if (normal_form_matrix > 0) { /* sht and hcm will store the union of the support * of all normal forms in tbr. */ diff --git a/src/neogb/sort_r.h b/src/neogb/sort_r.h index 147d3834..bc6842af 100644 --- a/src/neogb/sort_r.h +++ b/src/neogb/sort_r.h @@ -47,7 +47,7 @@ void sort_r(void *base, size_t nel, size_t width, #if (defined __APPLE__ || defined __MACH__ || defined __DARWIN__ || \ (defined __FreeBSD__ && !defined(qsort_r)) || defined __DragonFly__) # define _SORT_R_BSD -#elif (defined __MINGW32__ || defined __GLIBC__ || \ +#elif (defined __GLIBC__ || \ (defined (__FreeBSD__) && defined(qsort_r))) # define _SORT_R_LINUX diff --git a/src/usolve/refine.c b/src/usolve/refine.c index d07759bc..397a4dbf 100644 --- a/src/usolve/refine.c +++ b/src/usolve/refine.c @@ -38,6 +38,22 @@ static void make_exact_root(mpz_t *upol, unsigned long int *deg, } } +uint64_t mpz_get_ull(const mpz_t op) { + if (mp_bits_per_limb == 64) { + mp_limb_t limb = mpz_getlimbn(op, 0); + return limb; + } + if (mpz_size(op) == 0) { + return 0; + } + uint64_t limb_lo = (uint64_t)(mpz_getlimbn(op, 0)); + uint64_t limb_hi = 0; + if (mpz_size(op) >= 2) { + limb_hi = (uint64_t)(mpz_getlimbn(op, 1)); + } + return (limb_hi << 32) | limb_lo; +} + /* computes 2^logN f(a) / (f(a) - f(b)) where f(a) = vala, f(b) = valb */ static long long int index_linearinterp(mpz_t *vala, mpz_t *valb, mpz_t *q, long long int logN){ @@ -55,16 +71,31 @@ static long long int index_linearinterp(mpz_t *vala, mpz_t *valb, mpz_t *q, } return -1; } - long long int index = mpz_get_ui(*q); + long long int index = mpz_get_ull(*q); return index; } +void mpz_init_set_ull(mpz_t rop, uint64_t op) { + if (sizeof(mp_bitcnt_t) == 4) { + uint32_t hi = (uint32_t)(op >> 32); + uint32_t lo = (uint32_t)(op); + mpz_init_set_ui(rop, hi); + mpz_mul_2exp(rop, rop, 32); + mpz_add_ui(rop, rop, lo); + } else { + mpz_init_set_ui(rop, op); + } +} + static void getx_and_eval_2exp(mpz_t *upol, unsigned long int deg, mpz_t *a, mpz_t *x, mpz_t *value, mpz_t *q, - int k, long int index){ + int k, long long int index){ + mpz_t index_gmp; + mpz_init_set_ull(index_gmp, index); mpz_set(*x, *a); - mpz_add_ui(*x, *x, index); + mpz_add(*x, *x, index_gmp); mpz_poly_eval_2exp_naive(upol, deg, x, k, value, q); + mpz_clear(index_gmp); } static void getx_and_eval_2exp_mpzidx(mpz_t *upol, unsigned long int deg, @@ -77,21 +108,26 @@ static void getx_and_eval_2exp_mpzidx(mpz_t *upol, unsigned long int deg, static void getx_and_eval_2expleft(mpz_t *upol, unsigned long int deg, mpz_t *a, mpz_t *x, mpz_t *value, mpz_t *q, - int k, long int index){ + int k, long long int index){ + mpz_t index_gmp; + mpz_init_set_ull(index_gmp, index); mpz_set(*x, *a); - mpz_sub_ui(*x, *x, index); + mpz_sub(*x, *x, index_gmp); mpz_poly_eval_2exp_naive(upol, deg, x, k, value, q); + mpz_clear(index_gmp); } /* assumes k >= -2 */ static void getx_and_eval(mpz_t *upol, unsigned long int deg, mpz_t *a, mpz_t *x, mpz_t *value, mpz_t *q, - long int Nlog, long int k, long int index){ + long int Nlog, long int k, long long int index){ + mpz_t index_gmp; + mpz_init_set_ull(index_gmp, index); if(k-Nlog>=0){ mpz_set_ui(*x, 1); mpz_mul_2exp(*x, *x, k - Nlog); - mpz_mul_ui(*x,*x,index); + mpz_mul(*x,*x,index_gmp); mpz_add(*x, *x, *a); mpz_poly_eval_2exp_naive(upol, deg, x, 0, value, q); @@ -99,10 +135,11 @@ static void getx_and_eval(mpz_t *upol, unsigned long int deg, else{ mpz_set(*x, *a); mpz_mul_2exp(*x, *x, Nlog - k); - mpz_add_ui(*x, *x, index); + mpz_add(*x, *x, index_gmp); mpz_poly_eval_2exp_naive(upol, deg, x, Nlog - k, value, q); } + mpz_clear(index_gmp); } static void getx_and_eval_mpzidx(mpz_t *upol, unsigned long int deg, @@ -287,7 +324,7 @@ static void refine_root_by_N_positive_k(mpz_t *upol, unsigned long int *deg_ptr, mpz_set(*vala, *tmpvala); mpz_set(*valb, *tmpvalb); - int64_t maxindex = (1L<<(Nlog)); + int64_t maxindex = (1LL<<(Nlog)); if(index == -2 || index == 0 || (LOG2(index) > Nlog && index > 0) ){ if(Nlog == 2) index = 2;