123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296 |
- #include "stdafx.h"
- //-----------------------------------------------------------------------------
- //
- // Generate all divisors of a term
- //
- // Input: Term on stack (factor * factor * ...)
- //
- // Output: Divisors on stack
- //
- //-----------------------------------------------------------------------------
- #include "defs.h"
- static void gen(int, int);
- static void __factor_add(void);
- static int __cmp(const void *, const void *);
- void
- divisors(void)
- {
- int i, h, n;
- save();
- h = tos - 1;
- divisors_onstack();
- n = tos - h;
- qsort(stack + h, n, sizeof (U *), __cmp);
- p1 = alloc_tensor(n);
- p1->u.tensor->ndim = 1;
- p1->u.tensor->dim[0] = n;
- for (i = 0; i < n; i++)
- p1->u.tensor->elem[i] = stack[h + i];
- tos = h;
- push(p1);
- restore();
- }
- void
- divisors_onstack(void)
- {
- int h, i, k, n;
- save();
- p1 = pop();
- h = tos;
- // push all of the term's factors
- if (isnum(p1)) {
- push(p1);
- factor_small_number();
- } else if (car(p1) == symbol(ADD)) {
- push(p1);
- __factor_add();
- //printf(">>>");
- //for (i = h; i < tos; i++)
- //print(stdout, stack[i]);
- //printf("<<<");
- } else if (car(p1) == symbol(MULTIPLY)) {
- p1 = cdr(p1);
- if (isnum(car(p1))) {
- push(car(p1));
- factor_small_number();
- p1 = cdr(p1);
- }
- while (iscons(p1)) {
- p2 = car(p1);
- if (car(p2) == symbol(POWER)) {
- push(cadr(p2));
- push(caddr(p2));
- } else {
- push(p2);
- push(one);
- }
- p1 = cdr(p1);
- }
- } else if (car(p1) == symbol(POWER)) {
- push(cadr(p1));
- push(caddr(p1));
- } else {
- push(p1);
- push(one);
- }
- k = tos;
- // contruct divisors by recursive descent
- push(one);
- gen(h, k);
- // move
- n = tos - k;
- for (i = 0; i < n; i++)
- stack[h + i] = stack[k + i];
- tos = h + n;
- restore();
- }
- //-----------------------------------------------------------------------------
- //
- // Generate divisors
- //
- // Input: Base-exponent pairs on stack
- //
- // h first pair
- //
- // k just past last pair
- //
- // Output: Divisors on stack
- //
- // For example, factor list 2 2 3 1 results in 6 divisors,
- //
- // 1
- // 3
- // 2
- // 6
- // 4
- // 12
- //
- //-----------------------------------------------------------------------------
- #define ACCUM p1
- #define BASE p2
- #define EXPO p3
- static void
- gen(int h, int k)
- {
- int expo, i;
- save();
- ACCUM = pop();
- if (h == k) {
- push(ACCUM);
- restore();
- return;
- }
- BASE = stack[h + 0];
- EXPO = stack[h + 1];
- push(EXPO);
- expo = pop_integer();
- for (i = 0; i <= abs(expo); i++) {
- push(ACCUM);
- push(BASE);
- push_integer(sign(expo) * i);
- power();
- multiply();
- gen(h + 2, k);
- }
- restore();
- }
- //-----------------------------------------------------------------------------
- //
- // Factor ADD expression
- //
- // Input: Expression on stack
- //
- // Output: Factors on stack
- //
- // Each factor consists of two expressions, the factor itself followed
- // by the exponent.
- //
- //-----------------------------------------------------------------------------
- static void
- __factor_add(void)
- {
- save();
- p1 = pop();
- // get gcd of all terms
- p3 = cdr(p1);
- push(car(p3));
- p3 = cdr(p3);
- while (iscons(p3)) {
- push(car(p3));
- gcd();
- p3 = cdr(p3);
- }
- // check gcd
- p2 = pop();
- if (isplusone(p2)) {
- push(p1);
- push(one);
- restore();
- return;
- }
- // push factored gcd
- if (isnum(p2)) {
- push(p2);
- factor_small_number();
- } else if (car(p2) == symbol(MULTIPLY)) {
- p3 = cdr(p2);
- if (isnum(car(p3))) {
- push(car(p3));
- factor_small_number();
- } else {
- push(car(p3));
- push(one);
- }
- p3 = cdr(p3);
- while (iscons(p3)) {
- push(car(p3));
- push(one);
- p3 = cdr(p3);
- }
- } else {
- push(p2);
- push(one);
- }
- // divide each term by gcd
- push(p2);
- inverse();
- p2 = pop();
- push(zero);
- p3 = cdr(p1);
- while (iscons(p3)) {
- push(p2);
- push(car(p3));
- multiply();
- add();
- p3 = cdr(p3);
- }
- push(one);
- restore();
- }
- static int
- __cmp(const void *p1, const void *p2)
- {
- return cmp_expr(*((U **) p1), *((U **) p2));
- }
- #if SELFTEST
- static char *s[] = {
- "divisors(12)",
- "(1,2,3,4,6,12)",
- "divisors(-12)",
- "(1,2,3,4,6,12)",
- "divisors(a)",
- "(1,a)",
- "divisors(-a)",
- "(1,a)",
- "divisors(+3*x+3)",
- "(1,3,1+x,3+3*x)",
- "divisors(+3*x-3)",
- "(1,3,-3+3*x,-1+x)",
- "divisors(-3*x+3)",
- "(1,3,1-x,3-3*x)",
- "divisors(-3*x-3)",
- "(1,3,1+x,3+3*x)",
- };
- void
- test_divisors(void)
- {
- test(__FILE__, s, sizeof s / sizeof (char *));
- }
- #endif
|