tadd.c

a very popular packet of cryptography tools,it encloses the most common used algorithm and protocols

/* Test file for mpfr_add and mpfr_sub.

Copyright 1999, 2000, 2001, 2002 Free Software Foundation.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The MPFR Library 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 Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#define N 100000

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "gmp.h"
#include "mpfr.h"
#include "mpfr-impl.h"
#include "mpfr-test.h"

void checknan _PROTO((double, double, mp_rnd_t, unsigned int, unsigned int, unsigned int)); 
void check3 _PROTO((double, double, mp_rnd_t));
void check4 _PROTO((double, double, mp_rnd_t));
void check5 _PROTO((double, mp_rnd_t));
void check2 _PROTO((double, int, double, int, int, int)); 
void check2a _PROTO((double, int, double, int, int, int, char *)); 
void check64 _PROTO((void)); 
void check_same _PROTO((void)); 
void check_case_1b _PROTO((void)); 
void check_case_2 _PROTO((void));
void check_inexact _PROTO((void));


/* Parameter "z1" of check() used to be last in the argument list, but that
   tickled a bug in 32-bit sparc gcc 2.95.2.  A "double" in that position is
   passed on the stack at an address which is 4mod8, but the generated code
   didn't take into account that alignment, resulting in bus errors.  The
   easiest workaround is to move it to the start of the arg list (where it's
   passed in registers), this macro does that.  FIXME: Change the actual
   calls to check(), rather than using a macro.  */

#define check(x,y,rnd_mode,px,py,pz,z1)  _check(x,y,z1,rnd_mode,px,py,pz)

/* checks that x+y gives the same results in double
   and with mpfr with 53 bits of precision */
void
_check (double x, double y, double z1, mp_rnd_t rnd_mode, unsigned int px, 
        unsigned int py, unsigned int pz)
{
  double z2; mpfr_t xx,yy,zz; int cert=0;

  mpfr_init2(xx, px);
  mpfr_init2(yy, py);
  mpfr_init2(zz, pz);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  mpfr_add(zz, xx, yy, rnd_mode);
#ifdef MPFR_HAVE_FESETROUND
  mpfr_set_machine_rnd_mode(rnd_mode);
  if (px==53 && py==53 && pz==53) cert=1;
#endif
  if (z1==0.0) z1=x+y; else cert=1;
  z2 = mpfr_get_d1 (zz);
  mpfr_set_d (yy, z2, GMP_RNDN);
  if (!mpfr_cmp (zz, yy) && cert && z1!=z2 && !(isnan(z1) && isnan(z2))) {
    printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
    printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
	   x, y, mpfr_print_rnd_mode(rnd_mode));
    exit(1);
  }
  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}

void
checknan (double x, double y, mp_rnd_t rnd_mode, unsigned int px, 
          unsigned int py, unsigned int pz)
{
  double z2; mpfr_t xx, yy, zz;

  mpfr_init2(xx, px);
  mpfr_init2(yy, py);
  mpfr_init2(zz, pz);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  mpfr_add(zz, xx, yy, rnd_mode);
#ifdef MPFR_HAVE_FESETROUND
  mpfr_set_machine_rnd_mode(rnd_mode);
#endif
  if (MPFR_IS_NAN(zz) == 0) { printf("Error, not an MPFR_NAN for xx = %1.20e, y = %1.20e\n", x, y); exit(1); }
  z2 = mpfr_get_d1 (zz);
  if (!isnan(z2)) { printf("Error, not a NaN after conversion, xx = %1.20e yy = %1.20e, got %1.20e\n", x, y, z2); exit(1); }

  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}

#ifdef MPFR_HAVE_FESETROUND
/* idem than check for mpfr_add(x, x, y) */
void
check3 (double x, double y, mp_rnd_t rnd_mode)
{
  double z1,z2; mpfr_t xx,yy; int neg;

  neg = LONG_RAND() % 2;
  mpfr_init2(xx, 53);
  mpfr_init2(yy, 53);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  if (neg) mpfr_sub(xx, xx, yy, rnd_mode);
  else mpfr_add(xx, xx, yy, rnd_mode);
  mpfr_set_machine_rnd_mode(rnd_mode);
  z1 = (neg) ? x-y : x+y;
  z2 = mpfr_get_d1 (xx);
  mpfr_set_d (yy, z2, GMP_RNDN);
  if (!mpfr_cmp (xx, yy) && z1!=z2 && !(isnan(z1) && isnan(z2))) {
    printf("expected result is %1.20e, got %1.20e\n",z1,z2);
    printf("mpfr_%s(x,x,y) failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",
	   (neg) ? "sub" : "add",x,y,rnd_mode);
    exit(1);
  }
  mpfr_clear(xx); mpfr_clear(yy);
}

/* idem than check for mpfr_add(x, y, x) */
void
check4 (double x, double y, mp_rnd_t rnd_mode)
{
  double z1, z2;
  mpfr_t xx, yy;
  int neg;

  neg = LONG_RAND() % 2;
  mpfr_init2(xx, 53);
  mpfr_init2(yy, 53);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  if (neg) mpfr_sub(xx, yy, xx, rnd_mode);
  else mpfr_add(xx, yy, xx, rnd_mode);
  mpfr_set_machine_rnd_mode(rnd_mode);
  z1 = (neg) ? y-x : x+y;
  z2 = mpfr_get_d1 (xx);
  mpfr_set_d (yy, z2, GMP_RNDN);
  /* check that xx is representable as a double and no overflow occurred */
  if ((mpfr_cmp (xx, yy) == 0) && (z1 != z2)) {
    printf("expected result is %1.20e, got %1.20e\n", z1, z2);
    printf("mpfr_%s(x,y,x) failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
	   (neg) ? "sub" : "add", x, y, mpfr_print_rnd_mode(rnd_mode));
    exit(1);
  }
  mpfr_clear(xx); mpfr_clear(yy);
}

