newton_test.c

   1 // Test Gosper's example: find coth 1/2 via Newton's method.
   2
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6 #include <gmp.h>
   7 #include "cf.h"
   8 #include "test.h"
   9
  10 static void *cot1_fn(cf_t cf) {
  11   int n = 1;
  12   while(cf_wait(cf)) {
  13     cf_put_int(cf, n);
  14     n += 2;
  15   }
  16   return NULL;
  17 }
  18
  19 int main() {
  20   mpz_t a[6];
  21   int i;
  22   for (i = 0; i < 6; i++) mpz_init(a[i]);
  23   cf_t x, y;
  24
  25   // Gosper's example.
  26   mpz_set_si(a[0], 1);
  27   mpz_set_si(a[3], -1);
  28   mpz_set_si(a[5], 1);
  29   x = cf_new_const(cot1_fn);
  30   y = cf_new_newton(x, a, a[0]);
  31   CF_EXPECT_DEC(y, "2.16395341373865284877");
  32   cf_free(x);
  33   cf_free(y);
  34
  35   // Confirm sqrt(1-(sin 1)^2) = cos 1
  36   x = cf_new_sin1();
  37   cf_t s1[2];
  38   cf_tee(s1, x);
  39   mpz_t b[8];
  40   mpz8_init(b);
  41   mpz8_set_int(b,
  42       -1, 0, 0, 1,
  43        0, 0, 0, 1);
  44
  45   cf_t bi = cf_new_bihom(s1[0], s1[1], b);
  46   cf_t n = cf_new_sqrt(bi);
  47
  48   CF_EXPECT_DEC(n, "0.54030230586813971740");
  49   cf_free(x);
  50   cf_free(s1[0]);
  51   cf_free(s1[1]);
  52   cf_free(bi);
  53   cf_free(n);
  54
  55   x = cf_new_sqrt_int(355, 113);
  56   CF_EXPECT_DEC(x, "1.77245392615830279609");
  57   cf_free(x);
  58
  59   mpz8_clear(b);
  60   for (i = 0; i < 6; i++) mpz_clear(a[i]);
  61   return 0;
  62 }