1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_expand_mpfr_array: convert mpfr aray 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 mpfr array of size n and precision prec to
12 * an appropriate mpz array of size n. To do this, the number is multiplied by
13 * the appropriate power of 10 then the gcd is found. This function allows mpfr
14 * arrays to be used within SLIP LU.
15 */
16
17 #define SLIP_FREE_ALL \
18 SLIP_MPFR_CLEAR(expon); \
19 SLIP_MPZ_CLEAR(temp_expon); \
20 SLIP_MPZ_CLEAR(gcd); \
21 SLIP_MPZ_CLEAR(one); \
22 SLIP_MPQ_CLEAR(temp); \
23 SLIP_matrix_free(&x3, NULL); \
24
25 #include "slip_internal.h"
26
slip_expand_mpfr_array(mpz_t * x_out,mpfr_t * x,mpq_t scale,int64_t n,const SLIP_options * option)27 SLIP_info slip_expand_mpfr_array
28 (
29 mpz_t* x_out, // full precision mpz array
30 mpfr_t* x, // mpfr array to be expanded
31 mpq_t scale, // scaling factor used (x_out = scale*x)
32 int64_t n, // size of x
33 const SLIP_options *option // command options containing the prec
34 // and rounding for mpfr
35 )
36 {
37
38 //--------------------------------------------------------------------------
39 // Input has already been checked
40 //--------------------------------------------------------------------------
41
42 SLIP_info info ;
43
44 //--------------------------------------------------------------------------
45 // initializations
46 //--------------------------------------------------------------------------
47
48 int64_t i, k ;
49 int r1, r2 = 1 ;
50 bool nz_found = false;
51 mpfr_t expon; SLIP_MPFR_SET_NULL(expon);
52 mpz_t temp_expon, gcd, one;
53 SLIP_matrix* x3 = NULL;
54 SLIP_MPZ_SET_NULL(temp_expon);
55 SLIP_MPZ_SET_NULL(gcd);
56 SLIP_MPZ_SET_NULL(one);
57 mpq_t temp; SLIP_MPQ_SET_NULL(temp);
58
59 uint64_t prec = SLIP_OPTION_PREC (option) ;
60 mpfr_rnd_t round = SLIP_OPTION_ROUND (option) ;
61
62 SLIP_CHECK(SLIP_mpq_init(temp));
63 SLIP_CHECK(SLIP_mpfr_init2(expon, prec));
64 SLIP_CHECK(SLIP_mpz_init(temp_expon));
65 SLIP_CHECK(SLIP_mpz_init(gcd));
66 SLIP_CHECK(SLIP_mpz_init(one));
67
68 SLIP_CHECK (SLIP_matrix_allocate(&x3, SLIP_DENSE, SLIP_MPFR, n, 1, n,
69 false, true, option));
70
71 // expon = 2^prec (overestimate)
72 SLIP_CHECK(SLIP_mpfr_ui_pow_ui(expon, 2, prec, round)) ;
73
74 for (i = 0; i < n; i++)
75 {
76 // x3[i] = x[i]*2^prec
77 SLIP_CHECK(SLIP_mpfr_mul(x3->x.mpfr[i], x[i], expon, round));
78
79 // x_out[i] = x3[i]
80 SLIP_CHECK(SLIP_mpfr_get_z(x_out[i], x3->x.mpfr[i], round));
81 }
82 SLIP_CHECK(SLIP_mpfr_get_z(temp_expon, expon, round));
83 SLIP_CHECK(SLIP_mpq_set_z(scale, temp_expon));
84
85 //--------------------------------------------------------------------------
86 // Find the gcd to reduce scale
87 //--------------------------------------------------------------------------
88
89 SLIP_CHECK(SLIP_mpz_set_ui(one, 1));
90 // Find an initial GCD
91 for (i = 0; i < n; i++)
92 {
93 if (!nz_found)
94 {
95 SLIP_CHECK(SLIP_mpz_cmp_ui(&r1, x_out[i], 0));
96 if (r1 != 0)
97 {
98 nz_found = true;
99 k = i;
100 SLIP_CHECK(SLIP_mpz_set(gcd, x_out[i]));
101 }
102 }
103 else
104 {
105 // Compute the GCD of the numbers, stop if gcd == 1
106 SLIP_CHECK(SLIP_mpz_gcd(gcd, gcd, x_out[i]));
107 SLIP_CHECK(SLIP_mpz_cmp(&r2, gcd, one));
108 if (r2 == 0)
109 {
110 break;
111 }
112 }
113 }
114
115 if (!nz_found) // Array is all zeros
116 {
117 SLIP_FREE_ALL;
118 SLIP_mpq_set_z(scale, one);
119 return SLIP_OK;
120 }
121
122 //--------------------------------------------------------------------------
123 // Scale all entries to make as small as possible
124 //--------------------------------------------------------------------------
125
126 if (r2 != 0) // If gcd == 1 stop
127 {
128 for (i = k; i < n; i++)
129 {
130 SLIP_CHECK(SLIP_mpz_divexact(x_out[i],x_out[i],gcd));
131 }
132 SLIP_CHECK(SLIP_mpq_set_z(temp,gcd));
133 SLIP_CHECK(SLIP_mpq_div(scale,scale,temp));
134 }
135 SLIP_FREE_ALL;
136 return SLIP_OK;
137 }
138
139