src/utils_math.c

   1 /*
   2  *  Copyright (C) 1998,1999  Helmut Dersch <der@fh-furtwangen.de>
   3  *  Copyright (C) 2007  Florian octo Forster <octo at verplant.org>
   4  *
   5  *  This program is free software; you can redistribute it and/or modify it
   6  *  under the terms of the GNU General Public License as published by the Free
   7  *  Software Foundation; either version 2 of the License, or (at your option)
   8  *  any later version.
   9  *
  10  *  This program is distributed in the hope that it will be useful, but
  11  *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13  *  for more details.
  14  *
  15  *  You should have received a copy of the GNU General Public License along
  16  *  with this program; if not, write to the Free Software Foundation, Inc., 51
  17  *  Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  18  **/
  19
  20 #include <stdlib.h>
  21 #include <stdio.h>
  22 #include <math.h>
  23
  24 #include "utils_math.h"
  25
  26 // Lookup Tables for Trig-functions and interpolator
  27
  28 #define ATAN_TABLE_SIZE 4096
  29 #define SQRT_TABLE_SIZE 4096
  30
  31 static int *sqrt_table = NULL;
  32
  33 static int setup_sqrt_table (void)
  34 {
  35   double z;
  36   int i;
  37
  38   sqrt_table = (int *) malloc (SQRT_TABLE_SIZE * sizeof (int));
  39   if (sqrt_table == NULL)
  40     return (-1);
  41
  42   /* Calculate the table entries. */
  43   for (i = 0; i < SQRT_TABLE_SIZE; i++)
  44   {
  45     z = ((double) i) / ((double) (SQRT_TABLE_SIZE - 1));
  46     sqrt_table[i] = (int) (sqrt( 1.0 + z*z ) * 4096.0 + .5);
  47   }
  48
  49   return (0);
  50 } /* int setup_sqrt_table */
  51
  52 /*
  53  * Calculate
  54  *   256 * sqrt(x1^2 + x2^2)
  55  * efficiently
  56  */
  57 int utils_sqrt2 (int x1, int x2)
  58 {
  59   int ret;
  60
  61   if (sqrt_table == NULL)
  62     if (setup_sqrt_table () != 0)
  63       return (-1);
  64
  65   if ((x1 == 0) && (x2 == 0))
  66   {
  67     ret = 0;
  68   }
  69   else if (x1 > x2)
  70   {
  71     ret = x1 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x2 / x1] / 16;
  72   }
  73   else /* (x2 < x1) */
  74   {
  75     ret = x2 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x1 / x2] / 16;
  76   }
  77
  78   return (ret);
  79 } /* int utils_sqrt2 */
  80
  81 /*
  82  * vim: set tabstop=8 softtabstop=2 shiftwidth=2 :
  83  */