1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <float.h>
4 #include <math.h>
5 #include "eisl.h"
6
7
8 // data type transfer
9
10 int
exact_to_inexact(int x)11 exact_to_inexact(int x)
12 {
13 int res,
14 tag;
15 double val;
16
17 tag = GET_TAG(x);
18 switch (tag) {
19 case INTN:
20 res = freshcell();
21 SET_TAG(res, FLTN);
22 val = (double) GET_INT(x);
23 SET_FLT(res, val);
24 return (res);
25 case LONGN:
26 res = freshcell();
27 SET_TAG(res, FLTN);
28 val = (double) GET_LONG(x);
29 SET_FLT(res, val);
30 return (res);
31 case BIGX:
32 return (bigx_big_to_flt(x));
33 case FLTN:
34 return (x);
35 default:
36 return (x);
37 }
38 error(OUT_OF_RANGE, "exact_to_inexact", x);
39 return (UNDEF);
40 }
41
42 int
numeqp(int x,int y)43 numeqp(int x, int y)
44 {
45 int x1,
46 y1;
47
48 // if integer type and address are same, then true.
49 if (integerp(x) && integerp(y)) {
50 if (GET_INT(x) == GET_INT(y))
51 return (1);
52 else
53 return (0);
54 }
55 // if long integer data type and the values are same, then true.
56 else if (longnump(x) && longnump(y)) {
57 if (GET_LONG(x) == GET_LONG(y))
58 return (1);
59 else
60 return (0);
61 }
62
63 else if (bignump(x) && bignump(y)) {
64 if (bigx_eqp(x, y))
65 return (1);
66 else
67 return (0);
68 }
69 // if floating point number, consider EPSILON
70 else if (floatp(x) && floatp(y)) {
71 if (GET_FLT(x) - DBL_EPSILON <= GET_FLT(y) &&
72 GET_FLT(x) + DBL_EPSILON >= GET_FLT(y))
73 return (1);
74 else if (GET_FLT(x) >= GET_FLT(y) - DBL_EPSILON &&
75 GET_FLT(x) <= GET_FLT(y) + DBL_EPSILON)
76 return (1);
77 else
78 return (0);
79 } else {
80 x1 = exact_to_inexact(x);
81 y1 = exact_to_inexact(y);
82
83 if (GET_FLT(x1) - DBL_EPSILON <= GET_FLT(y1) &&
84 GET_FLT(x1) + DBL_EPSILON >= GET_FLT(y1))
85 return (1);
86 else if (GET_FLT(x1) >= GET_FLT(y1) - DBL_EPSILON &&
87 GET_FLT(x1) <= GET_FLT(y1) - DBL_EPSILON)
88 return (1);
89 else
90 return (0);
91 }
92 }
93
94 int
smallerp(int x1,int x2)95 smallerp(int x1, int x2)
96 {
97 if (IS_INTEGER(x1) && IS_INTEGER(x2)) {
98 if (GET_INT(x1) < GET_INT(x2))
99 return (1);
100 else
101 return (0);
102 } else if (integerp(x1) && bignump(x2))
103 return (1);
104 else if (longnump(x1) && bignump(x2))
105 return (1);
106 else if (longnump(x1) && longnump(x2)) {
107 if (GET_LONG(x1) < GET_LONG(x2))
108 return (1);
109 else
110 return (0);
111 } else if (bignump(x1) && bignump(x2)) {
112 if (bigx_smallerp(x1, x2))
113 return (1);
114 else
115 return (0);
116 } else {
117 if (GET_FLT(exact_to_inexact(x1)) < GET_FLT(exact_to_inexact(x2)))
118 return (1);
119 else
120 return (0);
121 }
122 }
123
124 int
eqsmallerp(int x1,int x2)125 eqsmallerp(int x1, int x2)
126 {
127 if (integerp(x1) && integerp(x2)) {
128 if (GET_INT(x1) <= GET_INT(x2))
129 return (1);
130 else
131 return (0);
132 } else if (integerp(x1) && bignump(x2))
133 return (1);
134 else if (longnump(x1) && bignump(x2))
135 return (1);
136 else if (longnump(x1) && longnump(x2)) {
137 if (GET_LONG(x1) <= GET_LONG(x2))
138 return (1);
139 else
140 return (0);
141 } else if (bignump(x1) && bignump(x2)) {
142 if (bigx_smallerp(x1, x2) || bigx_eqp(x1, x2))
143 return (1);
144 else
145 return (0);
146 } else {
147 if (GET_FLT(exact_to_inexact(x1)) <= GET_FLT(exact_to_inexact(x2)))
148 return (1);
149 else
150 return (0);
151 }
152 }
153
154 int
greaterp(int x1,int x2)155 greaterp(int x1, int x2)
156 {
157 if (smallerp(x2, x1))
158 return (1);
159 else
160 return (0);
161 }
162
163 int
eqgreaterp(int x1,int x2)164 eqgreaterp(int x1, int x2)
165 {
166 if (eqsmallerp(x2, x1))
167 return (1);
168 else
169 return (0);
170 }
171
172 int
positivep(int x)173 positivep(int x)
174 {
175
176 if (integerp(x) && GET_INT(x) > 0)
177 return (1);
178 else if (longnump(x) && GET_LONG(x) > 0)
179 return (1);
180 else if (bignump(x) && bigx_positivep(x))
181 return (1);
182 else if (floatp(x) && GET_FLT(x) > 0)
183 return (1);
184 else
185 return (0);
186 }
187
188 int
negativep(int x)189 negativep(int x)
190 {
191
192 if (integerp(x) && GET_INT(x) < 0)
193 return (1);
194 else if (longnump(x) && GET_LONG(x) < 0)
195 return (1);
196 else if (bignump(x) && bigx_negativep(x))
197 return (1);
198 else if (floatp(x) && GET_FLT(x) < 0)
199 return (1);
200 else
201 return (0);
202 }
203
204 int
zerop(int x)205 zerop(int x)
206 {
207 if (integerp(x) && GET_INT(x) == 0)
208 return (1);
209 else if (floatp(x) && GET_FLT(x) == 0.0)
210 return (1);
211 else
212 return (0);
213 }
214
215 int
positive_zerop(int x)216 positive_zerop(int x)
217 {
218 if (zerop(x) && floatp(x) && !signbit(GET_FLT(x)))
219 return (1);
220 else if (zerop(x) && floatp(x))
221 return (1);
222 else
223 return (0);
224 }
225
226 int
negative_zerop(int x)227 negative_zerop(int x)
228 {
229 if (zerop(x) && floatp(x) && signbit(GET_FLT(x)))
230 return (1);
231 else
232 return (0);
233 }
234
235
236 // -----------------------------------
237 // basic operation
238
239 int
plus(int arg1,int arg2)240 plus(int arg1, int arg2)
241 {
242 int tag1,
243 tag2;
244 double x1,
245 x2;
246
247 tag1 = GET_TAG(arg1);
248 tag2 = GET_TAG(arg2);
249 switch (tag1) {
250 case INTN:
251 switch (tag2) {
252 case INTN:{
253 long long int l;
254
255 l = (long long int) GET_INT(arg1) +
256 (long long int) GET_INT(arg2);
257 if (SMALL_INT_MIN < l && l < SMALL_INT_MAX)
258 return (makeint((int) l));
259 else
260 return (makelong(l));
261 }
262 case FLTN:{
263 int n;
264 double y1;
265
266 n = GET_INT(arg1);
267 x1 = (double) n;
268 y1 = GET_FLT(arg2);
269 return (makeflt(x1 + y1));
270 }
271 case LONGN:
272 return (bigx_plus
273 (bigx_int_to_big(arg1), bigx_long_to_big(arg2)));
274 case BIGX:
275 return (bigx_plus(bigx_int_to_big(arg1), arg2));
276 }
277 break;
278 case LONGN:
279 switch (tag2) {
280 case INTN:
281 return (bigx_plus
282 (bigx_long_to_big(arg1), bigx_int_to_big(arg2)));
283 case FLTN:
284 return (plus(exact_to_inexact(arg1), arg2));
285 case LONGN:
286 return (bigx_plus
287 (bigx_long_to_big(arg1), bigx_long_to_big(arg2)));
288 case BIGX:
289 return (bigx_plus(bigx_long_to_big(arg1), arg2));
290 }
291 break;
292 case BIGX:
293 switch (tag2) {
294 case INTN:
295 return (bigx_plus(arg1, bigx_int_to_big(arg2)));
296 case FLTN:
297 return (plus(exact_to_inexact(arg1), arg2));
298 case LONGN:
299 return (bigx_plus(arg1, bigx_long_to_big(arg2)));
300 case BIGX:
301 return (bigx_plus(arg1, arg2));
302 }
303 break;
304 case FLTN:
305 switch (tag2) {
306 case INTN:{
307 int s;
308
309 x1 = GET_FLT(arg1);
310 s = GET_INT(arg2);
311 x2 = (double) s;
312 return (makeflt(x1 + x2));
313 }
314 case FLTN:{
315 x1 = GET_FLT(arg1);
316 x2 = GET_FLT(arg2);
317 return (makeflt(x1 + x2));
318 }
319 case LONGN:
320 return (plus(arg1, exact_to_inexact(arg2)));
321 case BIGX:
322 return (plus(arg1, exact_to_inexact(arg2)));
323 }
324 break;
325 default:
326 break;
327 }
328 error(NOT_COMPUTABLE, "+", list2(arg1, arg2));
329 return (UNDEF);
330 }
331
332
333 int
minus(int arg1,int arg2)334 minus(int arg1, int arg2)
335 {
336 int tag1,
337 tag2;
338 double x1,
339 x2;
340
341 tag1 = GET_TAG(arg1);
342 tag2 = GET_TAG(arg2);
343 switch (tag1) {
344 case INTN:
345 switch (tag2) {
346 case INTN:{
347 long long int l;
348
349 l = (long long int) GET_INT(arg1) -
350 (long long int) GET_INT(arg2);
351 if (SMALL_INT_MIN < l && l < SMALL_INT_MAX)
352 return (makeint((int) l));
353 else
354 return (makelong(l));
355 }
356
357 case FLTN:{
358 int n;
359 double y1;
360
361 n = GET_INT(arg1);
362 x1 = (double) n;
363 y1 = GET_FLT(arg2);
364 return (makeflt(x1 - y1));
365 }
366
367 case LONGN:
368 return (bigx_minus(bigx_int_to_big(arg1),
369 bigx_long_to_big(arg2)));
370
371 case BIGX:
372 return (bigx_minus(bigx_int_to_big(arg1), arg2));
373 }
374 break;
375 case LONGN:
376 switch (tag2) {
377 case INTN:
378 return (bigx_minus
379 (bigx_long_to_big(arg1), bigx_int_to_big(arg2)));
380 case FLTN:
381 return (minus(exact_to_inexact(arg1), arg2));
382 case LONGN:
383 return (bigx_minus
384 (bigx_long_to_big(arg1), bigx_long_to_big(arg2)));
385 case BIGX:
386 return (bigx_minus(bigx_long_to_big(arg1), arg2));
387 }
388 break;
389 case BIGX:
390 switch (tag2) {
391 case INTN:
392 return (bigx_minus(arg1, bigx_int_to_big(arg2)));
393 case FLTN:
394 return (minus(exact_to_inexact(arg1), arg2));
395 case LONGN:
396 return (bigx_minus(arg1, bigx_long_to_big(arg2)));
397 case BIGX:
398 return (bigx_minus(arg1, arg2));
399 }
400 break;
401 case FLTN:
402 switch (tag2) {
403 case INTN:{
404 int s;
405
406 x1 = GET_FLT(arg1);
407 s = GET_INT(arg2);
408 x2 = (double) s;
409 return (makeflt(x1 - x2));
410 }
411 case FLTN:{
412 x1 = GET_FLT(arg1);
413 x2 = GET_FLT(arg2);
414 return (makeflt(x1 - x2));
415 }
416 case LONGN:
417 return (minus(arg1, exact_to_inexact(arg2)));
418 case BIGX:
419 return (minus(arg1, exact_to_inexact(arg2)));
420 }
421 break;
422 }
423 error(NOT_COMPUTABLE, "-", list2(arg1, arg2));
424 return (UNDEF);
425 }
426
427
428 int
mult(int arg1,int arg2)429 mult(int arg1, int arg2)
430 {
431 int tag1,
432 tag2;
433 double x1,
434 y1,
435 x2;
436
437 tag1 = GET_TAG(arg1);
438 tag2 = GET_TAG(arg2);
439
440 switch (tag1) {
441 case INTN:
442 switch (tag2) {
443 case INTN:{
444 long long int l1,
445 l2,
446 l;
447
448 l1 = (long long int) GET_INT(arg1);
449 l2 = (long long int) GET_INT(arg2);
450 l = l1 * l2;
451 if (l < SMALL_INT_MAX && l > SMALL_INT_MIN)
452 return (makeint((int) l));
453 else
454 return (bigx_mult
455 (bigx_int_to_big(arg1),
456 bigx_int_to_big(arg2)));
457 }
458
459 case FLTN:{
460 int n;
461
462 n = GET_INT(arg1);
463 x1 = (double) n;
464 y1 = GET_FLT(arg2);
465 return (makeflt(x1 * y1));
466 }
467
468 case LONGN:
469 if (GET_INT(arg1) != 0)
470 return (bigx_mult
471 (bigx_int_to_big(arg1), bigx_long_to_big(arg2)));
472 else
473 return (arg1); // int 0
474
475 case BIGX:
476 return (bigx_mult_i(arg2, arg1));
477 }
478 break;
479 case LONGN:
480 switch (tag2) {
481 case INTN:
482 if (GET_INT(arg2) != 0)
483 return (bigx_mult
484 (bigx_long_to_big(arg1), bigx_int_to_big(arg2)));
485 else
486 return (arg2); // int 0
487 case FLTN:
488 return (mult(exact_to_inexact(arg1), arg2));
489 case LONGN:
490 return (bigx_mult
491 (bigx_long_to_big(arg1), bigx_long_to_big(arg2)));
492 case BIGX:
493 return (bigx_mult(bigx_long_to_big(arg1), arg2));
494 }
495 break;
496 case BIGX:
497 switch (tag2) {
498 case INTN:
499 return (bigx_mult_i(arg1, arg2));
500 case FLTN:
501 return (mult(exact_to_inexact(arg1), arg2));
502 case LONGN:
503 return (bigx_mult(arg1, bigx_long_to_big(arg2)));
504 case BIGX:
505 return (bigx_mult(arg1, arg2));
506 }
507 break;
508 case FLTN:
509 switch (tag2) {
510 case INTN:{
511 int s;
512
513 x1 = GET_FLT(arg1);
514 s = GET_INT(arg2);
515 x2 = (double) s;
516 return (makeflt(x1 * x2));
517 }
518 case FLTN:{
519 x1 = GET_FLT(arg1);
520 x2 = GET_FLT(arg2);
521 y1 = x1 * x2;
522 if (y1 > DBL_MAX || y1 < -DBL_MAX)
523 error(FLT_OVERF, "*", list2(arg1, arg2));
524 if (x1 != 0.0 && x2 != 0.0 && y1 > -DBL_MIN
525 && y1 < DBL_MIN)
526 error(FLT_UNDERF, "*", list2(arg1, arg2));
527 return (makeflt(y1));
528 }
529 case LONGN:
530 case BIGX:
531 return (mult(arg1, exact_to_inexact(arg2)));
532 }
533 break;
534 }
535 error(NOT_COMPUTABLE, "*", list2(arg1, arg2));
536 return (UNDEF);
537 }
538
539 int
quotient(int arg1,int arg2)540 quotient(int arg1, int arg2)
541 {
542 int n,
543 s,
544 tag1,
545 tag2;
546 double x1,
547 y1,
548 x2;
549
550 tag1 = GET_TAG(arg1);
551 tag2 = GET_TAG(arg2);
552 switch (tag1) {
553 case INTN:
554 switch (tag2) {
555 case INTN:{
556 n = GET_INT(arg1);
557 if (n == 0)
558 return (arg1);
559 x1 = (double) n;
560 s = GET_INT(arg2);
561 y1 = (double) s;
562 x2 = x1 / y1;
563 if ((x2 - floor(x2)) == 0.0)
564 return (makeint((int) x2));
565 else
566 return (makeflt(x1 / y1));
567 }
568 case FLTN:{
569 n = GET_INT(arg1);
570 x1 = (double) n;
571 y1 = GET_FLT(arg2);
572 return (makeflt(x1 / y1));
573 }
574
575 case LONGN:
576 case BIGX:
577 if (GET_INT(arg1) == 0)
578 return (arg1);
579 else
580 return (quotient
581 (exact_to_inexact(arg1), exact_to_inexact(arg2)));
582 }
583 break;
584 case LONGN:
585 switch (tag2) {
586 case INTN:
587 case LONGN:
588 n = quotient(exact_to_inexact(arg1), exact_to_inexact(arg2));
589 x1 = GET_FLT(n);
590 x2 = x1 - ceil(x1);
591 if (x2 == 0.0) {
592 return (divide(arg1, arg2));
593 } else
594 return (n);
595 case BIGX:
596 return (quotient
597 (exact_to_inexact(arg1), exact_to_inexact(arg2)));
598 case FLTN:
599 return (quotient(exact_to_inexact(arg1), arg2));
600 }
601 break;
602 case BIGX:
603 switch (tag2) {
604 case INTN:
605 case LONGN:
606 case BIGX:
607 n = quotient(exact_to_inexact(arg1), exact_to_inexact(arg2));
608 x1 = GET_FLT(n);
609 x2 = x1 - ceil(x1);
610 if (x1 == 1.0)
611 return (makeint(1));
612 else if (x1 == -1.0)
613 return (makeint(-1));
614 else if (x2 == 0.0)
615 return (divide(arg1, arg2));
616 else
617 return (n);
618 case FLTN:
619 return (quotient(exact_to_inexact(arg1), arg2));
620 }
621 break;
622 case FLTN:
623 switch (tag2) {
624 case INTN:{
625 x1 = GET_FLT(arg1);
626 s = GET_INT(arg2);
627 x2 = (double) s;
628 return (makeflt(x1 / x2));
629 }
630 case FLTN:{
631 x1 = GET_FLT(arg1);
632 x2 = GET_FLT(arg2);
633 return (makeflt(x1 / x2));
634 }
635 case LONGN:
636 case BIGX:
637 return (quotient(arg1, exact_to_inexact(arg2)));
638
639 }
640 break;
641 }
642 error(NOT_COMPUTABLE, "/", list2(arg1, arg2));
643 return (UNDEF);
644 }
645
646
647 int
divide(int x,int y)648 divide(int x, int y)
649 {
650
651 if (integerp(x) && longnump(y))
652 return (makeint(0));
653 else if (integerp(x) && bignump(y))
654 return (makeint(0));
655 else if (longnump(x) && bignump(y))
656 return (makeint(0));
657 else if (integerp(x) && integerp(y))
658 return (makeint(GET_INT(x) / GET_INT(y)));
659 else if (longnump(x) && integerp(y))
660 return (long_int_div(x, y));
661 else if (longnump(x) && longnump(y))
662 return (long_long_div(x, y));
663 else if (bignump(x) && integerp(y))
664 return (bigx_div_i(x, y));
665 else if (bignump(x) && longnump(y))
666 return (bigx_div(x, bigx_long_to_big(y)));
667 else if (bignump(x) && bignump(y))
668 return (bigx_div(x, y));
669 else
670 error(ILLEGAL_ARGS, "div", list2(x, y));
671
672 return (UNDEF);
673 }
674
675 int
s_remainder(int x,int y)676 s_remainder(int x, int y)
677 {
678 int i;
679
680 if (integerp(x) && integerp(y))
681 return (makeint(GET_INT(x) % GET_INT(y)));
682 else if (integerp(x) && longnump(y))
683 return (x);
684 else if (integerp(x) && bignump(y))
685 return (x);
686 else if (longnump(x) && bignump(y))
687 return (x);
688 else if (longnump(x) && integerp(y))
689 return (long_int_remainder(x, y));
690 else if (longnump(x) && longnump(y))
691 return (long_long_remainder(x, y));
692 else if (bignump(x) && integerp(y))
693 return (bigx_remainder_i(x, y));
694 else if (bignump(x) && longnump(y)) {
695 i = bigx_long_to_big(y);
696 return (minus(x, mult(divide(x, i), i)));
697 } else if (bignump(x) && bignump(y))
698 return (minus(x, mult(divide(x, y), y)));
699
700 error(ILLEGAL_ARGS, "remainder", NIL);
701 return (UNDEF);
702 }
703
704
705 // remainder of longnum and int.
706 int
long_int_remainder(int x,int y)707 long_int_remainder(int x, int y)
708 {
709 long long int m,
710 n,
711 r;
712
713
714 m = GET_LONG(x);
715 n = (long long int) GET_INT(y);
716 r = m % n;
717
718 return (makeint((int) r));
719 }
720
721 // remainder of longnum and longnum.
722 int
long_long_remainder(int x,int y)723 long_long_remainder(int x, int y)
724 {
725 long long int m,
726 n,
727 r;
728
729
730 m = GET_LONG(x);
731 n = GET_LONG(y);
732 r = m % n;
733
734 if (llabs(r) < BIGNUM_BASE)
735 return (makeint((int) r));
736 else
737 return (makelong(r));
738 }
739
740 // divide of longnum and int.
741 int
long_int_div(int x,int y)742 long_int_div(int x, int y)
743 {
744 long long int m,
745 n,
746 q;
747
748
749 m = GET_LONG(x);
750 n = (long long int) GET_INT(y);
751 q = m / n;
752
753 if (llabs(q) < BIGNUM_BASE)
754 return (makeint((int) q));
755 else
756 return (makelong(q));
757 }
758
759 // divide of longnum and longnum
760 int
long_long_div(int x,int y)761 long_long_div(int x, int y)
762 {
763 long long int m,
764 n,
765 q;
766
767
768 m = GET_LONG(x);
769 n = GET_LONG(y);
770 q = m / n;
771
772 if (llabs(q) < BIGNUM_BASE)
773 return (makeint((int) q));
774 else
775 return (makelong(q));
776 }
777
778 int
absolute(int x)779 absolute(int x)
780 {
781 if (integerp(x))
782 return (makeint(abs(GET_INT(x))));
783 else if (longnump(x))
784 return (makelong(llabs(GET_LONG(x))));
785 else if (bignump(x)) {
786 return (bigx_abs(x));
787 } else if (floatp(x)) {
788 return (makeflt(fabs(GET_FLT(x))));
789 }
790 error(ILLEGAL_ARGS, "abs", x);
791 return (UNDEF);
792 }
793
794 int
angle(int y,int x)795 angle(int y, int x)
796 {
797
798 if (positive_zerop(y) && positivep(x))
799 return (makeflt(0.0));
800 else if (negative_zerop(y) && positivep(x))
801 return (makeflt(-0.0));
802 else if (positivep(y) && positivep(x))
803 return (makeflt
804 (atan
805 (GET_FLT(exact_to_inexact(y)) /
806 GET_FLT(exact_to_inexact(x)))));
807 else if (positivep(y) && zerop(x))
808 return (makeflt(M_PI_2));
809 else if (positivep(y) && negativep(x))
810 return (makeflt
811 (M_PI +
812 atan(GET_FLT(exact_to_inexact(y)) /
813 GET_FLT(exact_to_inexact(x)))));
814 else if (positive_zerop(y) && negativep(x))
815 return (makeflt(M_PI));
816 else if (negative_zerop(y) && negativep(x))
817 return (makeflt(-M_PI));
818 else if (negativep(y) && negativep(x))
819 return (makeflt
820 (-M_PI +
821 atan(GET_FLT(exact_to_inexact(y)) /
822 GET_FLT(exact_to_inexact(x)))));
823 else if (negativep(y) && zerop(x))
824 return (makeflt(-M_PI_2));
825 else if (negativep(y) && positivep(x))
826 return (makeflt
827 (atan
828 (GET_FLT(exact_to_inexact(y)) /
829 GET_FLT(exact_to_inexact(x)))));
830 else if (zerop(y) && positivep(x))
831 return (makeflt(0.0));
832 else if (zerop(y) && negativep(x))
833 return (makeflt(M_PI));
834 else if (positive_zerop(y) && positive_zerop(x))
835 return (makeflt(+0.0));
836 else if (negative_zerop(y) && positive_zerop(x))
837 return (makeflt(-0.0));
838 else if (positive_zerop(y) && negative_zerop(x))
839 return (makeflt(M_PI));
840 else if (negative_zerop(y) && negative_zerop(x))
841 return (makeflt(-M_PI));
842 else if (positive_zerop(y) && zerop(x))
843 return (makeflt(M_PI_2));
844 else if (negative_zerop(y) && zerop(x))
845 return (makeflt(-M_PI_2));
846 else {
847 error(ILLEGAL_ARGS, "angle", list2(x, y));
848 return (UNDEF);
849 }
850 }
851
852 // GCD of integer
853 int
int_gcd(int x,int y)854 int_gcd(int x, int y)
855 {
856 if (y == 0)
857 return (x);
858
859 while (y != 0) {
860 int r;
861
862 r = x % y;
863 x = y;
864 y = r;
865 }
866 return (x);
867 }
868
869 // GCD of number(include bignum)
870 int
gcd(int x,int y)871 gcd(int x, int y)
872 {
873 if (integerp(x) && integerp(y))
874 return (makeint(abs(int_gcd(GET_INT(x), GET_INT(y)))));
875
876 else if (floatp(x) && integerp(y))
877 return (makeflt
878 ((double) abs(int_gcd((int) GET_FLT(x), GET_INT(y)))));
879
880 else if (integerp(x) && floatp(y))
881 return (makeflt
882 ((double) abs(int_gcd(GET_INT(x), (int) GET_FLT(y)))));
883
884 else if (floatp(x) && floatp(y))
885 return (makeflt
886 ((double) abs(int_gcd(GET_FLT(x), (int) GET_FLT(y)))));
887
888 while (!zerop(y)) {
889 int r;
890
891 r = s_remainder(x, y);
892 x = y;
893 y = r;
894
895 }
896 return (absolute(x));
897 }
898
899 // LCM of integer
900 int
int_lcm(int m,int n)901 int_lcm(int m, int n)
902 {
903 if (m == 0 || n == 0)
904 return (0);
905
906 return ((m / int_gcd(m, n)) * n);
907 }
908
909 // LCM of number(include bignum)
910 int
lcm(int x,int y)911 lcm(int x, int y)
912 {
913 int g,
914 d,
915 res;
916 if (integerp(x) && integerp(y) && abs(GET_INT(x)) < 10000 && abs(GET_INT(y)) < 10000) // because
917 //
918 //
919 //
920 //
921 //
922 //
923 //
924 //
925 //
926 // x,y
927 // <
928 // sqrt(BIGNUM_BASE)
929 return (makeint(abs(int_lcm(GET_INT(x), GET_INT(y)))));
930
931 else if (floatp(x) && integerp(y))
932 return (makeflt
933 ((double) abs(int_lcm((int) GET_FLT(x), GET_INT(y)))));
934
935 else if (integerp(x) && floatp(y))
936 return (makeflt
937 ((double) abs(int_lcm(GET_INT(x), (int) GET_FLT(y)))));
938
939 else if (floatp(x) && floatp(y))
940 return (makeflt
941 ((double) abs(int_lcm(GET_FLT(x), (int) GET_FLT(y)))));
942
943 else {
944 g = gcd(x, y);
945 d = divide(absolute(x), g);
946 res = mult(d, absolute(y));
947 return (res);
948 }
949 }
950
951 int
isqrt(int x)952 isqrt(int x)
953 {
954 if (integerp(x))
955 return (makeint(floor(sqrt(GET_INT(x)))));
956 else if (floatp(x))
957 return (makeint(floor(sqrt(GET_FLT(x)))));
958 else
959 return (isqrt1(makeint(1), makeint(1), x));
960 }
961
962 int
isqrt1(int s,int s2,int x)963 isqrt1(int s, int s2, int x)
964 {
965 if (eqsmallerp(mult(s, s), x) && eqsmallerp(x, mult(s2, s2)))
966 return (s);
967 else
968 return (isqrt1(divide(plus(divide(x, s), s), makeint(2)), s, x));
969 }
970
971 /*
972 * try and error for OPEN-MP int mat_mult(int x, int y){ int
973 * i,j,k,m,n,o,p,m1,n1,o1,p1,res; double
974 * A[MATSIZE][MATSIZE],B[MATSIZE][MATSIZE],C[MATSIZE][MATSIZE];
975 *
976 * if(!(length(array_length(x)) == 2 && length(array_length(y)) == 2))
977 * error(ILLEGAL_ARGS, "*", list2(x,y)); if(GET_INT(cadr(array_length(x)))
978 * != GET_INT(car(array_length(y)))) error(ILLEGAL_ARGS, "*", list2(x,y));
979 *
980 * m1 = car(array_length(x)); n1 = cadr(array_length(x)); o1 =
981 * car(array_length(y)); p1 = cadr(array_length(y));
982 *
983 * m = GET_INT(m1); n = GET_INT(n1); o = GET_INT(o1); p = GET_INT(p1);
984 * if(m>MATSIZE || n>MATSIZE || o>MATSIZE || p>MATSIZE)
985 * error(OUT_OF_RANGE,"*",NIL); if(m*p+FREESIZE > fc){ shelterpush(x);
986 * shelterpush(y); gbc(); shelterpop(); shelterpop(); } res =
987 * makearray(list2(m1,p1),UNDEF); #pragma omp parallel for private(i,j)
988 * for(i=0; i<m; i++) for(j=0; j<n; j++) A[i][j] =
989 * GET_FLT(matrix_ref(x,n,i,j)); #pragma omp parallel for private(i,j)
990 * for(i=0; i<o; i++) for(j=0; j<p; j++) B[i][j] =
991 * GET_FLT(matrix_ref(y,p,i,j)); #pragma omp parallel for private(i,j)
992 * for(i=0; i<m; i++) for(j=0; j<p; j++) C[i][j] = 0.0; #pragma omp
993 * parallel for private(i,j,k) for(i=0; i<m; i++) for(j=0; j<p; j++)
994 * for(k=0; k<n; k++) C[i][j] = C[i][j] + A[i][k] * B[k][j];
995 *
996 *
997 * for(i=0; i<m; i++) for(j=0; j<p; j++)
998 * matrix_set(res,p,i,j,makeflt(C[i][j]));
999 *
1000 * return(NIL); }
1001 *
1002 * int mat_vec_mult(int x, int y){ int i,j,m,n,res; float
1003 * A[MATSIZE][MATSIZE],B[MATSIZE],C[MATSIZE];
1004 *
1005 * if(length(array_length(x)) != 2) error(ILLEGAL_ARGS, "*", x);
1006 * if(GET_INT(cadr(array_length(x))) != vector_length(y))
1007 * error(ILLEGAL_ARGS, "*", list2(x,y));
1008 *
1009 * m = GET_INT(car(array_length(x))); n = GET_INT(cadr(array_length(x)));
1010 * if(m>MATSIZE || n>MATSIZE) error(OUT_OF_RANGE,"*",NIL); if(m+FREESIZE >
1011 * fc) gbc(); res = makevec(m,UNDEF); #pragma omp parallel for
1012 * private(i,j) for(i=0; i<m; i++) for(j=0; j<n; j++) A[i][j] =
1013 * GET_FLT(matrix_ref(x,n,i,j)); #pragma omp parallel for private(i)
1014 * for(i=0; i<m; i++) B[i] = GET_FLT(vector_ref(y,i)); #pragma omp
1015 * parallel for private(i) for(i=0; i<m; i++) C[i] = 0.0; #pragma omp
1016 * parallel for private(i,j) for(i=0; i<m; i++) for(j=0; j<n; j++) C[i] =
1017 * C[i] + A[i][j] * B[j];
1018 *
1019 * for(i=0; i<m; i++) vector_set(res,i,makeflt(C[i]));
1020 *
1021 * return(res); }
1022 *
1023 * int vec_plus(int x, int y){ int i,size,res; double
1024 * A[MATSIZE],B[MATSIZE],C[MATSIZE];
1025 *
1026 * if(vector_length(x) != vector_length(y)) error(ILLEGAL_ARGS,"+",
1027 * list2(x,y));
1028 *
1029 * size = vector_length(x); res = makevec(size,UNDEF);
1030 *
1031 * #pragma omp parallel for private(i) for(i=0; i<size; i++) A[i] =
1032 * GET_FLT(vector_ref(x,i)); #pragma omp parallel for private(i) for(i=0;
1033 * i<size; i++) B[i] = GET_FLT(vector_ref(y,i)); #pragma omp parallel for
1034 * private(i) for(i=0; i<size; i++) C[i] = 0.0; #pragma omp parallel for
1035 * private(i) for(i=0; i<size; i++) C[i] = A[i] + B[i];
1036 *
1037 * for(i=0; i<size; i++) vector_set(res,i,makeflt(C[i]));
1038 *
1039 * return(res); }
1040 *
1041 * int vec_minus(int x, int y){ int i,size,res; double
1042 * A[MATSIZE],B[MATSIZE],C[MATSIZE];
1043 *
1044 * if(vector_length(x) != vector_length(y)) error(ILLEGAL_ARGS,"-",
1045 * list2(x,y));
1046 *
1047 * size = vector_length(x); res = makevec(size,UNDEF);
1048 *
1049 * #pragma omp parallel for private(i) for(i=0; i<size; i++) A[i] =
1050 * GET_FLT(vector_ref(x,i)); #pragma omp parallel for private(i) for(i=0;
1051 * i<size; i++) B[i] = GET_FLT(vector_ref(y,i)); #pragma omp parallel for
1052 * private(i) for(i=0; i<size; i++) C[i] = 0.0; #pragma omp parallel for
1053 * private(i) for(i=0; i<size; i++) C[i] = A[i] - B[i];
1054 *
1055 * for(i=0; i<size; i++) vector_set(res,i,makeflt(C[i]));
1056 *
1057 * return(res); }
1058 */
1059