src/test/prob_distr_mpfr_ref.c

   1 /* Copyright 2012-2021, The Tor Project, Inc
   2  * See LICENSE for licensing information */
   3
   4 /** prob_distr_mpfr_ref.c
   5  *
   6  * Example reference file for GNU MPFR vectors tested in test_prob_distr.c .
   7  * Code by Riastradh.
   8  */
   9
  10 #include <complex.h>
  11 #include <float.h>
  12 #include <math.h>
  13 #include <stdio.h>
  14
  15 /* Must come after <stdio.h> so we get mpfr_printf.  */
  16 #include <mpfr.h>
  17
  18 /*  gcc -o mpfr prob_distr_mpfr_ref.c -lmpfr -lm */
  19
  20 /* Computes logit(p) for p = .49999 */
  21 int
  22 main(void)
  23 {
  24   mpfr_t p, q, r;
  25   mpfr_init(p);
  26   mpfr_set_prec(p, 200);
  27   mpfr_init(q);
  28   mpfr_set_prec(q, 200);
  29   mpfr_init(r);
  30   mpfr_set_prec(r, 200);
  31   mpfr_set_d(p, .49999, MPFR_RNDN);
  32   mpfr_set_d(q, 1, MPFR_RNDN);
  33   /* r := q - p = 1 - p */
  34   mpfr_sub(r, q, p, MPFR_RNDN);
  35   /* q := p/r = p/(1 - p) */
  36   mpfr_div(q, p, r, MPFR_RNDN);
  37   /* r := log(q) = log(p/(1 - p)) */
  38   mpfr_log(r, q, MPFR_RNDN);
  39   mpfr_printf("mpfr 200-bit\t%.128Rg\n", r);
  40
  41   /*
  42    * Print a double approximation to logit three different ways.  All
  43    * three agree bit for bit on the libms I tried, with the nextafter
  44    * adjustment (which is well within the 10 eps relative error bound
  45    * advertised).  Apparently I must have used the Goldberg expression
  46    * for what I wrote down in the test case.
  47    */
  48   printf("mpfr 53-bit\t%.17g\n", nextafter(mpfr_get_d(r, MPFR_RNDN), 0), 0);
  49   volatile double p0 = .49999;
  50   printf("log1p\t\t%.17g\n", nextafter(-log1p((1 - 2*p0)/p0), 0));
  51   volatile double x = (1 - 2*p0)/p0;
  52   volatile double xp1 = x + 1;
  53   printf("Goldberg\t%.17g\n", -x*log(xp1)/(xp1 - 1));
  54
  55   /*
  56    * Print a bad approximation, using the naive expression, to see a
  57    * lot of wrong digits, far beyond the 10 eps relative error attained
  58    * by -log1p((1 - 2*p)/p).
  59    */
  60   printf("naive\t\t%.17g\n", log(p0/(1 - p0)));
  61
  62   fflush(stdout);
  63   return ferror(stdout);
  64 }