/* Calculate the balanced ternary approximation of Pi (or any other
   constant by first using mpfr to convert to base 3, then to
   balanced ternary.  If the base 3 representation ends in ...222,
   then there may be a carry from the digits that were stripped,
   leading to inaccuracy in the final digits.  The algorithm here
   simply calculates more digits.  If a carry could have changed the
   result, it aborts with an error message.

   Thomas König, 2020-01-14.  This file is in the public domain.  */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <mpfr.h>

int main()
{
  mpfr_t pi;
  mpfr_exp_t ex;
  size_t len;
  int *buf;
  char *str;
  long int i, bits;
  long int trits = 79;
  long int guard_trits = 20;
  char *digits;
  int fail;

  bits = (guard_trits + trits) * 1585/1000 + 40;  /* ln(3)/ln(2) == 1.58495... */
  digits = calloc (bits, 1);
  mpfr_init2 (pi, bits);
  mpfr_const_pi (pi, MPFR_RNDZ);
  str = mpfr_get_str (NULL, &ex, 3, trits + guard_trits, pi, MPFR_RNDZ);
  len = strlen (str);
  buf = malloc (sizeof(int) * (len+1));
  buf[0] = 0;

  fail = 1;
  for (i = trits; i < len; i++)
    {
      if (buf[i] != '2')
	{
	  fail = 0;
	  break;
	}
    }

  if (fail)
    {
      fprintf (stderr,"Insufficient guard trits!\n");
      return 1;
    }

  for (i = 0; i < len; i++)
    buf[i + 1] = str[i] - '0';

  for (i = len; i >= 1; i--)
    while (buf[i] > 1)
      {
	buf[i] -= 3;
	buf[i - 1] += 1;
      }

  if (buf[0])
    printf("%d, ", buf[0]);

  for (i=1; i<trits; i++)
    printf("%d, ", buf[i]);

  printf ("%d\n", buf[trits]);
  printf("Pi = ");
  if (buf[0])
    putchar ("T01"[buf[0]+1]);

  for (i = 1; i <= ex; i++)
    putchar ("T01"[buf[i]+1]);

  putchar ('.');
  for (i=ex+1; i<trits; i++)
    putchar ("T01"[buf[i]+1]);

  putchar ('\n');
  free (buf);
  free (digits);
  mpfr_clear (pi);
  mpfr_free_str (str);
  return 0;
}