/* zebre.c
 * Produce images of zebra numbers
 * (C) Eric Buchlin 2004
 * This program is free software; you can redistribute it and/or modify 
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation (http://www.fsf.org); either version 2 of 
 * the License, or (at your option) any later version.
 */

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <gmp.h>
#include <mpfr.h>
#include <png.h>
#include <math.h>

void
write_image (int sx, int sy, unsigned char * image, const char * filename);

void
usage (const char * argv0);

int
main (int argc, char ** argv) {
  int
    sx = 256,
    sy = 256;
  unsigned ibase = 10;
  unsigned bits_per_pix = ceil (log (ibase) / log (2)) + 1;
  mpfr_t a;

  float palette[16][3] = 
    {{0., 0., 0.},
     {.6, .6, .6},
     {0., 1., 0.},
     {.8, 0., .8},
     {.5, .5, 0.},
     {0., .8, .8},
     {1., 1., 1.},
     {1., 1., 0.},
     {1., .5, 0.},
     {0., 0., .7},
     {.1, .1, .1},
     {.2, .2, .2},
     {.3, .3, .3},
     {.4, .4, .4},
     {.5, .5, .5},
     {.8, .7, .6}};

  mpfr_t base;
  mpfr_t tmp, tmp2;
  int i, j;
  int digit;
  unsigned char r, g, b;
  unsigned char * image;
  char * digits;
  mp_exp_t exponent;
  int which_number = 1;
  int param = 49, param2 = 2; /* additional parameters for some numbers */

  assert (ibase <= 16); /* need a larger palette if this is not satisfied */

  if (argc < 4) {
    usage (argv[0]);
    exit (0);
  } else {
    sx = atoi (argv[1]);
    sy = atoi (argv[2]);
    which_number = atoi (argv[3]);
  }
  if (argc > 4) {
    param = atoi (argv[4]);
  }
  if (argc > 5) {
    param2 = atoi (argv[5]);
  }
  image = (unsigned char *)malloc (3 * sx * sy);
  digits = (char *)malloc (sx * sy + 2);

  mpfr_set_default_prec (sx * sy * bits_per_pix);
  mpfr_init (a);
  mpfr_init (tmp);
  mpfr_init (tmp2);
  mpfr_init_set_ui (base, ibase, GMP_RNDN);

  switch (which_number) {
    case 1:
     printf ("Yéléhada, n=%d, k=%d\n", param, param2);
     mpfr_set_ui (a, ibase, GMP_RNDN);
     mpfr_pow_ui (a, a, param, GMP_RNDN);
     mpfr_sub_ui (a, a, 1, GMP_RNDN);
     mpfr_set_ui (tmp, 1, GMP_RNDN);
     mpfr_div_ui (tmp, tmp, param2, GMP_RNDN);
     mpfr_pow (a, a, tmp, GMP_RNDN);
     break;
    case 2:
     printf ("355/ 113\n");
     mpfr_set_ui (a, 355, GMP_RNDN);
     mpfr_set_ui (tmp, 113, GMP_RNDN);
     mpfr_div (a, a, tmp, GMP_RNDN);
     break;
    case 3:
     printf ("Mad zebra (R. Israel), n=%d\n", param);
     {
       mpfr_set_ui (tmp, ibase * ibase, GMP_RNDN);
       mpfr_set_ui (tmp2, param, GMP_RNDN);
       mpfr_pow (a, tmp, tmp2, GMP_RNDN);
       mpfr_set_ui (tmp, 9, GMP_RNDN);
       mpfr_mul (a, a, tmp, GMP_RNDN);
       mpfr_set_si (tmp, 44, GMP_RNDN);
       mpfr_mul (tmp, tmp, tmp2, GMP_RNDN);
       mpfr_set_ui (tmp2, 112, GMP_RNDN);
       mpfr_sub (tmp, tmp2, tmp, GMP_RNDN);
       mpfr_add (a, a, tmp, GMP_RNDN);
       mpfr_set_ui (tmp, 121, GMP_RNDN);
       mpfr_div (a, a, tmp, GMP_RNDN);
       mpfr_sqrt (a, a, GMP_RNDN);
     }
     break;
    case 4:
     printf ("Brown, n=%d\n", param);
     {
       mpfr_set_ui (a, 0, GMP_RNDN);
       mpfr_set_ui (tmp, ibase, GMP_RNDN);
       for (i = 1; i <= param; i++) {
	 mpfr_mul (a, a, tmp, GMP_RNDN);
	 mpfr_set_ui (tmp2, i, GMP_RNDN);
	 mpfr_add (a, a, tmp2, GMP_RNDN);
       }
       mpfr_sqrt (a, a, GMP_RNDN);
     }
     break;
    case 5:
     printf ("sqrt(%d)\n", param);
     mpfr_set_ui (a, param, GMP_RNDN);
     mpfr_sqrt (a, a, GMP_RNDN);
     break;
    case 6:
     printf ("Pi\n");
     mpfr_const_pi (a, GMP_RNDN); /* uses Borwein Borwein Plouffe */
     break;
    default:
     printf ("Unknown type of number: %d.\n", which_number);
     usage (argv[0]);
     exit (1);
  }
  
  mpfr_get_str (digits, &exponent, ibase, sx * sy, a, GMP_RNDN);

  /*  printf ("%s\n", digits); */

  for (i = 0; i < sx * sy; i++) { 
    digit = digits[i] - (digits[i] >= 'a' ? 'a' - 10 : '0');
    for (j = 0; j < 3; j++) image[3 * i + j] = 255. * palette[digit][j];
  }

  write_image (sx, sy, image, "zebre.png");

  return 0;
}