/* idem than check for mpfr_add(x, x, x) */
void
check5 (double x, mp_rnd_t rnd_mode)
{
  double z1,z2; mpfr_t xx, yy; int neg;

  mpfr_init2(xx, 53);
  mpfr_init2(yy, 53);
  neg = LONG_RAND() % 2;
  mpfr_set_d(xx, x, rnd_mode);
  if (neg) mpfr_sub(xx, xx, xx, rnd_mode);
  else mpfr_add(xx, xx, xx, rnd_mode);
  mpfr_set_machine_rnd_mode(rnd_mode);
  z1 = (neg) ? x-x : x+x;
  z2 = mpfr_get_d1 (xx);
  mpfr_set_d (yy, z2, GMP_RNDN);
  /* check NaNs first since mpfr_cmp does not like them */
  if (!(isnan(z1) && isnan(z2)) && !mpfr_cmp (xx, yy) && z1!=z2)
    {
      printf ("expected result is %1.20e, got %1.20e\n",z1,z2);
      printf ("mpfr_%s(x,x,x) failed for x=%1.20e with rnd_mode=%s\n",
              (neg) ? "sub" : "add", x, mpfr_print_rnd_mode (rnd_mode));
      exit (1);
    }
  mpfr_clear(xx);
  mpfr_clear(yy);
}

void
check2 (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode)
{
  mpfr_t xx, yy, zz; double z,z2; int u;

  mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  mpfr_add(zz, xx, yy, rnd_mode);
  mpfr_set_machine_rnd_mode(rnd_mode);
  z = x+y; z2=mpfr_get_d1 (zz); u=ulp(z,z2);
  /* one ulp difference is possible due to composed rounding */
  if (px>=53 && py>=53 && pz>=53 && ABS(u)>1) { 
    printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
	   x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
    printf("got %1.20e\n",z2);
    printf("result should be %1.20e (diff=%d ulp)\n",z,u);
    mpfr_set_d(zz, z, rnd_mode);
    printf("i.e."); mpfr_print_binary(zz); putchar('\n');
    exit(1); }
  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}
#endif

void
check2a (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode,
	      char *res)
{
  mpfr_t xx, yy, zz;

  mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
  mpfr_set_d(xx, x, rnd_mode);
  mpfr_set_d(yy, y, rnd_mode);
  mpfr_add(zz, xx, yy, rnd_mode);
  mpfr_set_prec(xx, pz);
  mpfr_set_str(xx, res, 16, GMP_RNDN);
  if (mpfr_cmp(xx, zz)) { 
    printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
	   x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
    printf("got        "); mpfr_print_binary(zz); putchar('\n');
    printf("instead of "); mpfr_print_binary(xx); putchar('\n');
    exit(1); 
  }
  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}

void
check64 (void)
{
  mpfr_t x, t, u;

  mpfr_init (x);
  mpfr_init (t);
  mpfr_init (u);

  mpfr_set_prec (x, 29);
  mpfr_set_str_raw (x, "1.1101001000101111011010010110e-3");
  mpfr_set_prec (t, 58);
  mpfr_set_str_raw (t, "0.11100010011111001001100110010111110110011000000100101E-1");
  mpfr_set_prec (u, 29);
  mpfr_add (u, x, t, GMP_RNDD);
  mpfr_set_str_raw (t, "1.0101011100001000011100111110e-1");
  if (mpfr_cmp (u, t))
    {
      fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=29, prec(t)=58\n");
      printf ("expected "); mpfr_out_str (stdout, 2, 29, t, GMP_RNDN);
      putchar ('\n');
      printf ("got      "); mpfr_out_str (stdout, 2, 29, u, GMP_RNDN);
      putchar ('\n');
      exit(1);
    }

  mpfr_set_prec (x, 4);
  mpfr_set_str_raw (x, "-1.0E-2");
  mpfr_set_prec (t, 2);
  mpfr_set_str_raw (t, "-1.1e-2");
  mpfr_set_prec (u, 2);
  mpfr_add (u, x, t, GMP_RNDN);
  if (MPFR_MANT(u)[0] << 2)
    {
      fprintf (stderr, "result not normalized for prec=2\n");
      mpfr_print_binary (u); putchar ('\n');
      exit (1);
    }
  mpfr_set_str_raw (t, "-1.0e-1");
  if (mpfr_cmp (u, t))
    {
      fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=4, prec(t)=2\n");
      printf ("expected -1.0e-1\n");
      printf ("got      "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
      putchar ('\n');
      exit (1);
    }

  mpfr_set_prec (x, 8);
  mpfr_set_str_raw (x, "-0.10011010"); /* -77/128 */
  mpfr_set_prec (t, 4);
  mpfr_set_str_raw (t, "-1.110e-5"); /* -7/128 */
  mpfr_set_prec (u, 4);
  mpfr_add (u, x, t, GMP_RNDN); /* should give -5/8 */
  mpfr_set_str_raw (t, "-1.010e-1");
  if (mpfr_cmp (u, t)) {
    fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=8, prec(t)=4\n");
    printf ("expected -1.010e-1\n");
    printf ("got      "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
    putchar ('\n');
    exit (1);
  }

  mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54);
  mpfr_set_str_raw (x, "-0.11111100100000000011000011100000101101010001000111E-401");
  mpfr_set_str_raw (t, "0.10110000100100000101101100011111111011101000111000101E-464");