gcd.cpp

   1 // Greatest common denominator
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_gcd(void)
   8 {
   9         p1 = cdr(p1);
  10         push(car(p1));
  11         eval();
  12         p1 = cdr(p1);
  13         while (iscons(p1)) {
  14                 push(car(p1));
  15                 eval();
  16                 gcd();
  17                 p1 = cdr(p1);
  18         }
  19 }
  20
  21 void
  22 gcd(void)
  23 {
  24         int x = expanding;
  25         save();
  26         gcd_main();
  27         restore();
  28         expanding = x;
  29 }
  30
  31 void
  32 gcd_main(void)
  33 {
  34         expanding = 1;
  35
  36         p2 = pop();
  37         p1 = pop();
  38
  39         if (equal(p1, p2)) {
  40                 push(p1);
  41                 return;
  42         }
  43
  44         if (isrational(p1) && isrational(p2)) {
  45                 push(p1);
  46                 push(p2);
  47                 gcd_numbers();
  48                 return;
  49         }
  50
  51         if (car(p1) == symbol(ADD) && car(p2) == symbol(ADD)) {
  52                 gcd_expr_expr();
  53                 return;
  54         }
  55
  56         if (car(p1) == symbol(ADD)) {
  57                 gcd_expr(p1);
  58                 p1 = pop();
  59         }
  60
  61         if (car(p2) == symbol(ADD)) {
  62                 gcd_expr(p2);
  63                 p2 = pop();
  64         }
  65
  66         if (car(p1) == symbol(MULTIPLY) && car(p2) == symbol(MULTIPLY)) {
  67                 gcd_term_term();
  68                 return;
  69         }
  70
  71         if (car(p1) == symbol(MULTIPLY)) {
  72                 gcd_term_factor();
  73                 return;
  74         }
  75
  76         if (car(p2) == symbol(MULTIPLY)) {
  77                 gcd_factor_term();
  78                 return;
  79         }
  80
  81         // gcd of factors
  82
  83         if (car(p1) == symbol(POWER)) {
  84                 p3 = caddr(p1);
  85                 p1 = cadr(p1);
  86         } else
  87                 p3 = one;
  88
  89         if (car(p2) == symbol(POWER)) {
  90                 p4 = caddr(p2);
  91                 p2 = cadr(p2);
  92         } else
  93                 p4 = one;
  94
  95         if (!equal(p1, p2)) {
  96                 push(one);
  97                 return;
  98         }
  99
 100         // are both exponents numerical?
 101
 102         if (isnum(p3) && isnum(p4)) {
 103                 push(p1);
 104                 if (lessp(p3, p4))
 105                         push(p3);
 106                 else
 107                         push(p4);
 108                 power();
 109                 return;
 110         }
 111
 112         // are the exponents multiples of eah other?
 113
 114         push(p3);
 115         push(p4);
 116         divide();
 117
 118         p5 = pop();
 119
 120         if (isnum(p5)) {
 121
 122                 push(p1);
 123
 124                 // choose the smallest exponent
 125
 126                 if (car(p3) == symbol(MULTIPLY) && isnum(cadr(p3)))
 127                         p5 = cadr(p3);
 128                 else
 129                         p5 = one;
 130
 131                 if (car(p4) == symbol(MULTIPLY) && isnum(cadr(p4)))
 132                         p6 = cadr(p4);
 133                 else
 134                         p6 = one;
 135
 136                 if (lessp(p5, p6))
 137                         push(p3);
 138                 else
 139                         push(p4);
 140
 141                 power();
 142                 return;
 143         }
 144
 145         push(p3);
 146         push(p4);
 147         subtract();
 148
 149         p5 = pop();
 150
 151         if (!isnum(p5)) {
 152                 push(one);
 153                 return;
 154         }
 155
 156         // can't be equal because of test near beginning
 157
 158         push(p1);
 159
 160         if (isnegativenumber(p5))
 161                 push(p3);
 162         else
 163                 push(p4);
 164
 165         power();
 166 }
 167
 168 // in this case gcd is used as a composite function, i.e. gcd(gcd(gcd...
 169
 170 void
 171 gcd_expr_expr(void)
 172 {
 173         if (length(p1) != length(p2)) {
 174                 push(one);
 175                 return;
 176         }
 177
 178         p3 = cdr(p1);
 179         push(car(p3));
 180         p3 = cdr(p3);
 181         while (iscons(p3)) {
 182                 push(car(p3));
 183                 gcd();
 184                 p3 = cdr(p3);
 185         }
 186         p3 = pop();
 187
 188         p4 = cdr(p2);
 189         push(car(p4));
 190         p4 = cdr(p4);
 191         while (iscons(p4)) {
 192                 push(car(p4));
 193                 gcd();
 194                 p4 = cdr(p4);
 195         }
 196         p4 = pop();
 197
 198         push(p1);
 199         push(p3);
 200         divide();
 201         p5 = pop();
 202
 203         push(p2);
 204         push(p4);
 205         divide();
 206         p6 = pop();
 207
 208         if (equal(p5, p6)) {
 209                 push(p5);
 210                 push(p3);
 211                 push(p4);
 212                 gcd();
 213                 multiply();
 214         } else
 215                 push(one);
 216 }
 217
 218 void
 219 gcd_expr(U *p)
 220 {
 221         p = cdr(p);
 222         push(car(p));
 223         p = cdr(p);
 224         while (iscons(p)) {
 225                 push(car(p));
 226                 gcd();
 227                 p = cdr(p);
 228         }
 229 }
 230
 231 void
 232 gcd_term_term(void)
 233 {
 234         push(one);
 235         p3 = cdr(p1);
 236         while (iscons(p3)) {
 237                 p4 = cdr(p2);
 238                 while (iscons(p4)) {
 239                         push(car(p3));
 240                         push(car(p4));
 241                         gcd();
 242                         multiply();
 243                         p4 = cdr(p4);
 244                 }
 245                 p3 = cdr(p3);
 246         }
 247 }
 248
 249 void
 250 gcd_term_factor(void)
 251 {
 252         push(one);
 253         p3 = cdr(p1);
 254         while (iscons(p3)) {
 255                 push(car(p3));
 256                 push(p2);
 257                 gcd();
 258                 multiply();
 259                 p3 = cdr(p3);
 260         }
 261 }
 262
 263 void
 264 gcd_factor_term(void)
 265 {
 266         push(one);
 267         p4 = cdr(p2);
 268         while (iscons(p4)) {
 269                 push(p1);
 270                 push(car(p4));
 271                 gcd();
 272                 multiply();
 273                 p4 = cdr(p4);
 274         }
 275 }
 276
 277 #if SELFTEST
 278
 279 static char *s[] = {
 280
 281         "gcd(30,42)",
 282         "6",
 283
 284         "gcd(42,30)",
 285         "6",
 286
 287         "gcd(-30,42)",
 288         "6",
 289
 290         "gcd(42,-30)",
 291         "6",
 292
 293         "gcd(30,-42)",
 294         "6",
 295
 296         "gcd(-42,30)",
 297         "6",
 298
 299         "gcd(-30,-42)",
 300         "6",
 301
 302         "gcd(-42,-30)",
 303         "6",
 304
 305         "gcd(x,x)",
 306         "x",
 307
 308         "gcd(-x,x)",
 309         "x",
 310
 311         "gcd(x,-x)",
 312         "x",
 313
 314         "gcd(-x,-x)",
 315         "-x",
 316
 317         "gcd(x^2,x^3)",
 318         "x^2",
 319
 320         "gcd(x,y)",
 321         "1",
 322
 323         "gcd(y,x)",
 324         "1",
 325
 326         "gcd(x*y,y)",
 327         "y",
 328
 329         "gcd(x*y,y^2)",
 330         "y",
 331
 332         "gcd(x^2*y^2,x^3*y^3)",
 333         "x^2*y^2",
 334
 335         "gcd(x^2,x^3)",
 336         "x^2",
 337
 338         // gcd of expr
 339
 340         "gcd(x+y,x+z)",
 341         "1",
 342
 343         "gcd(x+y,x+y)",
 344         "x+y",
 345
 346         "gcd(x+y,2*x+2*y)",
 347         "x+y",
 348
 349         "gcd(-x-y,x+y)",
 350         "x+y",
 351
 352         "gcd(4*x+4*y,6*x+6*y)",
 353         "2*x+2*y",
 354
 355         "gcd(4*x+4*y+4,6*x+6*y+6)",
 356         "2+2*x+2*y",
 357
 358         "gcd(4*x+4*y+4,6*x+6*y+12)",
 359         "1",
 360
 361         "gcd(27*t^3+y^3+9*t*y^2+27*t^2*y,t+y)",
 362         "1",
 363
 364         // gcd expr factor
 365
 366         "gcd(2*a^2*x^2+a*x+a*b,a)",
 367         "a",
 368
 369         "gcd(2*a^2*x^2+a*x+a*b,a^2)",
 370         "a",
 371
 372         "gcd(2*a^2*x^2+2*a*x+2*a*b,a)",
 373         "a",
 374
 375         // gcd expr term
 376
 377         "gcd(2*a^2*x^2+2*a*x+2*a*b,2*a)",
 378         "2*a",
 379
 380         "gcd(2*a^2*x^2+2*a*x+2*a*b,3*a)",
 381         "a",
 382
 383         "gcd(2*a^2*x^2+2*a*x+2*a*b,4*a)",
 384         "2*a",
 385
 386         // gcd factor factor
 387
 388         "gcd(x,x^2)",
 389         "x",
 390
 391         "gcd(x,x^a)",
 392         "1",
 393
 394         // multiple arguments
 395
 396         "gcd(12,18,9)",
 397         "3",
 398 };
 399
 400 void
 401 test_gcd(void)
 402 {
 403         test(__FILE__, s, sizeof s / sizeof (char *));
 404 }
 405
 406 #endif