tadd.c

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

  mpfr_add (u, x, t, GMP_RNDN);
  if (mpfr_cmp (u, x)) {
    fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n");
    exit (1);
  }

  mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53);
  mpfr_set_d (x, -5.03525136761487735093e-74, GMP_RNDN);
  mpfr_set_d (t, 8.51539046314262304109e-91, GMP_RNDN);
  mpfr_add (u, x, t, GMP_RNDN);
  if (mpfr_get_d1 (u) != -5.0352513676148773509283672e-74) {
    fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n");
    exit (1);
  }

  mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76);
  mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
  mpfr_set_str_raw(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111");
  mpfr_sub(u, x, t, GMP_RNDU);
  mpfr_set_str_raw(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100");
  if (mpfr_cmp(u,t)) {
    printf("expect "); mpfr_print_binary(t); putchar('\n');
    fprintf (stderr, "mpfr_add failed for precisions 53-76\n"); exit(1);
  }
  mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108);
  mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
  mpfr_set_str_raw(t, "-0.101100010111001000010111111101111101000111001111011110011010101111001001111000111011001110011000000000111111");
  mpfr_sub(u, x, t, GMP_RNDU);
  mpfr_set_str_raw(t, "0.101100010111001000010111111101110011111110101001101111011010101110100000001011000010101110011000000000111111");
  if (mpfr_cmp(u,t)) {
    printf("expect "); mpfr_print_binary(t); putchar('\n');
    fprintf(stderr, "mpfr_add failed for precisions 53-108\n"); exit(1);
  }
  mpfr_set_prec(x, 97); mpfr_set_prec(t, 97); mpfr_set_prec(u, 97);
  mpfr_set_str_raw(x, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010000000000000000E-39");
  mpfr_set_ui(t, 1, GMP_RNDN);
  mpfr_add(u, x, t, GMP_RNDN);
  mpfr_set_str_raw(x, "0.1000000000000000000000000000000000000000111110110000100000000101100011011110100000101111100010001E1");
  if (mpfr_cmp(u,x)) {
    fprintf(stderr, "mpfr_add failed for precision 97\n"); exit(1);
  }
  mpfr_set_prec(x, 128); mpfr_set_prec(t, 128); mpfr_set_prec(u, 128);
  mpfr_set_str_raw(x, "0.10101011111001001010111011001000101100111101000000111111111011010100001100011101010001010111111101111010100110111111100101100010E-4");
  mpfr_set(t, x, GMP_RNDN);
  mpfr_sub(u, x, t, GMP_RNDN);
  mpfr_set_prec(x, 96); mpfr_set_prec(t, 96); mpfr_set_prec(u, 96);
  mpfr_set_str_raw(x, "0.111000000001110100111100110101101001001010010011010011100111100011010100011001010011011011000010E-4");
  mpfr_set(t, x, GMP_RNDN);
  mpfr_sub(u, x, t, GMP_RNDN);
  mpfr_set_prec(x, 85); mpfr_set_prec(t, 85); mpfr_set_prec(u, 85);
  mpfr_set_str_raw(x, "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4");
  mpfr_set_str_raw(t, "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4");
  mpfr_sub(u, x, t, GMP_RNDU);
  mpfr_sub(x, x, t, GMP_RNDU);
  if (mpfr_cmp(x, u) != 0) {
    printf("Error in mpfr_sub: u=x-t and x=x-t give different results\n");
    exit(1);
  }
  if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] & 
      ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
    printf("Error in mpfr_sub: result is not msb-normalized (1)\n"); exit(1);
  }
  mpfr_set_prec(x, 65); mpfr_set_prec(t, 65); mpfr_set_prec(u, 65);
  mpfr_set_str_raw(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
  mpfr_set_str_raw(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
  mpfr_sub(u, x, t, GMP_RNDU);
  if (mpfr_get_d1 (u) != 9.4349060620538533806e167) { /* 2^558 */
    printf("Error (1) in mpfr_sub\n"); exit(1);
  }

  mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64);
  mpfr_set_str_raw(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
  mpfr_set_str_raw(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
  mpfr_add(u, x, t, GMP_RNDU);
  if ((MPFR_MANT(u)[0] & 1) != 1) { 
    printf("error in mpfr_add with rnd_mode=GMP_RNDU\n");
    printf("b=  "); mpfr_print_binary(x); putchar('\n');
    printf("c=  "); mpfr_print_binary(t); putchar('\n');
    printf("b+c="); mpfr_print_binary(u); putchar('\n');
    exit(1);
  }

  /* bug found by Norbert Mueller, 14 Sep 2000 */
  mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10);
  mpfr_set_str_raw(x, "0.10001001011011001111101100110100000101111010010111010111E-7");
  mpfr_set_str_raw(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7");
  mpfr_sub(u, x, t, GMP_RNDU);

  /* array bound write found by Norbert Mueller, 26 Sep 2000 */
  mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95);
  mpfr_set_str_raw(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33");
  mpfr_set_str_raw(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33");
  mpfr_add(u, x, t, GMP_RNDN);

  /* array bound writes found by Norbert Mueller, 27 Sep 2000 */
  mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23);
  mpfr_set_str_raw(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59");
  mpfr_set_str_raw(t, "-0.10000111101011111110010100010001101100011100110100000E-59");
  mpfr_sub(u, x, t, GMP_RNDN);
  mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160);
  mpfr_set_str_raw(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35");
  mpfr_set_str_raw(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35");
  mpfr_add(u, x, t, GMP_RNDN);
  mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207);
  mpfr_set_str_raw(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66");
  mpfr_set_str_raw(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66");
  mpfr_add(u, x, t, GMP_RNDN);
  mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223);
  mpfr_set_str_raw(x, "0.10000000000000000000000000000000E1");
  mpfr_set_str_raw(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0");
  mpfr_sub(u, x, t, GMP_RNDN);
  if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] & 
      ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
    printf("Error in mpfr_sub: result is not msb-normalized (2)\n"); exit(1);
  }

  /* bug found by Nathalie Revol, 21 March 2001 */
  mpfr_set_prec (x, 65);
  mpfr_set_prec (t, 65);
  mpfr_set_prec (u, 65);
  mpfr_set_str_raw (x, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
  mpfr_set_str_raw (t, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
  mpfr_sub (u, t, x, GMP_RNDU);
  if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] & 
      ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
    fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (3)\n");
    exit (1);
  }

  /* bug found by Fabrice Rouillier, 27 Mar 2001 */
  mpfr_set_prec (x, 107);
  mpfr_set_prec (t, 107);
  mpfr_set_prec (u, 107);
  mpfr_set_str_raw (x, "0.10111001001111010010001000000010111111011011011101000001001000101000000000000000000000000000000000000000000E315");
  mpfr_set_str_raw (t, "0.10000000000000000000000000000000000101110100100101110110000001100101011111001000011101111100100100111011000E350");
  mpfr_sub (u, x, t, GMP_RNDU);
  if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] & 
      ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
    fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (4)\n");
    exit (1);
  }
  
  /* checks that NaN flag is correctly reset */
  mpfr_set_d (t, 1.0, GMP_RNDN);
  mpfr_set_d (u, 1.0, GMP_RNDN);
  mpfr_set_nan (x);
  mpfr_add (x, t, u, GMP_RNDN);
  if (mpfr_cmp_ui (x, 2)) {
    fprintf (stderr, "Error in mpfr_add: 1+1 gives %e\n", mpfr_get_d1 (x));
    exit (1);
  }

  mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
}

/* check case when c does not overlap with a, but both b and c count
   for rounding */
void
check_case_1b (void)
{
  mpfr_t a, b, c;
  unsigned int prec_a, prec_b, prec_c, dif;

  mpfr_init (a);
  mpfr_init (b);
  mpfr_init (c);

  for (prec_a = 2; prec_a <= 64; prec_a++)
    {
      mpfr_set_prec (a, prec_a);
      for (prec_b = prec_a + 2; prec_b <= 64; prec_b++)
	{
	  dif = prec_b - prec_a;
	  mpfr_set_prec (b, prec_b);
	  /* b = 1 - 2^(-prec_a) + 2^(-prec_b) */
	  mpfr_set_ui (b, 1, GMP_RNDN);
	  mpfr_div_2exp (b, b, dif, GMP_RNDN);
	  mpfr_sub_ui (b, b, 1, GMP_RNDN);
	  mpfr_div_2exp (b, b, prec_a, GMP_RNDN);
	  mpfr_add_ui (b, b, 1, GMP_RNDN);
	  for (prec_c = dif; prec_c <= 64; prec_c++)
	    {
	      /* c = 2^(-prec_a) - 2^(-prec_b) */
	      mpfr_set_prec (c, prec_c);
	      mpfr_set_si (c, -1, GMP_RNDN);
	      mpfr_div_2exp (c, c, dif, GMP_RNDN);
	      mpfr_add_ui (c, c, 1, GMP_RNDN);
	      mpfr_div_2exp (c, c, prec_a, GMP_RNDN);
	      mpfr_add (a, b, c, GMP_RNDN);
	      if (mpfr_cmp_ui (a, 1) != 0)
		{
		  fprintf (stderr, "case (1b) failed for prec_a=%u, prec_b=%u, prec_c=%u\n", prec_a, prec_b, prec_c);
		  printf("b="); mpfr_print_binary(b); putchar('\n');
		  printf("c="); mpfr_print_binary(c); putchar('\n');
		  printf("a="); mpfr_print_binary(a); putchar('\n');
		  exit (1);
		}
	    }
	}
    }

  mpfr_clear (a);
  mpfr_clear (b);
  mpfr_clear (c);
}

/* check case when c overlaps with a */
void
check_case_2 (void)
{
  mpfr_t a, b, c, d;

  mpfr_init2 (a, 300);
  mpfr_init2 (b, 800);
  mpfr_init2 (c, 500);
  mpfr_init2 (d, 800);

  mpfr_set_str_raw(a, "1E110");  /* a = 2^110 */
  mpfr_set_str_raw(b, "1E900");  /* b = 2^900 */
  mpfr_set_str_raw(c, "1E500");  /* c = 2^500 */
  mpfr_add(c, c, a, GMP_RNDZ);   /* c = 2^500 + 2^110 */
  mpfr_sub(d, b, c, GMP_RNDZ);   /* d = 2^900 - 2^500 - 2^110 */
  mpfr_add(b, b, c, GMP_RNDZ);   /* b = 2^900 + 2^500 + 2^110 */
  mpfr_add(a, b, d, GMP_RNDZ);   /* a = 2^901 */
  if (mpfr_cmp_ui_2exp (a, 1, 901))
    {
      fprintf (stderr, "b + d fails for b=2^900+2^500+2^110, d=2^900-2^500-2^110\n");
      fprintf (stderr, "expected 1.0e901, got ");
      mpfr_out_str (stderr, 2, 0, a, GMP_RNDN);
      fprintf (stderr, "\n");
      exit (1);
    }

  mpfr_clear (a);
  mpfr_clear (b);
  mpfr_clear (c);
  mpfr_clear (d);
}

/* checks when source and destination are equal */
void
check_same (void)
{
  mpfr_t x;

  mpfr_init(x); mpfr_set_d(x, 1.0, GMP_RNDZ);
  mpfr_add(x, x, x, GMP_RNDZ);
  if (mpfr_get_d1 (x) != 2.0) {
    printf("Error when all 3 operands are equal\n"); exit(1);
  }
  mpfr_clear(x);
}

#define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z)
#define check53nan(x, y, r) checknan(x, y, r, 53, 53, 53); 

#define MAX_PREC 100

void
check_inexact (void)
{
  mpfr_t x, y, z, u;
  mp_prec_t px, py, pu, pz;
  int inexact, cmp;
  mp_rnd_t rnd;
  
  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);
  mpfr_init (u);

  mpfr_set_prec (x, 2);
  mpfr_set_str_raw (x, "0.1E-4");
  mpfr_set_prec (u, 33);
  mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
  mpfr_set_prec (y, 31);
  if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
    {
      fprintf (stderr, "Wrong inexact flag (2): expected 0, got %d\n", inexact);
      exit (1);
    }

  mpfr_set_prec (x, 2);
  mpfr_set_str_raw (x, "0.1E-4");
  mpfr_set_prec (u, 33);
  mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
  mpfr_set_prec (y, 28);
  if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
    {
      fprintf (stderr, "Wrong inexact flag (1): expected 0, got %d\n", inexact);
      exit (1);
    }

  for (px=2; px<MAX_PREC; px++)
    {
      mpfr_set_prec (x, px);
      mpfr_random (x);
      for (pu=2; pu<MAX_PREC; pu++)
	{
	  mpfr_set_prec (u, pu);
	  mpfr_random (u);
	  for (py=2; py<MAX_PREC; py++)