1 /* -*- tab-width:4; -*- */
2 /*
3 * Numbers
4 */
5 #include "s.h"
6
7 #define SMALLNUMP(x) (fabs(x) <= SOBJ_INUM_MAX)
8 #define ISINT(x) (floor(x) == (x))
9
scm_newbnum()10 SOBJ scm_newbnum()
11 {
12 SOBJ z = scm_newcell(SOBJ_T_BNUM);
13 SCM_BNUM(z) = scm_must_alloc(sizeof(MP_INT));
14 mpz_init(SCM_BNUM(z));
15 return z;
16 }
17
scm_bnum_sweep(SOBJ obj)18 void scm_bnum_sweep(SOBJ obj)
19 {
20 fprintf(stderr, "freeing bignum at %p\n", obj);
21 mpz_clear(SCM_BNUM(obj));
22 scm_free(SCM_BNUM(obj));
23 }
24
25
scm_mkinum(long x)26 SOBJ scm_mkinum(long x)
27 {
28 return(SCM_MKINUM(x));
29 }
30
scm_mkfnum(double x)31 SOBJ scm_mkfnum(double x)
32 {
33 SOBJ z = scm_newcell(SOBJ_T_FNUM);
34 SCM_FNUM(z) = x;
35 return z;
36 }
37
38 /*-- conversion : The idea is to convert numbers if they are not of
39 * the same class to an apropriate class suitable for future
40 * operations */
41
42 /*-- convert to big */
scm_int2bnum(long n)43 SOBJ scm_int2bnum(long n)
44 {
45 SOBJ z = scm_newcell(SOBJ_T_BNUM);
46 SCM_BNUM(z) = scm_must_alloc(sizeof(MP_INT));
47 mpz_init_set_si(SCM_BNUM(z), n);
48 return z;
49 }
50
scm_uint2bnum(unsigned long n)51 SOBJ scm_uint2bnum(unsigned long n)
52 {
53 SOBJ z = scm_newcell(SOBJ_T_BNUM);
54 SCM_BNUM(z) = scm_must_alloc(sizeof(MP_INT));
55 mpz_init_set_ui(SCM_BNUM(z), n);
56 return z;
57 }
58
flt2bnum(double n)59 static SOBJ flt2bnum(double n)
60 {
61 SOBJ z = scm_newcell(SOBJ_T_BNUM);
62 SCM_BNUM(z) = scm_must_alloc(sizeof(MP_INT));
63 mpz_init_set_d(SCM_BNUM(z), n);
64 return(z);
65 }
66
67 #ifdef NOT_USED
cstr2bnum(char * s)68 static SOBJ cstr2bnum(char *s)
69 {
70 SOBJ z = scm_newcell(SOBJ_T_BNUM);
71 SCM_BNUM(z) = scm_must_alloc(sizeof(MP_INT));
72 mpz_init_set_str(SCM_BNUM(z), s, 10);
73 return z;
74 }
75 #endif
76
77 /*-- convert to flt */
int2fnum(long n)78 static SOBJ int2fnum(long n)
79 {
80 SOBJ z = scm_newcell(SOBJ_T_FNUM);
81 SCM_FNUM(z) = (double) n;
82 return(z);
83 }
84
big2fnum(MP_INT * n)85 static SOBJ big2fnum(MP_INT *n)
86 {
87 SOBJ z = scm_newcell(SOBJ_T_FNUM);
88 SCM_FNUM(z) = mpz_get_d(n);
89 return(z);
90 }
91
scm_flt2num(double n)92 SOBJ scm_flt2num(double n)
93 {
94 if (SCM_INUM_RANGE(n))
95 return(SCM_MKINUM((long) n));
96
97 return(flt2bnum(n));
98 }
99
scm_int2num(long n)100 SOBJ scm_int2num(long n)
101 {
102 if (SCM_INUM_RANGE(n))
103 return(SCM_MKINUM(n));
104 return(scm_int2bnum(n));
105 }
106
scm_uint2num(unsigned long n)107 SOBJ scm_uint2num(unsigned long n)
108 {
109 if (n <= SOBJ_INUM_MAX)
110 return(SCM_MKINUM(n));
111
112 return(scm_uint2bnum(n));
113 }
114
scm_number2double(SOBJ x)115 double scm_number2double(SOBJ x)
116 {
117 if (!SCM_NUMBERP(x)) SCM_ERR("scm_number2double: bad number", x);
118 switch(SCM_OBJTYPE(x)) {
119 case SOBJ_T_INUM: return( (double) SCM_INUM(x) );
120 case SOBJ_T_FNUM: return( SCM_FNUM(x) );
121 case SOBJ_T_BNUM: return( mpz_get_d(SCM_BNUM(x)) );
122 }
123 return(-1);
124 }
125
scm_number2long(SOBJ x)126 long scm_number2long(SOBJ x)
127 {
128 if (!SCM_NUMBERP(x)) SCM_ERR("scm_number2long: bad number", x);
129 switch(SCM_OBJTYPE(x)) {
130 case SOBJ_T_INUM: return( SCM_INUM(x) );
131 case SOBJ_T_FNUM: return( (long)SCM_FNUM(x) );
132 case SOBJ_T_BNUM: return( mpz_get_ui(SCM_BNUM(x)) );
133 }
134 return(-1);
135 }
136
137 /*-- convert to same type : return the type */
scm_convert(SOBJ * px,SOBJ * py)138 static int scm_convert(SOBJ *px, SOBJ *py)
139 {
140 SOBJ x = *px;
141 SOBJ y = *py;
142
143 if (SCM_OBJTYPE(x) == SCM_OBJTYPE(y)) return(SCM_OBJTYPE(x));
144
145 if (SCM_OBJTYPE(x) == SOBJ_T_FNUM && ISINT(SCM_FNUM(x))) {
146 *px = x = scm_flt2num(SCM_FNUM(x));
147 }
148 if (SCM_OBJTYPE(y) == SOBJ_T_FNUM && ISINT(SCM_FNUM(y))) {
149 *py = y = scm_flt2num(SCM_FNUM(y));
150 }
151
152 if (SCM_OBJTYPE(x) == SCM_OBJTYPE(y)) return(SCM_OBJTYPE(x));
153
154 switch(SCM_OBJTYPE(x)) {
155 case SOBJ_T_INUM:
156 switch(SCM_OBJTYPE(y)) {
157 case SOBJ_T_BNUM: *px = scm_int2bnum(SCM_INUM(x)); return(SOBJ_T_BNUM);
158 case SOBJ_T_FNUM: *px = int2fnum(SCM_INUM(x)); return(SOBJ_T_FNUM);
159 }
160 case SOBJ_T_BNUM:
161 switch(SCM_OBJTYPE(y)) {
162 case SOBJ_T_INUM: *py = scm_int2bnum(SCM_INUM(y)); return(SOBJ_T_BNUM);
163 case SOBJ_T_FNUM: *px = big2fnum(SCM_BNUM(x)); return(SOBJ_T_FNUM);
164 }
165 case SOBJ_T_FNUM:
166 switch(SCM_OBJTYPE(y)) {
167 case SOBJ_T_INUM: *py = int2fnum(SCM_INUM(y)); return(SOBJ_T_FNUM);
168 case SOBJ_T_BNUM: *py = big2fnum(SCM_BNUM(y)); return(SOBJ_T_FNUM);
169 }
170 }
171 return(SOBJ_T_UNDEFINED);
172 }
173
174 /*-- choose the most compact representation of the number */
scm_compact_number(SOBJ n)175 SOBJ scm_compact_number(SOBJ n)
176 {
177 if (SCM_INUMP(n)) return(n);
178 if (SCM_FNUMP(n) && ISINT(SCM_FNUM(n)))
179 n = scm_flt2num(SCM_FNUM(n));
180
181 if (SCM_INUMP(n)) return(n);
182 switch(n->type) {
183 case SOBJ_T_BNUM:
184 if (mpz_cmp_si(SCM_BNUM(n), SOBJ_INUM_MIN) >= 0 &&
185 mpz_cmp_si(SCM_BNUM(n), SOBJ_INUM_MAX) <= 0) {
186 n = SCM_MKINUM(mpz_get_si(SCM_BNUM(n)));
187 }
188 break;
189 case SOBJ_T_FNUM:
190 if (ISINT(SCM_FNUM(n))) n = scm_flt2num(SCM_FNUM(n));
191 break;
192 }
193 return(n);
194 }
195
196 /****************************************************************
197 * Arithmetic
198 ****************************************************************/
199
200 /*-- add 2 number */
scm_add2(SOBJ x,SOBJ y)201 SOBJ scm_add2(SOBJ x, SOBJ y)
202 {
203 SOBJ r;
204 switch(scm_convert(&x, &y)) {
205 case SOBJ_T_INUM: return(scm_int2num(SCM_INUM(x) + SCM_INUM(y)));
206 case SOBJ_T_FNUM: r = scm_mkfnum(SCM_FNUM(x) + SCM_FNUM(y)); break;
207 case SOBJ_T_BNUM:
208 r = scm_newbnum();
209 mpz_add(SCM_BNUM(r), SCM_BNUM(x), SCM_BNUM(y));
210 break;
211
212 default:
213 return(scm_undefined);
214 }
215 return(scm_compact_number(r));
216 }
217
218 /*-- mul 2 number */
scm_mul2(SOBJ x,SOBJ y)219 SOBJ scm_mul2(SOBJ x, SOBJ y)
220 {
221 SOBJ r;
222 switch(scm_convert(&x, &y)) {
223 case SOBJ_T_FNUM: r = scm_mkfnum(SCM_FNUM(x) * SCM_FNUM(y)); break;
224 case SOBJ_T_INUM:
225 {
226 long n;
227
228 if (SCM_INUM(x) == 0 || SCM_INUM(y) == 0) return(SCM_MKINUM(0));
229 n = SCM_INUM(x) * SCM_INUM(y);
230 if (SCM_INUM(x) == (n / SCM_INUM(y))) return(scm_int2num(n));
231
232 r = scm_mkfnum((double)SCM_INUM(x) * (double)SCM_INUM(y));
233 break;
234 }
235 case SOBJ_T_BNUM:
236 r = scm_newbnum();
237 mpz_mul(SCM_BNUM(r), SCM_BNUM(x), SCM_BNUM(y));
238 break;
239
240 default:
241 return(scm_undefined);
242 }
243 return(scm_compact_number(r));
244 }
245
246 /*-- sub 2 number */
scm_sub2(SOBJ x,SOBJ y)247 SOBJ scm_sub2(SOBJ x, SOBJ y)
248 {
249 SOBJ r;
250 switch(scm_convert(&x, &y)) {
251 case SOBJ_T_INUM: return(scm_int2num(SCM_INUM(x) - SCM_INUM(y)));
252 case SOBJ_T_FNUM: r = scm_mkfnum(SCM_FNUM(x) - SCM_FNUM(y)); break;
253 case SOBJ_T_BNUM:
254 r = scm_newbnum();
255 mpz_sub(SCM_BNUM(r), SCM_BNUM(x), SCM_BNUM(y));
256 break;
257
258 default:
259 return(scm_undefined);
260 }
261 return(scm_compact_number(r));
262 }
263
264 /*-- div 2 numbers */
scm_div2(SOBJ x,SOBJ y)265 SOBJ scm_div2(SOBJ x, SOBJ y)
266 {
267 SOBJ r;
268 switch(scm_convert(&x, &y)) {
269 case SOBJ_T_FNUM: r = scm_mkfnum(SCM_FNUM(x) / SCM_FNUM(y)); break;
270 case SOBJ_T_INUM:
271 {
272 double tmp = (double)SCM_INUM(x) / (double)SCM_INUM(y);
273 if (ISINT(tmp)) {
274 if (SCM_INUM_RANGE(tmp)) return(SCM_MKINUM((long)tmp));
275 r = scm_newbnum();
276 mpz_set_d(SCM_BNUM(r), tmp);
277 } else {
278 r = scm_mkfnum(tmp);
279 }
280 break;
281 }
282 case SOBJ_T_BNUM:
283 r = scm_newbnum();
284 mpz_div(SCM_BNUM(r), SCM_BNUM(x), SCM_BNUM(y));
285 break;
286
287 default:
288 return(scm_undefined);
289 }
290 return(scm_compact_number(r));
291 }
292
293 /*-- compare 2 numbers */
scm_cmpnum(SOBJ x,SOBJ y)294 int scm_cmpnum(SOBJ x, SOBJ y)
295 {
296 int t;
297
298 if (x == y) return(0); /* quick out when eq */
299
300 t = scm_convert(&x, &y);
301 if (x == y) return(0);
302 switch(t) {
303 case SOBJ_T_INUM:
304 return(((long)x > (long)y) ? 1 : -1 );
305 case SOBJ_T_FNUM:
306 if (SCM_FNUM(x) == SCM_FNUM(y)) return(0);
307 return((SCM_FNUM(x) > SCM_FNUM(y)) ? 1 : -1);
308
309 case SOBJ_T_BNUM:
310 return(mpz_cmp(SCM_BNUM(x), SCM_BNUM(y)));
311 }
312 return(SOBJ_INUM_MAX);
313 }
314
315 /*-- cmpops */
scm_lt2(SOBJ x,SOBJ y)316 SOBJ scm_lt2(SOBJ x, SOBJ y) { return(SCM_MKBOOL(scm_cmpnum(x,y) < 0)); }
scm_le2(SOBJ x,SOBJ y)317 SOBJ scm_le2(SOBJ x, SOBJ y) { return(SCM_MKBOOL(scm_cmpnum(x,y) <= 0)); }
scm_ge2(SOBJ x,SOBJ y)318 SOBJ scm_ge2(SOBJ x, SOBJ y) { return(SCM_MKBOOL(scm_cmpnum(x,y) >= 0)); }
scm_gt2(SOBJ x,SOBJ y)319 SOBJ scm_gt2(SOBJ x, SOBJ y) { return(SCM_MKBOOL(scm_cmpnum(x,y) > 0)); }
scm_eq2(SOBJ x,SOBJ y)320 SOBJ scm_eq2(SOBJ x, SOBJ y) { return(SCM_MKBOOL(scm_cmpnum(x,y) == 0)); }
321
scm_zerop(SOBJ x)322 SOBJ scm_zerop(SOBJ x)
323 {
324 return(SCM_MKBOOL(scm_cmpnum(x, SCM_MKINUM(0)) == 0));
325 }
326
scm_positivep(SOBJ x)327 SOBJ scm_positivep(SOBJ x)
328 {
329 return(SCM_MKBOOL(scm_cmpnum(x, SCM_MKINUM(0)) > 0));
330 }
331
scm_negativep(SOBJ x)332 SOBJ scm_negativep(SOBJ x)
333 {
334 return(SCM_MKBOOL(scm_cmpnum(x, SCM_MKINUM(0)) < 0));
335 }
336
scm_oddp(SOBJ x)337 SOBJ scm_oddp(SOBJ x)
338 {
339 if (!SCM_NUMBERP(x)) SCM_ERR("odd?: bad number", x);
340 switch(SCM_OBJTYPE(x)) {
341 case SOBJ_T_INUM: return( SCM_MKBOOL( (SCM_INUM(x) & 1) == 1 ));
342 case SOBJ_T_FNUM: return( SCM_MKBOOL( fmod(SCM_FNUM(x),2.0) == 1.0 ));
343 case SOBJ_T_BNUM:
344 {
345 MP_INT q, r;
346 SOBJ res;
347 mpz_init(&q); mpz_init(&r);
348 mpz_divmod_ui(&q, &r, SCM_BNUM(x), 2);
349 res = SCM_MKBOOL( mpz_cmp_ui(&r, 0) != 0);
350 mpz_clear(&q); mpz_clear(&r);
351 return(res);
352 }
353 }
354 SCM_ERR("odd?: bad number", x);
355 return(NULL);
356 }
357
scm_evenp(SOBJ x)358 SOBJ scm_evenp(SOBJ x)
359 {
360 if (!SCM_NUMBERP(x)) SCM_ERR("even?: bad number", x);
361 switch(SCM_OBJTYPE(x)) {
362 case SOBJ_T_INUM: return( SCM_MKBOOL( (SCM_INUM(x) & 1) == 0 ));
363 case SOBJ_T_FNUM: return( SCM_MKBOOL( fmod(SCM_FNUM(x),2.0) == 0.0 ));
364 case SOBJ_T_BNUM:
365 {
366 MP_INT q, r;
367 SOBJ res;
368 mpz_init(&q); mpz_init(&r);
369 mpz_divmod_ui(&q, &r, SCM_BNUM(x), 2);
370 res = SCM_MKBOOL( mpz_cmp_ui(&r, 0) == 0);
371 mpz_clear(&q); mpz_clear(&r);
372 return(res);
373 }
374 }
375 SCM_ERR("even?: bad number", x);
376 return(NULL);
377 }
378
scm_abs(SOBJ x)379 SOBJ scm_abs(SOBJ x)
380 {
381 if (!SCM_NUMBERP(x)) SCM_ERR("abs: bad number", x);
382 switch(SCM_OBJTYPE(x)) {
383 case SOBJ_T_INUM:
384 return( SCM_INUM(x) < 0 ? SCM_MKINUM(-(SCM_INUM(x))) : x);
385 case SOBJ_T_FNUM:
386 return( SCM_FNUM(x) < 0.0? scm_mkfnum( -(SCM_FNUM(x)) ) : x);
387 case SOBJ_T_BNUM:
388 if (mpz_cmp_ui(SCM_BNUM(x), 0) < 0) {
389 SOBJ new = scm_newbnum();
390 mpz_neg(SCM_BNUM(new), SCM_BNUM(x));
391 return(new);
392 }
393 return(x);
394 }
395 return(NULL);
396 }
397
scm_do_int_div(SOBJ x,SOBJ y,int quotient)398 SOBJ scm_do_int_div(SOBJ x, SOBJ y, int quotient)
399 {
400 MP_INT q, r;
401 SOBJ res;
402 int exact = TRUE;
403
404 if (SCM_FNUMP(x)) {
405 if (!ISINT(SCM_FNUM(x))) SCM_ERR("bad integer", x);
406 x = scm_flt2num(SCM_FNUM(x)); exact = FALSE;
407 }
408 if (SCM_FNUMP(y)) {
409 if (!ISINT(SCM_FNUM(y))) SCM_ERR("bad integer", y);
410 y = scm_flt2num(SCM_FNUM(y)); exact = FALSE;
411 }
412
413 if (SCM_INUMP(x) && SCM_INUMP(y)) {
414 long r = quotient ?
415 SCM_INUM(x) / SCM_INUM(y) :
416 SCM_INUM(x) % SCM_INUM(y) ;
417
418 return( exact ? SCM_MKINUM(r) : scm_mkfnum(r));
419 }
420 /* at least one of the 2 number is big */
421 if (SCM_INUMP(x)) x = scm_int2bnum(SCM_INUM(x));
422 if (SCM_INUMP(y)) y = scm_int2bnum(SCM_INUM(y));
423
424 res = scm_newbnum();
425 mpz_init(&q); mpz_init(&r);
426 /* mpz_divmod(&q, &r, SCM_BNUM(x), SCM_BNUM(y)); */
427 mpz_tdiv_qr(&q, &r, SCM_BNUM(x), SCM_BNUM(y));
428 mpz_set(SCM_BNUM(res), quotient ? &q : &r);
429 mpz_clear(&q); mpz_clear(&r);
430 return( exact ? res : big2fnum(SCM_BNUM(res)));
431 }
432
scm_quotient(SOBJ x,SOBJ y)433 SOBJ scm_quotient(SOBJ x, SOBJ y)
434 {
435 if (!SCM_NUMBERP(x) && !SCM_NUMBERP(y))
436 SCM_ERR("quotient: bad arguments", scm_cons(x, y));
437
438 if (scm_zerop(y) == scm_true) SCM_ERR("quotient: division by 0", NULL);
439
440 return(scm_do_int_div(x, y, TRUE));
441 }
442
scm_remainder(SOBJ x,SOBJ y)443 SOBJ scm_remainder(SOBJ x, SOBJ y)
444 {
445 if (!SCM_NUMBERP(x) && !SCM_NUMBERP(y))
446 SCM_ERR("remainder: bad arguments", scm_cons(x, y));
447
448 if (scm_zerop(y) == scm_true) SCM_ERR("remainder: division by 0", NULL);
449 return(scm_do_int_div(x, y, FALSE));
450 }
451
scm_modulo(SOBJ x,SOBJ y)452 SOBJ scm_modulo(SOBJ x, SOBJ y)
453 {
454 SOBJ z;
455
456 if (!SCM_NUMBERP(x) && !SCM_NUMBERP(y))
457 SCM_ERR("modulo: bad arguments", scm_cons(x, y));
458
459 if (scm_zerop(y) == scm_true) SCM_ERR("modulo: division by 0", NULL);
460
461 z = scm_do_int_div(x, y, FALSE);
462 if (scm_negativep(x) != scm_negativep(y) && scm_zerop(z) != scm_true) {
463 z = scm_add2(z, y);
464 }
465 return(z);
466 }
467
scm_gcd2(SOBJ x,SOBJ y)468 static SOBJ scm_gcd2(SOBJ x, SOBJ y)
469 {
470 return( scm_zerop(y) == scm_true ? x : scm_gcd2(y, scm_modulo(x,y)));
471 }
472
scm_gcd(int nargs,SOBJ * obj)473 SOBJ scm_gcd(int nargs, SOBJ *obj)
474 {
475 int i=0;
476 SOBJ res = SCM_MKINUM(0);
477 while(i < nargs) {
478 if (!SCM_NUMBERP(obj[i])) SCM_ERR("gcd: bad number", obj[i]);
479 res = (i == 0) ?
480 scm_abs(obj[i++]) :
481 scm_abs(scm_gcd2(res, obj[i++]));
482 }
483 return(res);
484 }
485
scm_lcm(int nargs,SOBJ * obj)486 SOBJ scm_lcm(int nargs, SOBJ *obj)
487 {
488 int i=0;
489 SOBJ res = SCM_MKINUM(1);
490 while(i < nargs) {
491 if (!SCM_NUMBERP(obj[i])) SCM_ERR("lcm: bad number", obj[i]);
492 if (i == 0) {
493 res = scm_abs(obj[i]);
494 } else {
495 if (scm_zerop(obj[i]) == scm_true) return(obj[i]);
496 res = scm_abs(scm_div2(scm_mul2(res, obj[i]), scm_gcd2(res, obj[i])));
497 }
498 i++;
499 }
500 return(res);
501 }
502
503 /*-- float specific */
504
scm_floor(SOBJ x)505 SOBJ scm_floor(SOBJ x)
506 {
507 if (!SCM_NUMBERP(x)) SCM_ERR("floor: bad number", x);
508 if (SCM_FNUMP(x)) return(scm_mkfnum(floor(SCM_FNUM(x))));
509 return(x);
510 }
511
scm_ceil(SOBJ x)512 SOBJ scm_ceil(SOBJ x)
513 {
514 if (!SCM_NUMBERP(x)) SCM_ERR("ceil: bad number", x);
515 if (SCM_FNUMP(x)) return(scm_mkfnum(ceil(SCM_FNUM(x))));
516 return(x);
517 }
518
scm_truncate(SOBJ x)519 SOBJ scm_truncate(SOBJ x)
520 {
521 if (!SCM_NUMBERP(x)) SCM_ERR("truncate: bad number", x);
522 if (SCM_FNUMP(x)) {
523 double d = SCM_FNUM(x);
524 return(scm_mkfnum( d < 0.0 ? ceil(d) : floor(d) ));
525 }
526 return(x);
527 }
528
scm_round(SOBJ x)529 SOBJ scm_round(SOBJ x)
530 {
531 if (!SCM_NUMBERP(x)) SCM_ERR("truncate: bad number", x);
532 if (SCM_FNUMP(x)) {
533 double v = SCM_FNUM(x);
534 double v_plus_0_5 = v + 0.5;
535 double res = floor(v_plus_0_5);
536
537 return(scm_mkfnum(
538 (v_plus_0_5==res && v_plus_0_5/2 != floor(v_plus_0_5/2)) ? res-1: res));
539 }
540 return(x);
541 }
542
exact2inexact(SOBJ x)543 static double exact2inexact(SOBJ x)
544 {
545 if (!SCM_NUMBERP(x)) SCM_ERR("exact->inexact: bad number", x);
546 switch(SCM_OBJTYPE(x)) {
547 case SOBJ_T_INUM: return( (double) SCM_INUM(x));
548 case SOBJ_T_FNUM: return( SCM_FNUM(x));
549 case SOBJ_T_BNUM: return( mpz_get_d(SCM_BNUM(x)));
550 }
551 return(0);
552 }
553
554 #define TRANS(fn) \
555 SOBJ scm_##fn(SOBJ x) { \
556 if (SCM_NUMBERP(x)) return(scm_mkfnum(fn(exact2inexact(x)))); \
557 SCM_ERR("trans: bad number", x); return(NULL);\
558 }
559
560 TRANS(exp);
561 TRANS(log);
562 TRANS(log10);
563 TRANS(sin);
564 TRANS(cos);
565 TRANS(tan);
566 TRANS(asin);
567 TRANS(acos);
568
scm_atan(SOBJ x,SOBJ y)569 SOBJ scm_atan(SOBJ x, SOBJ y)
570 {
571 double dx, dy;
572
573 if (!SCM_NUMBERP(x)) SCM_ERR("atan: bad number", x);
574 if (!SCM_NUMBERP(y)) return(scm_mkfnum(atan(exact2inexact(y))));
575
576 dx = exact2inexact(x); dy = exact2inexact(y);
577 return(scm_mkfnum( (dx == 0 && dy == 0) ? 0 : atan2(dx, dy)));
578 }
579
scm_sqrt(SOBJ x)580 SOBJ scm_sqrt(SOBJ x)
581 {
582 double d;
583
584 if (!SCM_NUMBERP(x)) SCM_ERR("sqrt: bad number", x);
585 switch(SCM_OBJTYPE(x)) {
586 case SOBJ_T_INUM:
587 if (SCM_INUM(x) < 0) goto errneg;
588 d = sqrt(SCM_INUM(x));
589 return( ISINT(d) ? SCM_MKINUM( (long)d) : scm_mkfnum(d) );
590
591 case SOBJ_T_FNUM:
592 if (SCM_FNUM(x) < 0.0) goto errneg;
593 return(scm_mkfnum(sqrt(SCM_FNUM(x))));
594
595 case SOBJ_T_BNUM: {
596 MP_INT root, rem;
597 SOBJ res = NULL;
598
599 if (mpz_cmp_si(SCM_BNUM(x), 0) < 0) goto errneg;
600 mpz_init(&root); mpz_init(&rem);
601 mpz_sqrtrem(&root, &rem, SCM_BNUM(x));
602 if (mpz_cmp_si(&rem, 0) == 0) {
603 res = scm_newbnum();
604 mpz_set(SCM_BNUM(res),&root);
605 } else {
606 res = big2fnum(&root);
607 }
608 mpz_clear(&root); mpz_clear(&rem);
609 return(res);
610 }
611 }
612 errneg:
613 SCM_ERR("sqrt: negative number", x);
614 return(NULL);
615 }
616
scm_exact_expt(SOBJ x,SOBJ y)617 SOBJ scm_exact_expt(SOBJ x, SOBJ y)
618 {
619 MP_INT mod;
620 SOBJ new = scm_newbnum();
621
622 mpz_init_set_ui(&mod, 1);
623 mpz_powm(SCM_BNUM(new), SCM_BNUM(x), SCM_BNUM(y), &mod);
624 mpz_clear(&mod);
625 return(new);
626 }
627
scm_expt(SOBJ x,SOBJ y)628 SOBJ scm_expt(SOBJ x, SOBJ y)
629 {
630 if (!SCM_NUMBERP(x)) SCM_ERR("expt: bad number", x);
631 if (!SCM_NUMBERP(y)) SCM_ERR("expt: bad number", y);
632
633 switch(scm_convert(&x, &y)) {
634 case SOBJ_T_INUM: {
635 SOBJ new = scm_newbnum();
636 mpz_ui_pow_ui(SCM_BNUM(new), SCM_INUM(x), SCM_INUM(y));
637 return(scm_compact_number(new));
638 }
639 case SOBJ_T_FNUM:
640 return(scm_mkfnum(pow(SCM_FNUM(x), SCM_FNUM(y))));
641
642 case SOBJ_T_BNUM: {
643 SOBJ new = scm_newbnum();
644 MP_INT mod;
645 mpz_init_set_ui(&mod, 1);
646 mpz_powm(SCM_BNUM(new), SCM_BNUM(x), SCM_BNUM(y), &mod);
647 mpz_clear(&mod);
648 return(new);
649 }
650 }
651 return(NULL);
652 }
653
scm_exact_to_inexact(SOBJ x)654 SOBJ scm_exact_to_inexact(SOBJ x)
655 {
656 return(scm_mkfnum(exact2inexact(x)));
657 }
658
scm_inexact_to_exact(SOBJ x)659 SOBJ scm_inexact_to_exact(SOBJ x)
660 {
661 if (!SCM_NUMBERP(x)) SCM_ERR("inexact->exact: bad number", x);
662 if (SCM_FNUMP(x)) return(scm_flt2num(SCM_FNUM(x)));
663 return(x);
664 }
665
666 /*-- return mallocated string */
scm_number2cstr(SOBJ n,long base)667 char *scm_number2cstr(SOBJ n, long base)
668 {
669 char *s;
670 char buf[128];
671
672 switch(SCM_OBJTYPE(n)) {
673 case SOBJ_T_FNUM:
674 if (base != 10) SCM_ERR("base must be 10 for this number", n);
675 sprintf(buf, "%.15g", SCM_FNUM(n));
676 return(scm_must_strdup(buf));
677
678 case SOBJ_T_INUM:
679 {
680 char *p;
681 long value;
682 int digit;
683 if (base < 1 || base > 36) SCM_ERR("bad base (valid range 1..36)", n);
684
685 p = buf + sizeof(buf);
686 *(--p) = 0;
687 value = SCM_INUM(n);
688 if (value < 0) value = -value;
689 do {
690 digit = value % base;
691 *(--p) = digit + ((digit < 10) ? '0' : 'a'-10);
692 value = value / base;
693 } while(value != 0);
694 if (SCM_INUM(n) < 0) {
695 *(--p) = '-';
696 }
697 return(scm_must_strdup(p));
698 }
699 case SOBJ_T_BNUM:
700 if (base < 1 || base > 36) SCM_ERR("bad base (valid range 1..36)", n);
701 s = mpz_get_str(NULL, base, SCM_BNUM(n));
702 return(s);
703 }
704 return NULL; /* never reached */
705 }
706
scm_cstr2number(char * s,int base)707 SOBJ scm_cstr2number(char *s, int base)
708 {
709 int exact, isint;
710 char *p;
711
712 exact = 0;
713 if (*s == '#') { /* got specs */
714 s++;
715 while(*s) {
716 switch(tolower(*s)) {
717 case 'e': exact = 1; s++; continue;
718 case 'i': exact = -1; s++; continue;
719 case 'b': base = 2; s++; continue;
720 case 'o': base = 8; s++; continue;
721 case 'd': base = 10; s++; continue;
722 case 'x': base = 16; s++; continue;
723 }
724 break;
725 }
726 }
727 if (*s == 0) return(scm_false);
728
729 if (exact == -1 && base != 10) SCM_ERR("inexact number, only base 10", NULL);
730
731 if (*s == '+') s++;
732
733 isint = 1; p = s;
734 while(*p) {
735 if ( *p == '-' || *p == '+') {
736 p++;
737 continue;
738 }
739 if ( (base==2 && (*p == '0' || *p == '1')) ||
740 (base==8 && (strchr("01234567", *p) != NULL)) ||
741 (base==10 && isdigit( *p)) ||
742 (base==16 && isxdigit(*p))) {
743 p++;
744 continue;
745 }
746 if (strchr("#.Ee", *p)) { /* float */
747 if (exact == 1) SCM_ERR("illegal . in exact number", NULL);
748 isint = 0;
749 if (*p == '#') *p = '0';
750 p++;
751 continue;
752 }
753 if (strchr("sSfFdDlL", *p)) { /* precision marker ignored */
754 *p = 0; break;
755 }
756 scm_puts("strange char '"); scm_putc(*p); scm_puts("' in number\n");
757 return(scm_false);
758 }
759 if (isint && exact != -1) {
760 MP_INT n;
761 SOBJ z;
762 if (mpz_init_set_str(&n, s, base) < 0) {
763 mpz_clear(&n);
764 return(scm_false);
765 }
766 if (mpz_cmp_si(&n, SOBJ_INUM_MIN) >= 0 &&
767 mpz_cmp_si(&n, SOBJ_INUM_MAX) <= 0) {
768 long num = mpz_get_si(&n);
769 mpz_clear(&n);
770 return(SCM_MKINUM(num));
771 }
772 z = scm_newbnum();
773 mpz_set(SCM_BNUM(z), &n);
774 return z;
775 } else {
776 return(scm_mkfnum(atof(s)));
777 }
778 }
779
scm_number_to_string(SOBJ x,SOBJ y)780 SOBJ scm_number_to_string(SOBJ x, SOBJ y)
781 {
782 int base;
783 char *p;
784
785 if (!SCM_NUMBERP(x)) SCM_ERR("number->string: bad number", x);
786 if (y == NULL) {
787 base = 10;
788 } else {
789 if (!SCM_INUMP(y)) SCM_ERR("number->string: bad base", y);
790 base = SCM_INUM(y);
791 }
792 p = scm_number2cstr(x, base);
793 x = scm_newcell(SOBJ_T_STRING);
794 SCM_STR_LEN(x) = strlen(p);
795 SCM_STR_VALUE(x) = p;
796 return(x);
797 }
798
scm_string_to_number(SOBJ x,SOBJ y)799 SOBJ scm_string_to_number(SOBJ x, SOBJ y)
800 {
801 int base;
802
803 if (!SCM_STRINGP(x)) SCM_ERR("string->number: bad string", x);
804 if (y == NULL) {
805 base = 10;
806 } else {
807 if (!SCM_INUMP(y)) SCM_ERR("string->number: bad base", y);
808 base = SCM_INUM(y);
809 }
810 return(scm_cstr2number(SCM_STR_VALUE(x), base));
811 }
812
813 /****************************************************************
814 * Bit operations
815 ****************************************************************/
816 /*E* (ashift NUMBER POS) => NUMBER */
817 /*D* Return NUMBER shifted by POS bits. If POS is positive, shift to
818 the left, else shift to the right. Equivalent to NUMBER * expt(2,
819 POS) */
scm_ashift(SOBJ val,SOBJ n)820 SOBJ scm_ashift(SOBJ val, SOBJ n)
821 {
822 long m, k;
823 m = scm_number2long(val);
824 k = scm_number2long(n);
825 if (k > 0) m = m << k;
826 else m = m >> (-k);
827 return(scm_int2num(m));
828 }
829
830 /*E* (bit-and N ...) => NUMBER */
831 /*D* Returns the logical and of it's argument. */
832
scm_bit_and(int narg,SOBJ arg[])833 SOBJ scm_bit_and(int narg, SOBJ arg[])
834 {
835 int i;
836 long r;
837
838 if (narg < 1) return(scm_undefined);
839 r = scm_number2long(arg[0]);
840 for (i = 1; i < narg; i++) r &= scm_number2long(arg[i]);
841 return(scm_int2num(r));
842 }
843
844 /*E* (bit-or N ...) => NUMBER */
845 /*D* Returns the logical or of it's argument. */
846
scm_bit_or(int narg,SOBJ arg[])847 SOBJ scm_bit_or(int narg, SOBJ arg[])
848 {
849 int i;
850 long r;
851
852 if (narg < 1) return(scm_undefined);
853 r = scm_number2long(arg[0]);
854 for (i = 1; i < narg; i++) r |= scm_number2long(arg[i]);
855 return(scm_int2num(r));
856 }
857
858 /*E* (bit-xor N ...) => NUMBER */
859 /*D* Returns the logical exclusive or of it's argument. */
860
scm_bit_xor(int narg,SOBJ arg[])861 SOBJ scm_bit_xor(int narg, SOBJ arg[])
862 {
863 int i;
864 long r;
865
866 if (narg < 1) return(scm_undefined);
867 r = scm_number2long(arg[0]);
868 for (i = 1; i < narg; i++) r ^= scm_number2long(arg[i]);
869 return(scm_int2num(r));
870 }
871
872 /*E* (bit-not N) => NUMBER */
873 /*D* Returns the logical not of N. */
874
scm_bit_not(SOBJ n)875 SOBJ scm_bit_not(SOBJ n)
876 {
877 return(scm_int2num(~(scm_number2long(n))));
878 }
879
880 /****************************************************************
881 * scheme interface
882 ****************************************************************/
883
scm_init_number()884 void scm_init_number()
885 {
886 scm_add_cprim("ashift", scm_ashift, 2);
887 scm_add_cprim("bit-and", scm_bit_and, -1);
888 scm_add_cprim("bit-or", scm_bit_or, -1);
889 scm_add_cprim("bit-xor", scm_bit_xor, -1);
890 scm_add_cprim("bit-not", scm_bit_not, 1);
891 }
892