void
usage (const char * argv0)
{
  printf ("Usage: %s sx sy type [param [param2]]\n"
	  "sx, sy: ouput size\n"
	  "type: 1: Yéléhada           [2]\n"
	  "      2: 355/113            [0]\n"
	  "      3: Mad zebra (Israël) [1]\n"
	  "      4: Brown              [1]\n"
	  "      5: sqrt(n)            [1]\n"
	  "      6: Pi                 [0]\n"
	  "The numbers in brackets indicate the number of parameters\n",
	  argv0);
}


void
user_error_fn (png_structp png_ptr, png_const_charp error_msg)
{
  fprintf (stderr, "%s\n", error_msg);
  exit (1);
}

void
user_warning_fn (png_structp png_ptr, png_const_charp error_msg)
{
  fprintf (stderr, "%s\n", error_msg);
}

/* Thanks to David Madore for helping me to write this function */
void
write_image (int sx, int sy, unsigned char * image, const char * filename)
{
  FILE * outf = fopen (filename, "wb");
  char errmsg[80];
  png_structp png_ptr;
  png_infop info_ptr;
  png_byte **row_pointers;
  int i;

  if (!outf) {
    snprintf (errmsg, 79, "Couldn't open %s", filename);
    perror (errmsg);
  }
  row_pointers = (png_byte **) malloc (sy * sizeof (png_byte *));
  
  for (i = 0; i < sy; i++ )
    row_pointers[i] = (png_byte *)(image + i * sx * 3);
  png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL,
                                     user_error_fn, user_warning_fn);
  if (!png_ptr) {
    fprintf (stderr, "Can't create png_ptr.\n");
    exit (1);
  }
  info_ptr = png_create_info_struct (png_ptr);
  if (!info_ptr ) {
    fprintf (stderr, "Can't create info_ptr.\n");
    exit (1);
  }
  png_set_IHDR (png_ptr, info_ptr, sx, sy,
		8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
  png_init_io (png_ptr, outf);
  png_write_info (png_ptr, info_ptr);
  png_write_image (png_ptr, row_pointers);
  png_write_end (png_ptr, info_ptr);
  png_write_flush (png_ptr);
  png_destroy_write_struct (&png_ptr, &info_ptr);
  
  fclose (outf);
}