1  /*
       2   * Copyright (c) 2008-2020 Stefan Krah. All rights reserved.
       3   *
       4   * Redistribution and use in source and binary forms, with or without
       5   * modification, are permitted provided that the following conditions
       6   * are met:
       7   *
       8   * 1. Redistributions of source code must retain the above copyright
       9   *    notice, this list of conditions and the following disclaimer.
      10   *
      11   * 2. Redistributions in binary form must reproduce the above copyright
      12   *    notice, this list of conditions and the following disclaimer in the
      13   *    documentation and/or other materials provided with the distribution.
      14   *
      15   * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
      16   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
      17   * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
      18   * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
      19   * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
      20   * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
      21   * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
      22   * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
      23   * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
      24   * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
      25   * SUCH DAMAGE.
      26   */
      27  
      28  
      29  #include "mpdecimal.h"
      30  
      31  #include <stdint.h>
      32  #include <stdio.h>
      33  #include <stdlib.h>
      34  #include <time.h>
      35  
      36  
      37  static void
      38  err_exit(const char *msg)
      39  {
      40      fprintf(stderr, "%s\n", msg);
      41      exit(1);
      42  }
      43  
      44  static mpd_t *
      45  new_mpd(void)
      46  {
      47      mpd_t *x = mpd_qnew();
      48      if (x == NULL) {
      49          err_exit("out of memory");
      50      }
      51  
      52      return x;
      53  }
      54  
      55  /* Nonsense version of escape-time algorithm for calculating a mandelbrot
      56   * set. Just for benchmarking. */
      57  static void
      58  color_point(mpd_t *x0, mpd_t *y0, long maxiter, mpd_context_t *ctx)
      59  {
      60      mpd_t *x, *y, *sq_x, *sq_y;
      61      mpd_t *two;
      62  
      63      x = new_mpd();
      64      y = new_mpd();
      65      mpd_set_u32(x, 0, ctx);
      66      mpd_set_u32(y, 0, ctx);
      67  
      68      sq_x = new_mpd();
      69      sq_y = new_mpd();
      70      mpd_set_u32(sq_x, 0, ctx);
      71      mpd_set_u32(sq_y, 0, ctx);
      72  
      73      two = new_mpd();
      74      mpd_set_u32(two, 2, ctx);
      75  
      76      for (long i = 0; i < maxiter; i++) {
      77          mpd_mul(y, x, y, ctx);
      78          mpd_mul(y, y, two, ctx);
      79          mpd_add(y, y, y0, ctx);
      80  
      81          mpd_sub(x, sq_x, sq_y, ctx);
      82          mpd_add(x, x, x0, ctx);
      83  
      84          mpd_mul(sq_x, x, x, ctx);
      85          mpd_mul(sq_y, y, y, ctx);
      86      }
      87  
      88      mpd_copy(x0, x, ctx);
      89  
      90      mpd_del(two);
      91      mpd_del(sq_y);
      92      mpd_del(sq_x);
      93      mpd_del(y);
      94      mpd_del(x);
      95  }
      96  
      97  
      98  int
      99  main(int argc, char **argv)
     100  {
     101      mpd_context_t ctx;
     102      mpd_t *x0, *y0;
     103      uint32_t prec = 19;
     104      long iter = 10000000;
     105      clock_t start_clock, end_clock;
     106  
     107      if (argc != 3) {
     108          err_exit("usage: bench prec iter\n");
     109      }
     110      prec = strtoul(argv[1], NULL, 10);
     111      iter = strtol(argv[2], NULL, 10);
     112  
     113      mpd_init(&ctx, prec);
     114      /* no more MPD_MINALLOC changes after here */
     115  
     116      x0 = new_mpd();
     117      y0 = new_mpd();
     118      mpd_set_string(x0, "0.222", &ctx);
     119      mpd_set_string(y0, "0.333", &ctx);
     120      if (ctx.status & MPD_Errors) {
     121          mpd_del(y0);
     122          mpd_del(x0);
     123          err_exit("unexpected error during conversion");
     124      }
     125  
     126      start_clock = clock();
     127      color_point(x0, y0, iter, &ctx);
     128      end_clock = clock();
     129  
     130      mpd_print(x0);
     131      fprintf(stderr, "time: %f\n\n", (double)(end_clock-start_clock)/(double)CLOCKS_PER_SEC);
     132  
     133      mpd_del(y0);
     134      mpd_del(x0);
     135  
     136      return 0;
     137  }