lemon-project-template-glpk: 4fc6ad2fb8a6 deps/glpk/src/glplib02.c

lemon-project-template-glpk

view deps/glpk/src/glplib02.c @ 11:4fc6ad2fb8a6

Test GLPK in src/main.cc

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 21:43:29 +0100
parents
children

line source

1 /* glplib02.c (64-bit arithmetic) */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

9 * E-mail: <mao@gnu.org>.

10 *

11 * GLPK is free software: you can redistribute it and/or modify it

12 * under the terms of the GNU General Public License as published by

13 * the Free Software Foundation, either version 3 of the License, or

14 * (at your option) any later version.

15 *

16 * GLPK is distributed in the hope that it will be useful, but WITHOUT

17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

19 * License for more details.

20 *

21 * You should have received a copy of the GNU General Public License

22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.

23 ***********************************************************************/

25 #include "glpenv.h"

26 #include "glplib.h"

28 /***********************************************************************

29 * NAME

30 *

31 * xlset - expand integer to long integer

32 *

33 * SYNOPSIS

34 *

35 * #include "glplib.h"

36 * glp_long xlset(int x);

37 *

38 * RETURNS

39 *

40 * The routine xlset returns x expanded to long integer. */

42 glp_long xlset(int x)

43 { glp_long t;

44 t.lo = x, t.hi = (x >= 0 ? 0 : -1);

45 return t;

46 }

48 /***********************************************************************

49 * NAME

50 *

51 * xlneg - negate long integer

52 *

53 * SYNOPSIS

54 *

55 * #include "glplib.h"

56 * glp_long xlneg(glp_long x);

57 *

58 * RETURNS

59 *

60 * The routine xlneg returns the difference 0 - x. */

62 glp_long xlneg(glp_long x)

63 { if (x.lo)

64 x.lo = - x.lo, x.hi = ~x.hi;

65 else

66 x.hi = - x.hi;

67 return x;

68 }

70 /***********************************************************************

71 * NAME

72 *

73 * xladd - add long integers

74 *

75 * SYNOPSIS

76 *

77 * #include "glplib.h"

78 * glp_long xladd(glp_long x, glp_long y);

79 *

80 * RETURNS

81 *

82 * The routine xladd returns the sum x + y. */

84 glp_long xladd(glp_long x, glp_long y)

85 { if ((unsigned int)x.lo <= 0xFFFFFFFF - (unsigned int)y.lo)

86 x.lo += y.lo, x.hi += y.hi;

87 else

88 x.lo += y.lo, x.hi += y.hi + 1;

89 return x;

90 }

92 /***********************************************************************

93 * NAME

94 *

95 * xlsub - subtract long integers

96 *

97 * SYNOPSIS

98 *

99 * #include "glplib.h"

100 * glp_long xlsub(glp_long x, glp_long y);

101 *

102 * RETURNS

103 *

104 * The routine xlsub returns the difference x - y. */

105

106 glp_long xlsub(glp_long x, glp_long y)

107 { return

108 xladd(x, xlneg(y));

109 }

110

111 /***********************************************************************

112 * NAME

113 *

114 * xlcmp - compare long integers

115 *

116 * SYNOPSIS

117 *

118 * #include "glplib.h"

119 * int xlcmp(glp_long x, glp_long y);

120 *

121 * RETURNS

122 *

123 * The routine xlcmp returns the sign of the difference x - y. */

124

125 int xlcmp(glp_long x, glp_long y)

126 { if (x.hi >= 0 && y.hi < 0) return +1;

127 if (x.hi < 0 && y.hi >= 0) return -1;

128 if ((unsigned int)x.hi < (unsigned int)y.hi) return -1;

129 if ((unsigned int)x.hi > (unsigned int)y.hi) return +1;

130 if ((unsigned int)x.lo < (unsigned int)y.lo) return -1;

131 if ((unsigned int)x.lo > (unsigned int)y.lo) return +1;

132 return 0;

133 }

134

135 /***********************************************************************

136 * NAME

137 *

138 * xlmul - multiply long integers

139 *

140 * SYNOPSIS

141 *

142 * #include "glplib.h"

143 * glp_long xlmul(glp_long x, glp_long y);

144 *

145 * RETURNS

146 *

147 * The routine xlmul returns the product x * y. */

148

149 glp_long xlmul(glp_long x, glp_long y)

150 { unsigned short xx[8], yy[4];

151 xx[4] = (unsigned short)x.lo;

152 xx[5] = (unsigned short)(x.lo >> 16);

153 xx[6] = (unsigned short)x.hi;

154 xx[7] = (unsigned short)(x.hi >> 16);

155 yy[0] = (unsigned short)y.lo;

156 yy[1] = (unsigned short)(y.lo >> 16);

157 yy[2] = (unsigned short)y.hi;

158 yy[3] = (unsigned short)(y.hi >> 16);

159 bigmul(4, 4, xx, yy);

160 x.lo = (unsigned int)xx[0] | ((unsigned int)xx[1] << 16);

161 x.hi = (unsigned int)xx[2] | ((unsigned int)xx[3] << 16);

162 return x;

163 }

164

165 /***********************************************************************

166 * NAME

167 *

168 * xldiv - divide long integers

169 *

170 * SYNOPSIS

171 *

172 * #include "glplib.h"

173 * glp_ldiv xldiv(glp_long x, glp_long y);

174 *

175 * RETURNS

176 *

177 * The routine xldiv returns a structure of type glp_ldiv containing

178 * members quot (the quotient) and rem (the remainder), both of type

179 * glp_long. */

180

181 glp_ldiv xldiv(glp_long x, glp_long y)

