modules/base/math/fixed_point/f32_sqrt.c

   1 /*
   2  *  Copyright Droids Corporation, Microb Technology, Eirbot (2005)
   3  *
   4  *  This program is free software; you can redistribute it and/or modify
   5  *  it under the terms of the GNU General Public License as published by
   6  *  the Free Software Foundation; either version 2 of the License, or
   7  *  (at your option) any later version.
   8  *
   9  *  This program is distributed in the hope that it will be useful,
  10  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  *  GNU General Public License for more details.
  13  *
  14  *  You should have received a copy of the GNU General Public License
  15  *  along with this program; if not, write to the Free Software
  16  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  17  *
  18  *  Revision : $Id: f32_sqrt.c,v 1.5.4.4 2009-02-02 22:21:20 zer0 Exp $
  19  *
  20  */
  21
  22 #include <stdio.h>
  23 #include <f32.h>
  24 #include <s32_to_f32.h>
  25
  26 #define NEXT(n, i)  (((n) + (i)/(n)) >> 1)
  27
  28 static uint32_t u32_sqrt(uint32_t number)
  29 {
  30         uint32_t n  = 1;
  31         uint32_t n1 = NEXT(n, number);
  32
  33         if (number == 0)
  34                 return 0;
  35
  36         while(ABS(n1 - n) > 1) {
  37                 n  = n1;
  38                 n1 = NEXT(n, number);
  39         }
  40         while((n1*n1) > number) {
  41                 n1 -= 1;
  42         }
  43         return n1;
  44 }
  45
  46 f32 f32_sqrt(f32 f)
  47 {
  48         uint32_t a,b,c,d;
  49
  50         if (F32_IS_ZERO(f))
  51                 return F32_ZERO;
  52
  53         if (F32_IS_NEG(f))
  54                 return F32_NAN;
  55
  56         if(f.f32_integer) {
  57                 /* sqrt(a+b) = sqrt(a)*sqrt(1+b/a) */
  58                 a=(uint32_t)(f.f32_integer) << 16 ;
  59                 b=(uint32_t)(f.f32_decimal) << 16 ;
  60                 c=u32_sqrt(a);
  61                 d=u32_sqrt(0x10000 + (b/a));
  62                 return f32_mul(s32_to_f32( c<<8 ), s32_to_f32( d<<8 ));
  63         }
  64         else {
  65                 b=(uint32_t)(f.f32_decimal) << 16 ;
  66                 return s32_to_f32(u32_sqrt(b));
  67         }
  68 }
  69