1module simpfact; 2 3% Redistribution and use in source and binary forms, with or without 4% modification, are permitted provided that the following conditions are met: 5% 6% * Redistributions of source code must retain the relevant copyright 7% notice, this list of conditions and the following disclaimer. 8% * Redistributions in binary form must reproduce the above copyright 9% notice, this list of conditions and the following disclaimer in the 10% documentation and/or other materials provided with the distribution. 11% 12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 16% CONTRIBUTORS 17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 23% POSSIBILITY OF SUCH DAMAGE. 24% 25 26 27% Simplification for quotients containing factorials 28 29% Matthew Rebbeck ( while in placement at ZIB) - March 1994. 30 31% The new 'really' improved version! Simplifies plain factorials as 32% well as those raised to integer powers and 1/2. 33% 34% Deals properly with the generalised factorial idea of simplifying 35% non integers, eg: (k+1/2)!/(k-1/2)! -> k+1/2. 36 37algebraic << 38operator simplify_factorial1; 39operator simplify_factorial; 40operator int_simplify_factorial; 41 42let simplify_factorial(~x) => simplify_factorial1(num x,den x); 43 44let { simplify_factorial1(~x,~z) => int_simplify_factorial(x/z)}; 45 46let { simplify_factorial1 (~x + ~y,~z) => 47 simplify_factorial1 (x,z) + simplify_factorial1(y,z)}; 48>>; 49 50symbolic procedure int_simplify_factorial (u); 51 begin 52 scalar minus_num,minus_denom,test_expt; 53 if not pairp u or car u neq 'quotient then u 54 else 55 << 56 % 57 % We firstly produce input of standard form. 58 % 59 if atom cadr u or atom caddr u then u 60 else << 61 % 62 % Remove 'minus if there. 63 % 64 if car cadr u eq 'minus then 65 << cadr u := cadr cadr u; 66 minus_num := t; >>; 67 if car caddr u eq 'minus then 68 << caddr u := cadr caddr u; 69 minus_denom := t; >>; 70 71 if car cadr u eq 'factorial then cadr u := {'times,cadr u}; 72 if car caddr u eq 'factorial then caddr u := {'times,caddr u}; 73 if car cadr u eq 'oddexpt or car cadr u eq 'expt 74 or car cadr u eq 'sqrt 75 then cadr u := {'times,cadr u}; 76 if car caddr u eq 'oddexpt or car caddr u eq 'expt or 77 car caddr u eq 'sqrt 78 then caddr u := {'times,caddr u}; 79 % 80 % Test to see if input contains any 'expt's. If it does 81 % then they are converted to 'oddexpts and re converted 82 % at the end. If not (ie: either contains 'oddexpt's or 83 % no powers at all), then no conversion is done and the 84 % output is left in this oddexpt form. 85 % 86 if test_for_expt(cadr u) or test_for_expt(caddr u) then 87 << test_expt := t; 88 convert_to_oddexpt(cadr u); 89 convert_to_oddexpt(caddr u); >>; 90 91 if test_for_facts(cadr u,caddr u) 92 then gothru_numerator(cadr u,caddr u); 93 94 if minus_num then cadr u := {'minus,cadr u}; 95 if minus_denom then caddr u := {'minus,caddr u}; 96 97 cadr u := reval cadr u; 98 caddr u := reval caddr u; >>; 99 % 100 % Output converted back to 'expt form regardless of the form 101 % of the input. For this conversion to occur only if input 102 % is in 'expt form (perhaps useful with Wolfram's input) 103 % then uncomment next line... 104 %if test_expt then 105 u := algebraic sub(oddexpt=expt,u); 106 >>; 107 return u; 108 end; 109 110flag('(int_simplify_factorial),'opfn); 111 112 113 114symbolic procedure test_for_expt(input); 115 % 116 % Tests to see if 'expt occurs anywhere. 117 % 118 begin 119 scalar found_expt,not_found; 120 not_found := t; 121 while input and not_found do 122 << 123 if pairp car input and (caar input = 'expt or caar input = 'sqrt) 124 then <<found_expt:=t; not_found:=nil;>>; 125 input := cdr input; 126 >>; 127 return found_expt; 128 end; 129 130flag('(test_for_expt),'boolean); 131 132 133 134symbolic procedure convert_to_oddexpt(input); 135 % 136 % Converts all expt's to standard form. ie: oddexpt(......,power). 137 % 138 begin 139 while input do 140 << 141 if pairp car input and caar input = 'expt 142 then caar input := 'oddexpt; 143 if pairp car input and caar input = 'sqrt then 144 << 145 caar input := 'oddexpt; 146 cdar input := {cadar input,{'quotient,1,2}}; 147 >>; 148 input := cdr input; 149 >>; 150 end; 151 152 153 154symbolic procedure gothru_numerator(num,denom); 155 % 156 % Go systematically through numerator, searching for factorials, and, 157 % when found, comparing with denominator. 'change' describes if 158 % simplifications have been made or not (ie:change eq 0). 159 % 160 begin 161 scalar change,orignum,origdenom; 162 change := 0; 163 orignum := num; 164 origdenom := denom; 165 % 166 % while in numerator. 167 % 168 while num do 169 << 170 if pairp car num and caar num eq 'oddexpt then 171 << 172 if pairp cadar num and caadar num eq 'factorial then 173 change := change + gothru_denominator(num,denom); 174 >> 175 else if pairp car num and caar num eq 'factorial then 176 << 177 change := change + gothru_denominator(num,denom); 178 >>; 179 num := cdr num; 180 >>; 181 % 182 % If at end of numerator but simplifications have been made, 183 % then repeat. 184 % 185 if not num and not eqn(change,0) then 186 << 187 if test_for_facts(orignum,origdenom) 188 then gothru_numerator(orignum,origdenom); % Beginning. 189 >>; 190 end; 191 192 193 194symbolic procedure gothru_denominator(num,denom); 195 % 196 % Systematically goes through denominator finding factorials and 197 % passing numerator and denom. factorials into oddexpt_test. There 198 % they are simplified if possible. 'Compared' describes if the 199 % factorials were simplified (ie: car test eq ok) or if it was not 200 % possible. 201 % 202 begin 203 scalar test,change; 204 change := 0; 205 206 while denom and change = 0 do 207 << 208 if pairp car denom and caar denom eq 'oddexpt then 209 << 210 if pairp cadar denom and caadar denom eq 'factorial then 211 << 212 test := oddexpt_test(num,denom,change); 213 change := change + test; 214 >>; 215 >> 216 else if pairp car denom and caar denom eq 'factorial then 217 << 218 test := oddexpt_test(num,denom,change); 219 change := change + test; 220 >>; 221 denom := cdr denom; 222 >>; 223 return change; 224 end; 225 226 227 228symbolic procedure oddexpt_test(num,denom,change); 229 % 230 % Tests which parts of quotient, (if any), are exponentials, passing 231 % the quotient onto the relevant simplifying function. 232 % 233 begin 234 scalar test; 235 236 if caar num eq 'oddexpt and caar denom neq 'oddexpt then 237 << test := compare_numoddexptfactorial(num,denom,change); >> 238 else if caar num neq 'oddexpt and caar denom eq 'oddexpt then 239 << test := compare_denomoddexptfactorial(num,denom,change); >> 240 else if caar num eq 'oddexpt and caar denom eq 'oddexpt then 241 << test := compare_bothoddexptfactorial(num,denom,change); >> 242 else test := compare_factorial(num,denom,change); 243 244 return test; 245 end; 246 247 248 249symbolic procedure compare_factorial (num,denom,change); 250 % 251 % Compares factorials, simplifying if possible. 252 % 253 begin 254 scalar numsimp,denomsimp,diff; 255 256 % If both factorial arguments are of the same form. 257 if numberp (reval list('difference,cadar (num),cadar(denom))) then 258 << 259 change := change + 1; 260 % Difference between num. and denom. factorial arguments. 261 diff :=(reval list('difference,cadar (num),cadar(denom))); 262 % If argument of num. factorial > argument of denom. factorial. 263 if diff >0 then 264 << 265 % numsimp collects simplified numerator arguments. 266 numsimp := for i := 1:diff collect reval {'plus,cadar denom,i}; 267 % Remove num. factorial and replace with simplification. 268 car num := 'times.numsimp; 269 % Remove denom. factorial. 270 car denom := 1; 271 >> 272 else % if diff <= 0 then 273 << 274 diff := -diff; 275 denomsimp := for i := 1:diff collect reval {'plus,cadar num,i}; 276 car denom := 'times.denomsimp; 277 car num := 1; 278 >>; 279 >>; 280 return change; 281 end; 282 283 284 285symbolic procedure compare_numoddexptfactorial (num,denom,change); 286 % 287 % Compares factorials with oddexpt num., simplifying if possible.See 288 % compare_factorial for more detailed comments. 289 % 290 begin 291 scalar diff; 292 293 if numberp (reval list('difference,car cdadar num,cadar denom)) 294 then 295 << 296 % New sqrt additions... 297 if sqrt_test(num) then 298 << 299 << 300 diff :=(reval list('difference, car cdadar num,cadar denom)); 301 change := change+1; 302 if diff > 0 then simplify_sqrt1(num,denom,diff) 303 else simplify_sqrt2(num,denom,diff); 304 >>; 305 >> 306 % If power is not integer or 1/2 then can't simplify. 307 else if not_int_or_sqrt(num) then <<>> 308 % If oddexpt. of power 2. 309 else if eqn(caddar num-1,1) then 310 << 311 % Remove oddexpt. 312 car num := car {cadar num}; 313 diff := (reval list('difference,cadar num,cadar denom)); 314 change := change +1; 315 if diff > 0 then << simplify1(num,denom,diff); >> 316 else simplify2(num,denom,diff); 317 >> 318 else 319 << 320 % Reduce oddexpt by one. 321 car num := {caar num,cadar num,car cddar num -1}; 322 diff :=(reval list('difference,car cdadar num,cadar denom)); 323 change := change + 1; 324 if diff >0 then << simplify1(num,denom,diff); >> 325 else simplify2(cdar num,denom,diff); 326 >>; 327 >>; 328 return change; 329 end; 330 331 332 333symbolic procedure simplify_sqrt1(num,denom,diff); 334 begin 335 scalar numsimp; 336 numsimp := for i := 1:diff collect reval {'plus,cadar denom,i}; 337 cadar num := car{'times.numsimp}; 338 car denom := {'oddexpt,car denom,{'quotient,1,2}}; 339 end; 340 341 342 343symbolic procedure simplify_sqrt2(num,denom,diff); 344 begin 345 scalar denomsimp; 346 diff := -diff; 347 denomsimp := for i := 1:diff collect reval {'plus,car cdadar num,i}; 348 car denom := reval {'times,car num,car{'times.denomsimp}}; 349 car num := 1; 350 end; 351 352 353 354symbolic procedure simplify1(num,denom,diff); 355 begin 356 scalar numsimp; 357 numsimp := for i := 1:diff collect reval {'plus,cadar denom,i}; 358 cdr num := car{'times.numsimp}.cdr num; 359 car denom := 1; 360 end; 361 362 363 364symbolic procedure simplify2(num,denom,diff); 365 begin 366 scalar denomsimp; 367 diff := -diff; 368 denomsimp := for i := 1:diff collect reval {'plus,cadar num,i}; 369 cdr denom := car{'times.denomsimp}.cdr denom; 370 car denom := 1; 371 end; 372 373 374 375symbolic procedure compare_denomoddexptfactorial (num,denom,change); 376 % 377 % Compares factorials with oddexpt denom, simplifying if possible.See 378 % compare_factorial and compare_numoddexptfactorial for more detailed 379 % comments. 380 % 381 begin 382 scalar change,diff; 383 384 if numberp (reval list('difference, cadar num,car cdadar denom)) 385 then 386 << 387 % New sqrt additions... 388 if sqrt_test(denom) then 389 << 390 << 391 diff :=(reval list('difference, cadar num,car cdadar denom)); 392 change := change+1; 393 if diff > 0 then simplify_sqrt3(num,denom,diff) 394 else % if diff <= 0 395 simplify_sqrt4(num,denom,diff); 396 >>; 397 >> 398 else if not_int_or_sqrt(denom) then <<>> 399 else if eqn(caddar denom-1,1) then 400 << 401 car denom := car {cadar denom}; 402 diff := (reval list('difference,cadar num,cadar denom)); 403 change := change +1; 404 if diff > 0 then simplify3(num,denom,diff) 405 else % if diff <= 0 then 406 simplify4(num,denom,diff); 407 >> 408 else 409 << 410 car denom := {caar denom,cadar denom,car cddar denom -1}; 411 diff :=(reval list('difference, cadar num,car cdadar denom)); 412 change := change + 1; 413 if diff >0 then simplify3(num,cdar denom,diff) 414 else simplify4(num,denom,diff); 415 >>; 416 >>; 417 return change; 418 end; 419 420 421 422symbolic procedure sqrt_test(input); 423 % 424 % tests if the expt power is 1/2. (boolean) 425 % 426 begin 427 if caddar input = '(quotient 1 2) then return t 428 else return nil; 429 end; 430 431flag('(sqrt_test),'boolean); 432 433 434 435symbolic procedure not_int_or_sqrt(input); 436 % 437 % tests if the expt power is neither int or 1/2. (boolean) 438 % 439 begin 440 if pairp caddar input and car caddar input = 'quotient 441 and cdr caddar input neq '(1 2) then return t 442 else return nil; 443 end; 444 445flag('(not_int_or_sqrt),'boolean); 446 447 448 449symbolic procedure simplify_sqrt3(num,denom,diff); 450 begin 451 scalar numsimp; 452 numsimp := for i := 1:diff collect reval {'plus,car cdadar denom,i}; 453 car num := reval{'times,car denom,car{'times.numsimp}}; 454 car denom := 1; 455 end; 456 457 458 459symbolic procedure simplify_sqrt4(num,denom,diff); 460 begin 461 scalar denomsimp; 462 diff := -diff; 463 denomsimp := for i := 1:diff collect reval {'plus,cadar num,i}; 464 if diff = 0 then car denom := 1 465 else cadar denom := car{'times.denomsimp}; 466 car num := {'oddexpt,car num,{'quotient,1,2}}; 467 end; 468 469 470 471symbolic procedure simplify3(num,denom,diff); 472 begin 473 scalar numsimp; 474 numsimp := for i := 1:diff collect reval {'plus,cadar denom,i}; 475 cdr num := car{'times.numsimp}.cdr num; 476 car num := 1; 477 end; 478 479 480 481symbolic procedure simplify4(num,denom,diff); 482 begin 483 scalar denomsimp; 484 diff := -diff; 485 denomsimp := for i := 1:diff collect reval {'plus,cadar num,i}; 486 cdr denom := car{'times.denomsimp}.cdr denom; 487 car num := 1; 488 end; 489 490 491 492symbolic procedure compare_bothoddexptfactorial (num,denom,change); 493 % 494 % Compares factorials with both oddexpt num. & denom., simplifying if 495 % possible. See previous compare_...... functions for more detailed 496 % comments. 497 % 498 begin 499 scalar change,diff; 500 501 if numberp(reval list('difference,car cdadar num,car cdadar denom)) 502 then 503 << 504 % New sqrt additions... 505 if sqrt_test(num) and sqrt_test(denom) then 506 << 507 << 508 diff :=(reval list('difference, car cdadar num,car cdadar 509 denom)); 510 change := change+1; 511 if diff > 0 then simplify_sqrt5(num,denom,diff) 512 else % if diff <= 0 513 simplify_sqrt6(num,denom,diff); 514 >>; 515 >> 516 else if not_int_or_sqrt(num) or not_int_or_sqrt(denom) then <<>> 517 % If denom is sqrt but num is not. 518 else if sqrt_test(denom) then 519 << 520 diff := reval list('difference,cadr cadar num,cadr cadar denom); 521 if diff > 0 then simplify_sqrt5(num,denom,diff) 522 else % if diff <= 0 then 523 simplify_sqrt6(num,denom,diff); 524 >> 525 % If num is sqrt but denom is not. 526 else if sqrt_test(num) then 527 << 528 diff := reval list('difference,cadr cadar num,cadr cadar denom); 529 if diff > 0 then simplify_sqrt7(num,denom,diff) 530 else % if diff <= 0 then 531 simplify_sqrt8(num,denom,diff); 532 >> 533 else if eqn(caddar num-1,1) and eqn(caddar denom-1,1) then 534 << 535 car num := car {cadar num}; 536 car denom := car {cadar denom}; 537 diff := (reval list('difference,cadar num,cadar denom)); 538 change := change +1; 539 if diff > 0 then simplify5(num,denom,diff) 540 else % if diff <= 0 then 541 simplify6(num,denom,diff); 542 >> 543 else if eqn(caddar num-1,1) and not eqn(caddar denom-1,1) then 544 << 545 car num := car {cadar num}; 546 car denom := {caar denom,cadar denom,car cddar denom-1}; 547 diff := (reval list('difference,cadar num,car cdadar denom)); 548 change := change +1; 549 if diff >0 then simplify5(num,cdar denom,diff) 550 else % if diff <= 0 then 551 simplify6(num,denom,diff); 552 >> 553 else if caddar num-1 neq 1 and caddar denom-1 eq 1 then 554 << 555 car num := {caar num,cadar num,car cddar num-1}; 556 car denom := car {cadar denom}; 557 diff := (reval list('difference,car cdadar num,cadar denom)); 558 change := change +1; 559 if diff >0 then simplify5(num,denom,diff) 560 else simplify6(cdar num,denom,diff); 561 >> 562 else if caddar num-1 neq 1 and caddar denom-1 neq 1 then 563 << 564 car num := {caar num,cadar num,car cddar num-1}; 565 car denom := {caar denom,cadar denom,car cddar denom-1}; 566 diff:=(reval list('difference,car cdadar num,car cdadar denom)); 567 change := change +1; 568 if diff >0 then simplify5(num,cdar denom,diff) 569 else simplify6(cdar num,denom,diff); 570 >>; 571 >>; 572 return change; 573 end; 574 575 576 577symbolic procedure simplify_sqrt5(num,denom,diff); 578 begin 579 scalar numsimp; 580 numsimp := for i := 1:diff collect reval {'plus,car cdadar denom,i}; 581 car num := {'times,{'oddexpt,cadar denom,{'plus,caddar num, 582 {'minus,{'quotient,1,2}}}},{'oddexpt,car{'times.numsimp}, 583 caddar num}}; 584 car denom := 1; 585 end; 586 587 588 589symbolic procedure simplify_sqrt6(num,denom,diff); 590 begin 591 scalar denomsimp; 592 diff := -diff; 593 denomsimp := for i := 1:diff collect reval {'plus,car cdadar num,i}; 594 car denom := {'oddexpt,car{'times.denomsimp},{'quotient,1,2}}; 595 caddar num := {'plus,caddar num,{'minus,{'quotient,1,2}}}; 596 end; 597 598 599 600symbolic procedure simplify_sqrt7(num,denom,diff); 601 begin 602 scalar numsimp; 603 numsimp := for i := 1:diff collect reval {'plus,car cdadar denom,i}; 604 car num := {'oddexpt,car{'times.numsimp},{'quotient,1,2}}; 605 caddar denom := {'plus,caddar denom,{'minus,{'quotient,1,2}}}; 606 end; 607 608 609 610symbolic procedure simplify_sqrt8(num,denom,diff); 611 begin 612 scalar denomsimp; 613 diff := -diff; 614 denomsimp := for i := 1:diff collect reval {'plus,car cdadar num,i}; 615 car denom:= {'times,{'oddexpt, cadar num,{'plus,caddar denom, 616 {'minus,{'quotient,1,2}}}},{'oddexpt,car{'times.denomsimp}, 617 caddar denom}}; 618 car num := 1; 619 end; 620 621 622 623symbolic procedure simplify5(num,denom,diff); 624 begin 625 scalar numsimp; 626 numsimp := for i := 1:diff collect reval {'plus,cadar denom,i}; 627 cdr num := car{'times.numsimp}.cdr num; 628 end; 629 630 631 632symbolic procedure simplify6(num,denom,diff); 633 begin 634 scalar denomsimp; 635 diff := -diff; 636 denomsimp := for i := 1:diff collect reval {'plus,cadar num,i}; 637 cdr denom := car{'times.denomsimp}.cdr denom; 638 end; 639 640 641 642symbolic procedure test_for_facts(num,denom); 643 % 644 % Systematically goes through numerator and then denom. looking for 645 % factorials. 646 % (boolean). 647 % 648 begin 649 scalar test; 650 if test_num(num) and test_denom(denom) then test := t; 651 return test 652 end; 653 654flag('(test_for_facts),'boolean); 655 656 657 658symbolic procedure test_num(num); 659 % 660 % Systematically goes through num., looking for factorials. 661 % (boolean). 662 % 663 begin 664 scalar test; 665 test := nil; 666 if eqcar (num ,'times) or eqcar (num ,'oddexpt) then 667 while num and not test do 668 << 669 if pairp car num and caar num eq 'factorial then test := t 670 else if pairp car num and caar num eq 'oddexpt 671 then if pairp cadar num and caadar num eq 'factorial 672 then test := t; 673 num := cdr num; 674 >>; 675 return test; 676 end; 677 678flag ('(test_num),'boolean); 679 680 681 682symbolic procedure test_denom(denom); 683 % 684 % Systematically goes through denominator, looking for factorials. 685 % (boolean). 686 % 687 begin 688 scalar test; 689 test := nil; 690 if eqcar (denom ,'times) or eqcar (denom ,'oddexpt) then 691 while denom and not test do 692 << 693 if pairp car denom and caar denom eq 'factorial then test := t 694 else if pairp car denom and caar denom eq 'oddexpt 695 then if pairp cadar denom and caadar denom eq 'factorial 696 then test := t; 697 denom:= cdr denom; 698 >>; 699 return test; 700 end; 701 702flag ('(test_denom),'boolean); 703 704endmodule; 705 706end; 707 708