182 { glp_ldiv t;

183 int m, sx, sy;

184 unsigned short xx[8], yy[4];

185 /* sx := sign(x) */

186 sx = (x.hi < 0);

187 /* sy := sign(y) */

188 sy = (y.hi < 0);

189 /* x := |x| */

190 if (sx) x = xlneg(x);

191 /* y := |y| */

192 if (sy) y = xlneg(y);

193 /* compute x div y and x mod y */

194 xx[0] = (unsigned short)x.lo;

195 xx[1] = (unsigned short)(x.lo >> 16);

196 xx[2] = (unsigned short)x.hi;

197 xx[3] = (unsigned short)(x.hi >> 16);

198 yy[0] = (unsigned short)y.lo;

199 yy[1] = (unsigned short)(y.lo >> 16);

200 yy[2] = (unsigned short)y.hi;

201 yy[3] = (unsigned short)(y.hi >> 16);

202 if (yy[3])

203 m = 4;

204 else if (yy[2])

205 m = 3;

206 else if (yy[1])

207 m = 2;

208 else if (yy[0])

209 m = 1;

210 else

211 xerror("xldiv: divide by zero\n");

212 bigdiv(4 - m, m, xx, yy);

213 /* remainder in x[0], x[1], ..., x[m-1] */

214 t.rem.lo = (unsigned int)xx[0], t.rem.hi = 0;

215 if (m >= 2) t.rem.lo |= (unsigned int)xx[1] << 16;

216 if (m >= 3) t.rem.hi = (unsigned int)xx[2];

217 if (m >= 4) t.rem.hi |= (unsigned int)xx[3] << 16;

218 if (sx) t.rem = xlneg(t.rem);

219 /* quotient in x[m], x[m+1], ..., x[4] */

220 t.quot.lo = (unsigned int)xx[m], t.quot.hi = 0;

221 if (m <= 3) t.quot.lo |= (unsigned int)xx[m+1] << 16;

222 if (m <= 2) t.quot.hi = (unsigned int)xx[m+2];

223 if (m <= 1) t.quot.hi |= (unsigned int)xx[m+3] << 16;

224 if (sx ^ sy) t.quot = xlneg(t.quot);

225 return t;

226 }

227

228 /***********************************************************************

229 * NAME

230 *

231 * xltod - convert long integer to double

232 *

233 * SYNOPSIS

234 *

235 * #include "glplib.h"

236 * double xltod(glp_long x);

237 *

238 * RETURNS

239 *

240 * The routine xltod returns x converted to double. */

241

242 double xltod(glp_long x)

243 { double s, z;

244 if (x.hi >= 0)

245 s = +1.0;

246 else

247 s = -1.0, x = xlneg(x);

248 if (x.hi >= 0)

249 z = 4294967296.0 * (double)x.hi + (double)(unsigned int)x.lo;

250 else

251 { xassert(x.hi == 0x80000000 && x.lo == 0x00000000);

252 z = 9223372036854775808.0; /* 2^63 */

253 }

254 return s * z;

255 }

256

257 char *xltoa(glp_long x, char *s)

258 { /* convert long integer to character string */

259 static const char *d = "0123456789";

260 glp_ldiv t;

261 int neg, len;

262 if (x.hi >= 0)

263 neg = 0;

264 else

265 neg = 1, x = xlneg(x);

266 if (x.hi >= 0)

267 { len = 0;

268 while (!(x.hi == 0 && x.lo == 0))

269 { t = xldiv(x, xlset(10));

270 xassert(0 <= t.rem.lo && t.rem.lo <= 9);

271 s[len++] = d[t.rem.lo];

272 x = t.quot;

273 }

274 if (len == 0) s[len++] = d[0];

275 if (neg) s[len++] = '-';

276 s[len] = '\0';

277 strrev(s);

278 }

279 else

280 strcpy(s, "-9223372036854775808"); /* -2^63 */

281 return s;

282 }

283

284 /**********************************************************************/

285

286 #if 0

287 #include "glprng.h"

288

289 #define N_TEST 1000000

290 /* number of tests */

291

292 static glp_long myrand(RNG *rand)

293 { glp_long x;

294 int k;

295 k = rng_unif_rand(rand, 4);

296 xassert(0 <= k && k <= 3);

297 x.lo = rng_unif_rand(rand, 65536);

298 if (k == 1 || k == 3)

299 { x.lo <<= 16;

300 x.lo += rng_unif_rand(rand, 65536);

301 }

302 if (k <= 1)

303 x.hi = 0;

304 else

305 x.hi = rng_unif_rand(rand, 65536);

306 if (k == 3)

307 { x.hi <<= 16;

308 x.hi += rng_unif_rand(rand, 65536);

309 }

310 if (rng_unif_rand(rand, 2)) x = xlneg(x);

311 return x;

312 }

313

314 int main(void)

315 { RNG *rand;

316 glp_long x, y;

317 glp_ldiv z;

318 int test;

319 rand = rng_create_rand();

320 for (test = 1; test <= N_TEST; test++)

321 { x = myrand(rand);

322 y = myrand(rand);

323 if (y.lo == 0 && y.hi == 0) y.lo = 1;

324 /* z.quot := x div y, z.rem := x mod y */

325 z = xldiv(x, y);

326 /* x must be equal to y * z.quot + z.rem */

327 xassert(xlcmp(x, xladd(xlmul(y, z.quot), z.rem)) == 0);

328 }

329 xprintf("%d tests successfully passed\n", N_TEST);

330 rng_delete_rand(rand);

331 return 0;

332 }

333 #endif

334

335 /* eof */