Mozzi  alpha 0.01.1n
sound synthesis library for Arduino
 All Classes Functions
fastSqrt.cpp
00001 /*
00002 // Integer Square Root function
00003 // Contributors include Arne Steinarson for the basic approximation idea,
00004 // Dann Corbit and Mathew Hendry for the first cut at the algorithm,
00005 // Lawrence Kirby for the rearrangement, improvments and range optimization
00006 // and Paul Hsieh for the round-then-adjust idea.
00007 */
00008 
00009 #ifndef FASTSQRT_H_
00010 #define FASTSQRT_H_
00011 
00012 static unsigned fastSqrt(unsigned long x)
00013 {
00014                 static const unsigned char sqq_table[] =
00015                         {
00016                                 0,  16,  22,  27,  32,  35,  39,  42,  45,  48,  50,  53,  55,  57,
00017                                 59,  61,  64,  65,  67,  69,  71,  73,  75,  76,  78,  80,  81,  83,
00018                                 84,  86,  87,  89,  90,  91,  93,  94,  96,  97,  98,  99, 101, 102,
00019                                 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
00020                                 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
00021                                 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
00022                                 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
00023                                 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
00024                                 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
00025                                 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
00026                                 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
00027                                 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
00028                                 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
00029                                 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
00030                                 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
00031                                 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
00032                                 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
00033                                 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
00034                                 253, 254, 254, 255
00035                         };
00036 
00037                 unsigned long xn;
00038 
00039                 if (x >= 0x10000)
00040                                 if (x >= 0x1000000)
00041                                                 if (x >= 0x10000000)
00042                                                                 if (x >= 0x40000000)
00043                                                                 {
00044                                                                                 if (x >= 65535UL*65535UL)
00045                                                                                                 return 65535;
00046                                                                                 xn = sqq_table[x>>24] << 8;
00047                                                                 }
00048                                                                 else
00049                                                                                 xn = sqq_table[x>>22] << 7;
00050                                                 else if (x >= 0x4000000)
00051                                                                 xn = sqq_table[x>>20] << 6;
00052                                                 else
00053                                                                 xn = sqq_table[x>>18] << 5;
00054                                 else
00055                                 {
00056                                                 if (x >= 0x100000)
00057                                                                 if (x >= 0x400000)
00058                                                                                 xn = sqq_table[x>>16] << 4;
00059                                                                 else
00060                                                                                 xn = sqq_table[x>>14] << 3;
00061                                                 else if (x >= 0x40000)
00062                                                                 xn = sqq_table[x>>12] << 2;
00063                                                 else
00064                                                                 xn = sqq_table[x>>10] << 1;
00065 
00066                                                 goto nr1;
00067                                 }
00068                 else if (x >= 0x100)
00069                 {
00070                                 if (x >= 0x1000)
00071                                                 if (x >= 0x4000)
00072                                                                 xn = (sqq_table[x>>8] >> 0) + 1;
00073                                                 else
00074                                                                 xn = (sqq_table[x>>6] >> 1) + 1;
00075                                 else if (x >= 0x400)
00076                                                 xn = (sqq_table[x>>4] >> 2) + 1;
00077                                 else
00078                                                 xn = (sqq_table[x>>2] >> 3) + 1;
00079 
00080                                 goto adj;
00081                 }
00082                 else
00083                                 return sqq_table[x] >> 4;
00084 
00085                 /* Run two iterations of the standard convergence formula */
00086 
00087                 xn = (xn + 1 + x / xn) / 2;
00088 nr1:
00089                 xn = (xn + 1 + x / xn) / 2;
00090 adj:
00091 
00092                 if (xn * xn > x) /* Correct rounding if necessary */
00093                                 xn--;
00094 
00095                 return xn;
00096 }
00097 
00098 
00099 #endif /* FASTSQRT_H_ */