1/* 2 3Copyright 2012, 2016 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library test suite. 6 7The GNU MP Library test suite is free software; you can redistribute it 8and/or modify it under the terms of the GNU General Public License as 9published by the Free Software Foundation; either version 3 of the License, 10or (at your option) any later version. 11 12The GNU MP Library test suite is distributed in the hope that it will be 13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15Public License for more details. 16 17You should have received a copy of the GNU General Public License along with 18the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ 19 20#include <assert.h> 21#include <limits.h> 22#include <stdlib.h> 23#include <stdio.h> 24 25#include "testutils.h" 26 27#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) 28 29#define COUNT 10000 30 31static void 32test_2by1(const mpz_t u) 33{ 34 mpz_t m, p, t; 35 36 mpz_init (m); 37 mpz_init (p); 38 mpz_init (t); 39 40 assert (mpz_size (u) == 1); 41 42 mpz_set_ui (m, mpn_invert_limb (u->_mp_d[0])); 43 mpz_setbit (m, GMP_LIMB_BITS); 44 45 mpz_mul (p, m, u); 46 47 mpz_set_ui (t, 0); 48 mpz_setbit (t, 2* GMP_LIMB_BITS); 49 mpz_sub (t, t, p); 50 51 /* Should have 0 < B^2 - m u <= u */ 52 if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0) 53 { 54 fprintf (stderr, "mpn_invert_limb failed:\n"); 55 dump ("u", u); 56 dump ("m", m); 57 dump ("p", p); 58 dump ("t", t); 59 abort (); 60 } 61 mpz_clear (m); 62 mpz_clear (p); 63 mpz_clear (t); 64} 65 66static void 67test_3by2(const mpz_t u) 68{ 69 mpz_t m, p, t; 70 71 mpz_init (m); 72 mpz_init (p); 73 mpz_init (t); 74 75 assert (mpz_size (u) == 2); 76 77 mpz_set_ui (m, mpn_invert_3by2 (u->_mp_d[1], u[0]._mp_d[0])); 78 79 mpz_setbit (m, GMP_LIMB_BITS); 80 81 mpz_mul (p, m, u); 82 83 mpz_set_ui (t, 0); 84 mpz_setbit (t, 3 * GMP_LIMB_BITS); 85 mpz_sub (t, t, p); 86 87 /* Should have 0 < B^3 - m u <= u */ 88 if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0) 89 { 90 fprintf (stderr, "mpn_invert_3by2 failed:\n"); 91 dump ("u", u); 92 dump ("m", m); 93 dump ("p", p); 94 dump ("t", t); 95 abort (); 96 } 97 mpz_clear (m); 98 mpz_clear (p); 99 mpz_clear (t); 100} 101 102void 103testmain (int argc, char **argv) 104{ 105 unsigned i; 106 mpz_t u, m, p, t; 107 108 mpz_init (u); 109 mpz_init (m); 110 mpz_init (p); 111 mpz_init (t); 112 113 /* These values trigger 32-bit overflow of ql in mpn_invert_3by2. */ 114 if (GMP_LIMB_BITS == 64) 115 { 116 mpz_set_str (u, "80007fff3ffe0000", 16); 117 test_2by1 (u); 118 mpz_set_str (u, "80007fff3ffe000000000000000003ff", 16); 119 test_3by2 (u); 120 } 121 122 for (i = 0; i < COUNT; i++) 123 { 124 mini_urandomb (u, GMP_LIMB_BITS); 125 mpz_setbit (u, GMP_LIMB_BITS -1); 126 127 test_2by1 (u); 128 } 129 130 for (i = 0; i < COUNT; i++) 131 { 132 mini_urandomb (u, 2*GMP_LIMB_BITS); 133 mpz_setbit (u, 2*GMP_LIMB_BITS -1); 134 135 test_3by2 (u); 136 } 137 138 mpz_clear (u); 139} 140