1/* Test file for mpfr_pow. 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include <stdio.h> 24#include <limits.h> 25#include <stdlib.h> 26 27#include "mpfr-test.h" 28 29 30int 31main (int argc, char *argv[]) 32{ 33 mpfr_t x, y, z; 34 35 mpfr_prec_t prec, yprec; 36 mpfr_t t, s; 37 mpfr_rnd_t rnd; 38 int inexact, compare, compare2; 39 unsigned int n, err; 40 41 mpfr_prec_t p0=2, p1=100; 42 unsigned int N=25; 43 44 tests_start_mpfr (); 45 46 mpfr_init (x); 47 mpfr_init2 (y,sizeof(unsigned long int)*CHAR_BIT); 48 mpfr_init (z); 49 50 mpfr_init (s); 51 mpfr_init (t); 52 53 /* generic test */ 54 for (prec = p0; prec <= p1; prec++) 55 { 56 mpfr_set_prec (x, prec); 57 mpfr_set_prec (s, sizeof(unsigned long int)*CHAR_BIT); 58 mpfr_set_prec (z, prec); 59 mpfr_set_prec (t, prec); 60 yprec = prec + 10; 61 62 for (n=0; n<N; n++) 63 { 64 mpfr_urandomb (x, RANDS); 65 mpfr_urandomb (s, RANDS); 66 if (randlimb () % 2) 67 mpfr_neg (s, s, MPFR_RNDN); 68 rnd = RND_RAND (); 69 mpfr_set_prec (y, yprec); 70 compare = mpfr_pow (y, x, s, rnd); 71 err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec; 72 if (mpfr_can_round (y, err, rnd, rnd, prec)) 73 { 74 mpfr_set (t, y, rnd); 75 inexact = mpfr_pow (z, x, s, rnd); 76 if (mpfr_cmp (t, z)) 77 { 78 printf ("results differ for x^y with x="); 79 mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN); 80 printf (" y="); 81 mpfr_out_str (stdout, 2, 0, s, MPFR_RNDN); 82 printf (" prec=%u rnd_mode=%s\n", (unsigned int) prec, 83 mpfr_print_rnd_mode (rnd)); 84 printf ("got "); 85 mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN); 86 puts (""); 87 printf ("expected "); 88 mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN); 89 puts (""); 90 printf ("approx "); 91 mpfr_print_binary (y); 92 puts (""); 93 exit (1); 94 } 95 compare2 = mpfr_cmp (t, y); 96 /* if rounding to nearest, cannot know the sign of t - f(x) 97 because of composed rounding: y = o(f(x)) and t = o(y) */ 98 if ((rnd != MPFR_RNDN) && (compare * compare2 >= 0)) 99 compare = compare + compare2; 100 else 101 compare = inexact; /* cannot determine sign(t-f(x)) */ 102 if (((inexact == 0) && (compare != 0)) || 103 ((inexact > 0) && (compare <= 0)) || 104 ((inexact < 0) && (compare >= 0))) 105 { 106 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d" 107 "\n", mpfr_print_rnd_mode (rnd), compare, inexact); 108 printf ("x="); mpfr_print_binary (x); puts (""); 109 printf ("y="); mpfr_print_binary (y); puts (""); 110 printf ("t="); mpfr_print_binary (t); puts (""); 111 exit (1); 112 } 113 } 114 } 115 } 116 117 mpfr_clear (s); 118 mpfr_clear (t); 119 120 mpfr_clear (x); 121 mpfr_clear (y); 122 mpfr_clear (z); 123 124 tests_end_mpfr (); 125 return 0; 126} 127