<commit_before>/*=============================================================================

    This file is part of FLINT.

    FLINT 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.

    FLINT 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 FLINT; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA

=============================================================================*/
/******************************************************************************

    Copyright (C) 2011 <NAME>

******************************************************************************/

#include <stdlib.h>
#include "flint.h"
#include "fmpz.h"
#include "fmpz_vec.h"
#include "fmpz_mat.h"
#include "nmod_mat.h"
#include "ulong_extras.h"


int
fmpz_mat_solve_dixon(fmpz_mat_t X, fmpz_t mod,
                        const fmpz_mat_t A, const fmpz_mat_t B)
{
    long n, cols;
    mp_limb_t p;
    fmpz_t N, D, bound, ppow, tested_p;

    fmpz_mat_t x, y, d, Ay;
    nmod_mat_t Ainv, dmod, ymod;

    int success = 0;

    n = A->r;
    cols = B->c;

    if (!fmpz_mat_is_square(A))
    {
        printf("fmpz_mat_solve_dixon: nonsquare system matrix");
        abort();
    }

    if (fmpz_mat_is_empty(A) || fmpz_mat_is_empty(B))
        return 1;

    fmpz_init(bound);
    fmpz_init(N);
    fmpz_init(D);
    fmpz_init(tested_p);
    fmpz_init(ppow);

    fmpz_mat_init(x, n, cols);
    fmpz_mat_init(y, n, cols);
    fmpz_mat_init(Ay, n, cols);
    fmpz_mat_init_set(d, B);

    fmpz_mat_solve_bound(N, D, A, B);

    /* Compute inverse */
    p = n_nextprime(1UL << (FLINT_BITS - FLINT_BIT_COUNT(n) - 2), 0);
    fmpz_set_ui(tested_p, 1UL);
    while (1)
    {
        nmod_mat_init(Ainv, n, n, p);
        fmpz_mat_get_nmod_mat(Ainv, A);

        if (nmod_mat_inv(Ainv, Ainv))
        {
            success = 1;
            break;
        }

        fmpz_mul_ui(tested_p, tested_p, p);

        if (fmpz_cmp(tested_p, D) > 0)
        {
            success = 0;
            goto cleanup;
        }

        p = n_nextprime(p, 0);
        nmod_mat_clear(Ainv);
    }

    if (success)
    {
        nmod_mat_init(dmod, n, cols, p);
        nmod_mat_init(ymod, n, cols, p);

        /* TODO: if one of N and D is much smaller, we could use a tighter<commit_msg>afford a slightly larger modulus in fmpz_mat_solve_dixon (todo: move this computation to a function nmod_mat or whatever)<commit_after>/*=============================================================================

    This file is part of FLINT.

    FLINT 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.

    FLINT 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 FLINT; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA

=============================================================================*/
/******************************************************************************

    Copyright (C) 2011 <NAME>

******************************************************************************/

#include <stdlib.h>
#include "flint.h"
#include "fmpz.h"
#include "fmpz_vec.h"
#include "fmpz_mat.h"
#include "nmod_mat.h"
#include "ulong_extras.h"


int
fmpz_mat_solve_dixon(fmpz_mat_t X, fmpz_t mod,
                        const fmpz_mat_t A, const fmpz_mat_t B)
{
    long n, cols;
    mp_limb_t p;
    fmpz_t N, D, bound, ppow, tested_p;

    fmpz_mat_t x, y, d, Ay;
    nmod_mat_t Ainv, dmod, ymod;

    int success = 0;

    n = A->r;
    cols = B->c;

    if (!fmpz_mat_is_square(A))
    {
        printf("fmpz_mat_solve_dixon: nonsquare system matrix");
        abort();
    }

    if (fmpz_mat_is_empty(A) || fmpz_mat_is_empty(B))
        return 1;

    fmpz_init(bound);
    fmpz_init(N);
    fmpz_init(D);
    fmpz_init(tested_p);
    fmpz_init(ppow);

    fmpz_mat_init(x, n, cols);
    fmpz_mat_init(y, n, cols);
    fmpz_mat_init(Ay, n, cols);
    fmpz_mat_init_set(d, B);

    fmpz_mat_solve_bound(N, D, A, B);

    /* Compute inverse */
    p = n_nextprime(1UL << (FLINT_BITS - 4), 0);
    fmpz_set_ui(tested_p, 1UL);
    while (1)
    {
        nmod_mat_init(Ainv, n, n, p);
        fmpz_mat_get_nmod_mat(Ainv, A);

        if (nmod_mat_inv(Ainv, Ainv))
        {
            success = 1;
            break;
        }

        fmpz_mul_ui(tested_p, tested_p, p);

        if (fmpz_cmp(tested_p, D) > 0)
        {
            success = 0;
            goto cleanup;
        }

        p = n_nextprime(p, 0);
        nmod_mat_clear(Ainv);
    }

    if (success)
    {
        nmod_mat_init(dmod, n, cols, p);
        nmod_mat_init(ymod, n, cols, p);

        /* TODO: if one of N and D is much smaller, we could use a tighter
{'fraction_non_alphanumeric': 0.19719808476680262, 'fraction_numerical': 0.009576165986877105, 'mean_word_length': 2.836734693877551, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1421779', 'n_tokens_mistral': 2034, 'n_tokens_neox': 1903, 'n_words': 552}