sys/src/libc/power/sqrt.c

   1 #include <u.h>
   2 #include <libc.h>
   3
   4 static  long    sqtab[64] =
   5 {
   6         0x6cdb2, 0x726d4, 0x77ea3, 0x7d52f, 0x82a85, 0x87eb1, 0x8d1c0, 0x923bd,
   7         0x974b2, 0x9c4a8, 0xa13a9, 0xa61be, 0xaaeee, 0xafb41, 0xb46bf, 0xb916e,
   8         0xbdb55, 0xc247a, 0xc6ce3, 0xcb495, 0xcfb95, 0xd41ea, 0xd8796, 0xdcca0,
   9         0xe110c, 0xe54dd, 0xe9818, 0xedac0, 0xf1cd9, 0xf5e67, 0xf9f6e, 0xfdfef,
  10         0x01fe0, 0x05ee6, 0x09cfd, 0x0da30, 0x11687, 0x1520c, 0x18cc8, 0x1c6c1,
  11         0x20000, 0x2388a, 0x27068, 0x2a79e, 0x2de32, 0x3142b, 0x3498c, 0x37e5b,
  12         0x3b29d, 0x3e655, 0x41989, 0x44c3b, 0x47e70, 0x4b02b, 0x4e16f, 0x51241,
  13         0x542a2, 0x57296, 0x5a220, 0x5d142, 0x60000, 0x62e5a, 0x65c55, 0x689f2,
  14 };
  15
  16 double
  17 sqrt(double arg)
  18 {
  19         int e, ms;
  20         double a, t;
  21         union
  22         {
  23                 double  d;
  24                 struct
  25                 {
  26                         long    ms;
  27                         long    ls;
  28                 };
  29         } u;
  30
  31         u.d = arg;
  32         ms = u.ms;
  33
  34         /*
  35          * sign extend the mantissa with
  36          * exponent. result should be > 0 for
  37          * normal case.
  38          */
  39         e = ms >> 20;
  40         if(e <= 0) {
  41                 if(e == 0)
  42                         return 0;
  43                 return NaN();
  44         }
  45
  46         /*
  47          * pick up arg/4 by adjusting exponent
  48          */
  49         u.ms = ms - (2 << 20);
  50         a = u.d;
  51
  52         /*
  53          * use 5 bits of mantissa and 1 bit
  54          * of exponent to form table index.
  55          * insert exponent/2 - 1.
  56          */
  57         e = (((e - 1023) >> 1) + 1022) << 20;
  58         u.ms = *(long*)((char*)sqtab + ((ms >> 13) & 0xfc)) | e;
  59         u.ls = 0;
  60
  61         /*
  62          * three laps of newton
  63          */
  64         e = 1 << 20;
  65         t = u.d;
  66         u.d = t + a/t;
  67         u.ms -= e;              /* u.d /= 2; */
  68         t = u.d;
  69         u.d = t + a/t;
  70         u.ms -= e;              /* u.d /= 2; */
  71         t = u.d;
  72
  73         return t + a/t;
  74 }
  75
  76 /*
  77  * this is the program that generated the table.
  78  * it calls sqrt by some other means.
  79  *
  80  * void
  81  * main(void)
  82  * {
  83  *      int i;
  84  *      union   U
  85  *      {
  86  *              double  d;
  87  *              struct
  88  *              {
  89  *                      long    ms;
  90  *                      long    ls;
  91  *              };
  92  *      } u;
  93  *
  94  *      for(i=0; i<64; i++) {
  95  *              u.ms = (i<<15) | 0x3fe04000;
  96  *              u.ls = 0;
  97  *              u.d = sqrt(u.d);
  98  *              print(" 0x%.5lux,", u.ms & 0xfffff);
  99  *      }
 100  *      print("\n");
 101  *      exits(0);
 102  * }
 103  */