1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_expand_double_array: convert double vector to mpz
3 //------------------------------------------------------------------------------
4
5 // SLIP_LU: (c) 2019-2020, Chris Lourenco, Jinhao Chen, Erick Moreno-Centeno,
6 // Timothy A. Davis, Texas A&M University. All Rights Reserved. See
7 // SLIP_LU/License for the license.
8
9 //------------------------------------------------------------------------------
10
11 /* Purpose: This function converts a double array of size n to an appropriate
12 * mpz array of size n. To do this, the number is multiplied by an appropriate
13 * power, then, the GCD is found. This function allows the use of matrices in
14 * double precision to work with SLIP LU.
15 *
16 */
17
18 #define SLIP_FREE_ALL \
19 SLIP_MPZ_CLEAR(gcd); \
20 SLIP_MPZ_CLEAR(one); \
21 SLIP_MPQ_CLEAR(temp); \
22 SLIP_matrix_free(&x3, NULL); \
23
24 #include "slip_internal.h"
25
slip_expand_double_array(mpz_t * x_out,double * x,mpq_t scale,int64_t n,const SLIP_options * option)26 SLIP_info slip_expand_double_array
27 (
28 mpz_t* x_out, // integral final array
29 double* x, // double array that needs to be made integral
30 mpq_t scale, // the scaling factor used (x_out = scale * x)
31 int64_t n, // size of x
32 const SLIP_options* option // Command options
33 )
34 {
35
36 //--------------------------------------------------------------------------
37 // check inputs
38 //--------------------------------------------------------------------------
39
40 // inputs have been checked in the only caller slip_cast_array
41
42 //--------------------------------------------------------------------------
43
44 int64_t i, k ;
45 int r1, r2 = 1;
46 bool nz_found = false;
47 SLIP_info info ;
48 // Double precision accurate to about 2e-16. We multiply by 10e17 to convert
49 // (overestimate to be safe)
50 double expon = pow(10, 17);
51 // Quad precision in case input is huge
52 SLIP_matrix* x3 = NULL;
53 mpz_t gcd, one; SLIP_MPZ_SET_NULL(gcd); SLIP_MPZ_SET_NULL(one);
54 mpq_t temp; SLIP_MPQ_SET_NULL(temp);
55
56 mpfr_rnd_t round = SLIP_OPTION_ROUND (option) ;
57
58 SLIP_CHECK(SLIP_mpq_init(temp));
59 SLIP_CHECK(SLIP_mpz_init(gcd));
60 SLIP_CHECK(SLIP_mpz_init(one));
61
62 SLIP_CHECK (SLIP_matrix_allocate(&x3, SLIP_DENSE, SLIP_MPFR, n, 1, n,
63 false, true, option));
64
65 SLIP_CHECK(SLIP_mpq_set_d(scale, expon)); // scale = 10^17
66 for (i = 0; i < n; i++)
67 {
68 // Set x3[i] = x[i]
69 SLIP_CHECK(SLIP_mpfr_set_d(x3->x.mpfr[i], x[i], round));
70
71 // x3[i] = x[i] * 10^17
72 SLIP_CHECK(SLIP_mpfr_mul_d(x3->x.mpfr[i], x3->x.mpfr[i], expon, round));
73
74 // x_out[i] = x3[i]
75 SLIP_CHECK(SLIP_mpfr_get_z(x_out[i], x3->x.mpfr[i], round));
76 }
77
78 //--------------------------------------------------------------------------
79 // Compute the GCD to reduce the size of scale
80 //--------------------------------------------------------------------------
81
82 SLIP_CHECK(SLIP_mpz_set_ui(one, 1));
83 // Find an initial GCD
84 for (i = 0; i < n; i++)
85 {
86 if (!nz_found)
87 {
88 SLIP_CHECK(SLIP_mpz_cmp_ui(&r1, x_out[i], 0)); // Check if x[i] == 0
89 if (r1 != 0)
90 {
91 nz_found = true;
92 k = i;
93 SLIP_CHECK(SLIP_mpz_set(gcd, x_out[i]));
94 }
95 }
96 else
97 {
98 // Compute the GCD, stop if gcd == 1
99 SLIP_CHECK(SLIP_mpz_gcd(gcd, gcd, x_out[i]));
100 SLIP_CHECK(SLIP_mpz_cmp(&r2, gcd, one));
101 if (r2 == 0)
102 {
103 break;
104 }
105 }
106 }
107
108 if (!nz_found) // Array is all zeros
109 {
110 SLIP_FREE_ALL;
111 SLIP_mpq_set_z(scale, one);
112 return SLIP_OK;
113 }
114
115 //--------------------------------------------------------------------------
116 // Scale all entries to make as small as possible
117 //--------------------------------------------------------------------------
118
119 if (r2 != 0) // If gcd == 1 then stop
120 {
121 for (i = k; i < n; i++)
122 {
123 SLIP_CHECK(SLIP_mpz_divexact(x_out[i], x_out[i], gcd));
124 }
125 SLIP_CHECK(SLIP_mpq_set_z(temp, gcd));
126 SLIP_CHECK(SLIP_mpq_div(scale, scale, temp));
127 }
128 SLIP_FREE_ALL;
129 return SLIP_OK;
130 }
131
132