divisors.cpp

   1 #include "stdafx.h"
   2
   3 //-----------------------------------------------------------------------------
   4 //
   5 //      Generate all divisors of a term
   6 //
   7 //      Input:          Term on stack (factor * factor * ...)
   8 //
   9 //      Output:         Divisors on stack
  10 //
  11 //-----------------------------------------------------------------------------
  12
  13 #include "defs.h"
  14
  15 static void gen(int, int);
  16 static void __factor_add(void);
  17 static int __cmp(const void *, const void *);
  18
  19 void
  20 divisors(void)
  21 {
  22         int i, h, n;
  23         save();
  24         h = tos - 1;
  25         divisors_onstack();
  26         n = tos - h;
  27         qsort(stack + h, n, sizeof (U *), __cmp);
  28         p1 = alloc_tensor(n);
  29         p1->u.tensor->ndim = 1;
  30         p1->u.tensor->dim[0] = n;
  31         for (i = 0; i < n; i++)
  32                 p1->u.tensor->elem[i] = stack[h + i];
  33         tos = h;
  34         push(p1);
  35         restore();
  36 }
  37
  38 void
  39 divisors_onstack(void)
  40 {
  41         int h, i, k, n;
  42
  43         save();
  44
  45         p1 = pop();
  46
  47         h = tos;
  48
  49         // push all of the term's factors
  50
  51         if (isnum(p1)) {
  52                 push(p1);
  53                 factor_small_number();
  54         } else if (car(p1) == symbol(ADD)) {
  55                 push(p1);
  56                 __factor_add();
  57 //printf(">>>");
  58 //for (i = h; i < tos; i++)
  59 //print(stdout, stack[i]);
  60 //printf("<<<");
  61         } else if (car(p1) == symbol(MULTIPLY)) {
  62                 p1 = cdr(p1);
  63                 if (isnum(car(p1))) {
  64                         push(car(p1));
  65                         factor_small_number();
  66                         p1 = cdr(p1);
  67                 }
  68                 while (iscons(p1)) {
  69                         p2 = car(p1);
  70                         if (car(p2) == symbol(POWER)) {
  71                                 push(cadr(p2));
  72                                 push(caddr(p2));
  73                         } else {
  74                                 push(p2);
  75                                 push(one);
  76                         }
  77                         p1 = cdr(p1);
  78                 }
  79         } else if (car(p1) == symbol(POWER)) {
  80                 push(cadr(p1));
  81                 push(caddr(p1));
  82         } else {
  83                 push(p1);
  84                 push(one);
  85         }
  86
  87         k = tos;
  88
  89         // contruct divisors by recursive descent
  90
  91         push(one);
  92
  93         gen(h, k);
  94
  95         // move
  96
  97         n = tos - k;
  98
  99         for (i = 0; i < n; i++)
 100                 stack[h + i] = stack[k + i];
 101
 102         tos = h + n;
 103
 104         restore();
 105 }
 106
 107 //-----------------------------------------------------------------------------
 108 //
 109 //      Generate divisors
 110 //
 111 //      Input:          Base-exponent pairs on stack
 112 //
 113 //                      h       first pair
 114 //
 115 //                      k       just past last pair
 116 //
 117 //      Output:         Divisors on stack
 118 //
 119 //      For example, factor list 2 2 3 1 results in 6 divisors,
 120 //
 121 //              1
 122 //              3
 123 //              2
 124 //              6
 125 //              4
 126 //              12
 127 //
 128 //-----------------------------------------------------------------------------
 129
 130 #define ACCUM p1
 131 #define BASE p2
 132 #define EXPO p3
 133
 134 static void
 135 gen(int h, int k)
 136 {
 137         int expo, i;
 138
 139         save();
 140
 141         ACCUM = pop();
 142
 143         if (h == k) {
 144                 push(ACCUM);
 145                 restore();
 146                 return;
 147         }
 148
 149         BASE = stack[h + 0];
 150         EXPO = stack[h + 1];
 151
 152         push(EXPO);
 153         expo = pop_integer();
 154
 155         for (i = 0; i <= abs(expo); i++) {
 156                 push(ACCUM);
 157                 push(BASE);
 158                 push_integer(sign(expo) * i);
 159                 power();
 160                 multiply();
 161                 gen(h + 2, k);
 162         }
 163
 164         restore();
 165 }
 166
 167 //-----------------------------------------------------------------------------
 168 //
 169 //      Factor ADD expression
 170 //
 171 //      Input:          Expression on stack
 172 //
 173 //      Output:         Factors on stack
 174 //
 175 //      Each factor consists of two expressions, the factor itself followed
 176 //      by the exponent.
 177 //
 178 //-----------------------------------------------------------------------------
 179
 180 static void
 181 __factor_add(void)
 182 {
 183         save();
 184
 185         p1 = pop();
 186
 187         // get gcd of all terms
 188
 189         p3 = cdr(p1);
 190         push(car(p3));
 191         p3 = cdr(p3);
 192         while (iscons(p3)) {
 193                 push(car(p3));
 194                 gcd();
 195                 p3 = cdr(p3);
 196         }
 197
 198         // check gcd
 199
 200         p2 = pop();
 201         if (isplusone(p2)) {
 202                 push(p1);
 203                 push(one);
 204                 restore();
 205                 return;
 206         }
 207
 208         // push factored gcd
 209
 210         if (isnum(p2)) {
 211                 push(p2);
 212                 factor_small_number();
 213         } else if (car(p2) == symbol(MULTIPLY)) {
 214                 p3 = cdr(p2);
 215                 if (isnum(car(p3))) {
 216                         push(car(p3));
 217                         factor_small_number();
 218                 } else {
 219                         push(car(p3));
 220                         push(one);
 221                 }
 222                 p3 = cdr(p3);
 223                 while (iscons(p3)) {
 224                         push(car(p3));
 225                         push(one);
 226                         p3 = cdr(p3);
 227                 }
 228         } else {
 229                 push(p2);
 230                 push(one);
 231         }
 232
 233         // divide each term by gcd
 234
 235         push(p2);
 236         inverse();
 237         p2 = pop();
 238
 239         push(zero);
 240         p3 = cdr(p1);
 241         while (iscons(p3)) {
 242                 push(p2);
 243                 push(car(p3));
 244                 multiply();
 245                 add();
 246                 p3 = cdr(p3);
 247         }
 248
 249         push(one);
 250
 251         restore();
 252 }
 253
 254 static int
 255 __cmp(const void *p1, const void *p2)
 256 {
 257         return cmp_expr(*((U **) p1), *((U **) p2));
 258 }
 259
 260 #if SELFTEST
 261
 262 static char *s[] = {
 263
 264         "divisors(12)",
 265         "(1,2,3,4,6,12)",
 266
 267         "divisors(-12)",
 268         "(1,2,3,4,6,12)",
 269
 270         "divisors(a)",
 271         "(1,a)",
 272
 273         "divisors(-a)",
 274         "(1,a)",
 275
 276         "divisors(+3*x+3)",
 277         "(1,3,1+x,3+3*x)",
 278
 279         "divisors(+3*x-3)",
 280         "(1,3,-3+3*x,-1+x)",
 281
 282         "divisors(-3*x+3)",
 283         "(1,3,1-x,3-3*x)",
 284
 285         "divisors(-3*x-3)",
 286         "(1,3,1+x,3+3*x)",
 287 };
 288
 289 void
 290 test_divisors(void)
 291 {
 292         test(__FILE__, s, sizeof s / sizeof (char *));
 293 }
 294
 295 #endif