1:   #include <windows.h>
2:   #include <stdio.h>
3:   #include <math.h>
4:   #include <unistd.h>
5:   #include <limits.h>
6:   #include <time.h>
7:  
8:   #undef NDEBUG
9:   #include <assert.h>
10:   #include <string.h>
11:   #include <mmsystem.h>
12:   #include <float.h>
13:   #include <stdlib.h>
14:  
15:  
16:   /*
17:   NVidia corporation provides sample code at their website 
18:   for *their* fast square root implementation. 
19:   It is basically a 256K table look up. 
20:  
21:   This is the link to the fastsqrt function source
22:   http://www.azillionmonkeys.com/qed/fastmath.cpp
23:  
24:   */
25:   float fastsqrt (float n);
26:  
27:  
28:   /*
29:  
30:   The isqrt function is a integer square root implementation
31:    , it was retrieved from http://www.azillionmonkeys.com/qed/sqroot.html
32:  
33:   It is fast but not faster than the table lookup Via implementation
34:  
35:   Since it provides only integer square root and is slower than
36:   the table lookup implementation , there is no reason to use this method
37:  
38:   Use it only if you cannot use a table lookup implementation
39:  
40:   */
41:   unsigned int isqrt (unsigned int n);
42:  
43:   double
44:   sqrt2 (double val)
45:   {
46:  
47:  
48:      register double temp2 0;
49:      register double temp3;
50:      register double lastval;
51:  
52:  
53:      if (== val)
54:        {
55:       return 1;
56:        }
57:  
58:      lastval = -val;
59:  
60:  
61:   // Requires the precomputed sqrt table in the fast sqrt Via implementation
62:  
63:   // It gives a very good approximation of the sqrt using 32 bit floating point
64:  
65:   // There is no reason to compute something that you can keep in a precomputed table
66:      temp3 fastsqrt (val);
67:  
68:   // if you prefer to use as approximation the fast integer square root 
69:      // temp3 = (double) isqrt ((unsigned int) val);
70:  
71:  
72:      while (1)
73:        {
74:  
75:   //Now with a 32 bit floating point square root we only 
76:   //need to extend it to a 64 bit floating point implementation
77:  
78:   //This is the slower point of the code , since it requires a division
79:   //and a multiplication
80:  
81:   //It spend 20 percent of the time in the table lookup and the
82:   //rest of the time to make the double type approximation
83:  
84:   // Notice that yet I think that there is too much room to
85:   // improve the square root algorithm
86:  
87:       if (temp2 == temp3)
88:         {
89:  
90:            goto achou;
91:  
92:         }
93:  
94:       temp2 = (val temp3);
95:  
96:       temp3 = (temp2 temp3) * 0.5;
97:  
98:       if (lastval == temp3)
99:         {
100:  
101:            goto achou;
102:  
103:         }
104:  
105:       lastval temp3;
106:  
107:        }
108:      ;
109:    achou:
110:      ;
111:  
112:  
113:      return temp2;
114:  
115:   }
116:  
117:  
118:   /*
119:   this is the source to the integer square root implementation
120:    , it was retrieved from http://www.azillionmonkeys.com/qed/sqroot.html
121:  
122:   It is fast but not faster than the table lookup Via implementation
123:  
124:   Since it provides only integer square root and is slower than
125:   the table lookup implementation , there is no reason to use this
126:  
127:   Use it only if you cannot use a table lookup implementation
128:   */
129:  
130:  
131:   unsigned int
132:   isqrt (unsigned int n)
133:   {
134:      unsigned int root 0try;
135:  
136:      try root + (<< (15));
137:  
138:      if (>= try << (15))
139:        {
140:       -= try << (15);
141:       root |= << (15);
142:        };
143:  
144:      try root + (<< (14));
145:  
146:      if (>= try << (14))
147:        {
148:       -= try << (14);
149:       root |= << (14);
150:        };
151:  
152:      try root + (<< (13));
153:  
154:      if (>= try << (13))
155:        {
156:       -= try << (13);
157:       root |= << (13);
158:        };
159:  
160:      try root + (<< (12));
161:  
162:      if (>= try << (12))
163:        {
164:       -= try << (12);
165:       root |= << (12);
166:        };
167:  
168:      try root + (<< (11));
169:  
170:      if (>= try << (11))
171:        {
172:       -= try << (11);
173:       root |= << (11);
174:        };
175:  
176:      try root + (<< (10));
177:  
178:      if (>= try << (10))
179:        {
180:       -= try << (10);
181:       root |= << (10);
182:        };
183:  
184:      try root + (<< (9));
185:  
186:      if (>= try << (9))
187:        {
188:       -= try << (9);
189:       root |= << (9);
190:        };
191:  
192:      try root + (<< (8));
193:  
194:      if (>= try << (8))
195:        {
196:       -= try << (8);
197:       root |= << (8);
198:        };
199:  
200:      try root + (<< (7));
201:  
202:      if (>= try << (7))
203:        {
204:       -= try << (7);
205:       root |= << (7);
206:        };
207:  
208:      try root + (<< (6));
209:  
210:      if (>= try << (6))
211:        {
212:       -= try << (6);
213:       root |= << (6);
214:        };
215:  
216:      try root + (<< (5));
217:  
218:      if (>= try << (5))
219:        {
220:       -= try << (5);
221:       root |= << (5);
222:        };
223:  
224:      try root + (<< (4));
225:  
226:      if (>= try << (4))
227:        {
228:       -= try << (4);
229:       root |= << (4);
230:        };
231:  
232:      try root + (<< (3));
233:  
234:      if (>= try << (3))
235:        {
236:       -= try << (3);
237:       root |= << (3);
238:        };
239:  
240:      try root + (<< (2));
241:  
242:      if (>= try << (2))
243:        {
244:       -= try << (2);
245:       root |= << (2);
246:        };
247:  
248:      try root + (<< (1));
249:  
250:      if (>= try << (1))
251:        {
252:       -= try << (1);
253:       root |= << (1);
254:        };
255:  
256:      try root + (<< (0));
257:  
258:      if (>= try << (0))
259:        {
260:       -= try << (0);
261:       root |= << (0);
262:        };
263:  
264:      return root >> 1;
265:   }
266:  
267:  
268:   int
269:   main ()
270:   {
271:   #define dprintf printf
272:      int i;
273:      int p;
274:  
275:      volatile double sqrtval 0;
276:      //to avoid the optimizer to remove this reference
277:  
278:  
279:      // this call will initiate the table on the fastsqrt function
280:      build_sqrt_table ();
281:  
282:      GetTickCount ();
283:  
284:      for (130000000i++)
285:        {
286:       sqrtval sqrt (i);
287:        }
288:  
289:      GetTickCount () - p;
290:      dprintf ("sqrt       = %5d milliseconds %f\n"psqrtval);
291:  
292:      GetTickCount ();
293:  
294:      for (130000000i++)
295:        {
296:       sqrtval sqrt2 (i);
297:        }
298:  
299:      GetTickCount () - p;
300:      dprintf ("sqrt2      = %5d milliseconds %f\n"psqrtval);
301:  
302:      GetTickCount ();
303:  
304:      for (130000000i++)
305:        {
306:       sqrtval fastsqrt (i);
307:        }
308:  
309:      GetTickCount () - p;
310:      dprintf ("fastsqrt   = %5d milliseconds %f\n"psqrtval);
311:  
312:      GetTickCount ();
313:  
314:      for (130000000i++)
315:        {
316:       sqrtval = (doubleisqrt ((inti);
317:        }
318:  
319:      GetTickCount () - p;
320:      dprintf ("isqrt      = %5d milliseconds %f\n"psqrtval);
321:  
322:      return 0;
323:  
324:  
325:   }