1if lisp !*rounded then rounded_was_on := t 2 else rounded_was_on := nil; 3 4 5 6mat1 := mat((1,2,3,4,5),(2,3,4,5,6),(3,4,5,6,7),(4,5,6,7,8),(5,6,7,8,9)); 7 8 9 [1 2 3 4 5] 10 [ ] 11 [2 3 4 5 6] 12 [ ] 13mat1 := [3 4 5 6 7] 14 [ ] 15 [4 5 6 7 8] 16 [ ] 17 [5 6 7 8 9] 18 19 20mat2 := mat((1,1,1,1),(2,2,2,2),(3,3,3,3),(4,4,4,4)); 21 22 23 [1 1 1 1] 24 [ ] 25 [2 2 2 2] 26mat2 := [ ] 27 [3 3 3 3] 28 [ ] 29 [4 4 4 4] 30 31 32mat3 := mat((x),(x),(x),(x)); 33 34 35 [x] 36 [ ] 37 [x] 38mat3 := [ ] 39 [x] 40 [ ] 41 [x] 42 43 44mat4 := mat((3,3),(4,4),(5,5),(6,6)); 45 46 47 [3 3] 48 [ ] 49 [4 4] 50mat4 := [ ] 51 [5 5] 52 [ ] 53 [6 6] 54 55 56mat5 := mat((1,2,1,1),(1,2,3,1),(4,5,1,2),(3,4,5,6)); 57 58 59 [1 2 1 1] 60 [ ] 61 [1 2 3 1] 62mat5 := [ ] 63 [4 5 1 2] 64 [ ] 65 [3 4 5 6] 66 67 68mat6 := mat((i+1,i+2,i+3),(4,5,2),(1,i,0)); 69 70 71 [i + 1 i + 2 i + 3] 72 [ ] 73mat6 := [ 4 5 2 ] 74 [ ] 75 [ 1 i 0 ] 76 77 78mat7 := mat((1,1,0),(1,3,1),(0,1,1)); 79 80 81 [1 1 0] 82 [ ] 83mat7 := [1 3 1] 84 [ ] 85 [0 1 1] 86 87 88mat8 := mat((1,3),(-4,3)); 89 90 91 [1 3] 92mat8 := [ ] 93 [-4 3] 94 95 96mat9 := mat((1,2,3,4),(9,8,7,6)); 97 98 99 [1 2 3 4] 100mat9 := [ ] 101 [9 8 7 6] 102 103 104poly := x^7+x^5+4*x^4+5*x^3+12; 105 106 107 7 5 4 3 108poly := x + x + 4*x + 5*x + 12 109 110poly1 := x^2+x*y^3+x*y*z^3+y*x+2+y*3; 111 112 113 2 3 3 114poly1 := x + x*y + x*y*z + x*y + 3*y + 2 115 116 117on errcont; 118 119 120 121 122% Basis matrix manipulations. 123 124add_columns(mat1,1,2,5*y); 125 126 127[1 5*y + 2 3 4 5] 128[ ] 129[2 10*y + 3 4 5 6] 130[ ] 131[3 15*y + 4 5 6 7] 132[ ] 133[4 5*(4*y + 1) 6 7 8] 134[ ] 135[5 25*y + 6 7 8 9] 136 137 138add_rows(mat1,1,2,x); 139 140 141[ 1 2 3 4 5 ] 142[ ] 143[x + 2 2*x + 3 3*x + 4 4*x + 5 5*x + 6] 144[ ] 145[ 3 4 5 6 7 ] 146[ ] 147[ 4 5 6 7 8 ] 148[ ] 149[ 5 6 7 8 9 ] 150 151 152 153add_to_columns(mat1,3,1000); 154 155 156[1 2 1003 4 5] 157[ ] 158[2 3 1004 5 6] 159[ ] 160[3 4 1005 6 7] 161[ ] 162[4 5 1006 7 8] 163[ ] 164[5 6 1007 8 9] 165 166 167add_to_columns(mat1,{1,2,3},y); 168 169 170[y + 1 y + 2 y + 3 4 5] 171[ ] 172[y + 2 y + 3 y + 4 5 6] 173[ ] 174[y + 3 y + 4 y + 5 6 7] 175[ ] 176[y + 4 y + 5 y + 6 7 8] 177[ ] 178[y + 5 y + 6 y + 7 8 9] 179 180 181add_to_rows(mat1,2,1000); 182 183 184[ 1 2 3 4 5 ] 185[ ] 186[1002 1003 1004 1005 1006] 187[ ] 188[ 3 4 5 6 7 ] 189[ ] 190[ 4 5 6 7 8 ] 191[ ] 192[ 5 6 7 8 9 ] 193 194 195add_to_rows(mat1,{1,2,3},x); 196 197 198[x + 1 x + 2 x + 3 x + 4 x + 5] 199[ ] 200[x + 2 x + 3 x + 4 x + 5 x + 6] 201[ ] 202[x + 3 x + 4 x + 5 x + 6 x + 7] 203[ ] 204[ 4 5 6 7 8 ] 205[ ] 206[ 5 6 7 8 9 ] 207 208 209 210augment_columns(mat1,2); 211 212 213[2] 214[ ] 215[3] 216[ ] 217[4] 218[ ] 219[5] 220[ ] 221[6] 222 223 224augment_columns(mat1,{1,2,5}); 225 226 227[1 2 5] 228[ ] 229[2 3 6] 230[ ] 231[3 4 7] 232[ ] 233[4 5 8] 234[ ] 235[5 6 9] 236 237 238stack_rows(mat1,3); 239 240 241[3 4 5 6 7] 242 243 244stack_rows(mat1,{1,3,5}); 245 246 247[1 2 3 4 5] 248[ ] 249[3 4 5 6 7] 250[ ] 251[5 6 7 8 9] 252 253 254 255char_poly(mat1,x); 256 257 258 3 2 259x *(x - 25*x - 50) 260 261 262column_dim(mat2); 263 264 2654 266 267row_dim(mat1); 268 269 2705 271 272 273copy_into(mat7,mat1,2,3); 274 275 276[1 2 3 4 5] 277[ ] 278[2 3 1 1 0] 279[ ] 280[3 4 1 3 1] 281[ ] 282[4 5 0 1 1] 283[ ] 284[5 6 7 8 9] 285 286 287copy_into(mat7,mat1,5,5); 288 289***** Error in copy_into: the matrix 290 291[1 1 0] 292[ ] 293[1 3 1] 294[ ] 295[0 1 1] 296 297 does not fit into 298 299[1 2 3 4 5] 300[ ] 301[2 3 4 5 6] 302[ ] 303[3 4 5 6 7] 304[ ] 305[4 5 6 7 8] 306[ ] 307[5 6 7 8 9] 308 309 at position 5,5. 310 311 312diagonal(3); 313 314 315[3] 316 317 318% diagonal can take both a list of arguments or just the arguments. 319diagonal({mat2,mat6}); 320 321 322[1 1 1 1 0 0 0 ] 323[ ] 324[2 2 2 2 0 0 0 ] 325[ ] 326[3 3 3 3 0 0 0 ] 327[ ] 328[4 4 4 4 0 0 0 ] 329[ ] 330[0 0 0 0 i + 1 i + 2 i + 3] 331[ ] 332[0 0 0 0 4 5 2 ] 333[ ] 334[0 0 0 0 1 i 0 ] 335 336 337diagonal(mat1,mat2,mat5); 338 339 340[1 2 3 4 5 0 0 0 0 0 0 0 0] 341[ ] 342[2 3 4 5 6 0 0 0 0 0 0 0 0] 343[ ] 344[3 4 5 6 7 0 0 0 0 0 0 0 0] 345[ ] 346[4 5 6 7 8 0 0 0 0 0 0 0 0] 347[ ] 348[5 6 7 8 9 0 0 0 0 0 0 0 0] 349[ ] 350[0 0 0 0 0 1 1 1 1 0 0 0 0] 351[ ] 352[0 0 0 0 0 2 2 2 2 0 0 0 0] 353[ ] 354[0 0 0 0 0 3 3 3 3 0 0 0 0] 355[ ] 356[0 0 0 0 0 4 4 4 4 0 0 0 0] 357[ ] 358[0 0 0 0 0 0 0 0 0 1 2 1 1] 359[ ] 360[0 0 0 0 0 0 0 0 0 1 2 3 1] 361[ ] 362[0 0 0 0 0 0 0 0 0 4 5 1 2] 363[ ] 364[0 0 0 0 0 0 0 0 0 3 4 5 6] 365 366 367 368extend(mat1,3,2,x); 369 370 371[1 2 3 4 5 x x] 372[ ] 373[2 3 4 5 6 x x] 374[ ] 375[3 4 5 6 7 x x] 376[ ] 377[4 5 6 7 8 x x] 378[ ] 379[5 6 7 8 9 x x] 380[ ] 381[x x x x x x x] 382[ ] 383[x x x x x x x] 384[ ] 385[x x x x x x x] 386 387 388 389find_companion(mat5,x); 390 391 392 2 393x - 2*x - 2 394 395 396get_columns(mat1,1); 397 398 399{ 400 401 [1] 402 [ ] 403 [2] 404 [ ] 405 [3] 406 [ ] 407 [4] 408 [ ] 409 [5] 410 411 } 412 413get_columns(mat1,{1,2}); 414 415 416{ 417 418 [1] 419 [ ] 420 [2] 421 [ ] 422 [3] 423 [ ] 424 [4] 425 [ ] 426 [5] 427 428 , 429 430 [2] 431 [ ] 432 [3] 433 [ ] 434 [4] 435 [ ] 436 [5] 437 [ ] 438 [6] 439 440 } 441 442get_rows(mat1,3); 443 444 445{ 446 447 [3 4 5 6 7] 448 449 } 450 451get_rows(mat1,{1,3}); 452 453 454{ 455 456 [1 2 3 4 5] 457 458 , 459 460 [3 4 5 6 7] 461 462 } 463 464 465hermitian_tp(mat6); 466 467 468[ - i + 1 4 1 ] 469[ ] 470[ - i + 2 5 - i] 471[ ] 472[ - i + 3 2 0 ] 473 474 475 476% matrix_augment and matrix_stack can take both a list of arguments 477% or just the arguments. 478matrix_augment({mat1,mat2}); 479 480***** Error in matrix_augment: 481***** all input matrices must have the same row dimension. 482 483matrix_augment(mat4,mat2,mat4); 484 485 486[3 3 1 1 1 1 3 3] 487[ ] 488[4 4 2 2 2 2 4 4] 489[ ] 490[5 5 3 3 3 3 5 5] 491[ ] 492[6 6 4 4 4 4 6 6] 493 494 495matrix_stack(mat1,mat2); 496 497***** Error in matrix_stack: 498***** all input matrices must have the same column dimension. 499 500matrix_stack({mat6,mat((z,z,z)),mat7}); 501 502 503[i + 1 i + 2 i + 3] 504[ ] 505[ 4 5 2 ] 506[ ] 507[ 1 i 0 ] 508[ ] 509[ z z z ] 510[ ] 511[ 1 1 0 ] 512[ ] 513[ 1 3 1 ] 514[ ] 515[ 0 1 1 ] 516 517 518 519minor(mat1,2,3); 520 521 522[1 2 4 5] 523[ ] 524[3 4 6 7] 525[ ] 526[4 5 7 8] 527[ ] 528[5 6 8 9] 529 530 531 532mult_columns(mat1,3,y); 533 534 535[1 2 3*y 4 5] 536[ ] 537[2 3 4*y 5 6] 538[ ] 539[3 4 5*y 6 7] 540[ ] 541[4 5 6*y 7 8] 542[ ] 543[5 6 7*y 8 9] 544 545 546mult_columns(mat1,{2,3,4},100); 547 548 549[1 200 300 400 5] 550[ ] 551[2 300 400 500 6] 552[ ] 553[3 400 500 600 7] 554[ ] 555[4 500 600 700 8] 556[ ] 557[5 600 700 800 9] 558 559 560mult_rows(mat1,2,x); 561 562 563[ 1 2 3 4 5 ] 564[ ] 565[2*x 3*x 4*x 5*x 6*x] 566[ ] 567[ 3 4 5 6 7 ] 568[ ] 569[ 4 5 6 7 8 ] 570[ ] 571[ 5 6 7 8 9 ] 572 573 574mult_rows(mat1,{1,3,5},10); 575 576 577[10 20 30 40 50] 578[ ] 579[2 3 4 5 6 ] 580[ ] 581[30 40 50 60 70] 582[ ] 583[4 5 6 7 8 ] 584[ ] 585[50 60 70 80 90] 586 587 588 589pivot(mat1,3,3); 590 591 592[ - 4 - 2 2 4 ] 593[------ ------ 0 --- --- ] 594[ 5 5 5 5 ] 595[ ] 596[ - 2 - 1 1 2 ] 597[------ ------ 0 --- --- ] 598[ 5 5 5 5 ] 599[ ] 600[ 3 4 5 6 7 ] 601[ ] 602[ 2 1 - 1 - 2 ] 603[ --- --- 0 ------ ------] 604[ 5 5 5 5 ] 605[ ] 606[ 4 2 - 2 - 4 ] 607[ --- --- 0 ------ ------] 608[ 5 5 5 5 ] 609 610 611rows_pivot(mat1,3,3,{1,5}); 612 613 614[ - 4 - 2 2 4 ] 615[------ ------ 0 --- --- ] 616[ 5 5 5 5 ] 617[ ] 618[ 2 3 4 5 6 ] 619[ ] 620[ 3 4 5 6 7 ] 621[ ] 622[ 4 5 6 7 8 ] 623[ ] 624[ 4 2 - 2 - 4 ] 625[ --- --- 0 ------ ------] 626[ 5 5 5 5 ] 627 628 629 630remove_columns(mat1,3); 631 632 633[1 2 4 5] 634[ ] 635[2 3 5 6] 636[ ] 637[3 4 6 7] 638[ ] 639[4 5 7 8] 640[ ] 641[5 6 8 9] 642 643 644remove_columns(mat1,{2,3,4}); 645 646 647[1 5] 648[ ] 649[2 6] 650[ ] 651[3 7] 652[ ] 653[4 8] 654[ ] 655[5 9] 656 657 658remove_rows(mat1,2); 659 660 661[1 2 3 4 5] 662[ ] 663[3 4 5 6 7] 664[ ] 665[4 5 6 7 8] 666[ ] 667[5 6 7 8 9] 668 669 670remove_rows(mat1,{1,3}); 671 672 673[2 3 4 5 6] 674[ ] 675[4 5 6 7 8] 676[ ] 677[5 6 7 8 9] 678 679 680remove_rows(mat1,{1,2,3,4,5}); 681 682***** Warning in remove_rows: 683 all the rows have been removed. Returning nil. 684 685 686swap_columns(mat1,2,4); 687 688 689[1 4 3 2 5] 690[ ] 691[2 5 4 3 6] 692[ ] 693[3 6 5 4 7] 694[ ] 695[4 7 6 5 8] 696[ ] 697[5 8 7 6 9] 698 699 700swap_rows(mat1,1,2); 701 702 703[2 3 4 5 6] 704[ ] 705[1 2 3 4 5] 706[ ] 707[3 4 5 6 7] 708[ ] 709[4 5 6 7 8] 710[ ] 711[5 6 7 8 9] 712 713 714swap_entries(mat1,{1,1},{5,5}); 715 716 717[9 2 3 4 5] 718[ ] 719[2 3 4 5 6] 720[ ] 721[3 4 5 6 7] 722[ ] 723[4 5 6 7 8] 724[ ] 725[5 6 7 8 1] 726 727 728 729 730% Constructors - functions that create matrices. 731 732band_matrix(x,5); 733 734 735[x 0 0 0 0] 736[ ] 737[0 x 0 0 0] 738[ ] 739[0 0 x 0 0] 740[ ] 741[0 0 0 x 0] 742[ ] 743[0 0 0 0 x] 744 745 746band_matrix({x,y,z},6); 747 748 749[y z 0 0 0 0] 750[ ] 751[x y z 0 0 0] 752[ ] 753[0 x y z 0 0] 754[ ] 755[0 0 x y z 0] 756[ ] 757[0 0 0 x y z] 758[ ] 759[0 0 0 0 x y] 760 761 762 763block_matrix(1,2,{mat1,mat2}); 764 765***** Error in block_matrix: row dimensions of 766***** matrices into block_matrix are not compatible 767 768block_matrix(2,3,{mat2,mat3,mat2,mat3,mat2,mat2}); 769 770 771[1 1 1 1 x 1 1 1 1] 772[ ] 773[2 2 2 2 x 2 2 2 2] 774[ ] 775[3 3 3 3 x 3 3 3 3] 776[ ] 777[4 4 4 4 x 4 4 4 4] 778[ ] 779[x 1 1 1 1 1 1 1 1] 780[ ] 781[x 2 2 2 2 2 2 2 2] 782[ ] 783[x 3 3 3 3 3 3 3 3] 784[ ] 785[x 4 4 4 4 4 4 4 4] 786 787 788 789char_matrix(mat1,x); 790 791 792[x - 1 -2 -3 -4 -5 ] 793[ ] 794[ -2 x - 3 -4 -5 -6 ] 795[ ] 796[ -3 -4 x - 5 -6 -7 ] 797[ ] 798[ -4 -5 -6 x - 7 -8 ] 799[ ] 800[ -5 -6 -7 -8 x - 9] 801 802 803 804cfmat := coeff_matrix({x+y+4*z=10,y+x-z=20,x+y+4}); 805 806 807cfmat := { 808 809 [4 1 1] 810 [ ] 811 [-1 1 1] 812 [ ] 813 [0 1 1] 814 815 , 816 817 [z] 818 [ ] 819 [y] 820 [ ] 821 [x] 822 823 , 824 825 [10] 826 [ ] 827 [20] 828 [ ] 829 [-4] 830 831 } 832 833first cfmat * second cfmat; 834 835 836[x + y + 4*z] 837[ ] 838[ x + y - z ] 839[ ] 840[ x + y ] 841 842 843third cfmat; 844 845 846[10] 847[ ] 848[20] 849[ ] 850[-4] 851 852 853 854companion(poly,x); 855 856 857[0 0 0 0 0 0 -12] 858[ ] 859[1 0 0 0 0 0 0 ] 860[ ] 861[0 1 0 0 0 0 0 ] 862[ ] 863[0 0 1 0 0 0 -5 ] 864[ ] 865[0 0 0 1 0 0 -4 ] 866[ ] 867[0 0 0 0 1 0 -1 ] 868[ ] 869[0 0 0 0 0 1 0 ] 870 871 872 873hessian(poly1,{w,x,y,z}); 874 875 876[0 0 0 0 ] 877[ ] 878[ 2 3 2 ] 879[0 2 3*y + z + 1 3*y*z ] 880[ ] 881[ 2 3 2 ] 882[0 3*y + z + 1 6*x*y 3*x*z ] 883[ ] 884[ 2 2 ] 885[0 3*y*z 3*x*z 6*x*y*z] 886 887 888 889hilbert(4,1); 890 891 892[ 1 1 1 ] 893[ 1 --- --- ---] 894[ 2 3 4 ] 895[ ] 896[ 1 1 1 1 ] 897[--- --- --- ---] 898[ 2 3 4 5 ] 899[ ] 900[ 1 1 1 1 ] 901[--- --- --- ---] 902[ 3 4 5 6 ] 903[ ] 904[ 1 1 1 1 ] 905[--- --- --- ---] 906[ 4 5 6 7 ] 907 908 909hilbert(3,y+x); 910 911 912[ - 1 - 1 - 1 ] 913[----------- ----------- -----------] 914[ x + y - 2 x + y - 3 x + y - 4 ] 915[ ] 916[ - 1 - 1 - 1 ] 917[----------- ----------- -----------] 918[ x + y - 3 x + y - 4 x + y - 5 ] 919[ ] 920[ - 1 - 1 - 1 ] 921[----------- ----------- -----------] 922[ x + y - 4 x + y - 5 x + y - 6 ] 923 924 925 926% NOTE WELL. The function tested here used to be called just "jacobian" 927% however us of that name was in conflict with another Reduce package so 928% now it is called mat_jacobian. 929mat_jacobian({x^4,x*y^2,x*y*z^3},{w,x,y,z}); 930 931 932[ 3 ] 933[0 4*x 0 0 ] 934[ ] 935[ 2 ] 936[0 y 2*x*y 0 ] 937[ ] 938[ 3 3 2] 939[0 y*z x*z 3*x*y*z ] 940 941 942 943jordan_block(x,5); 944 945 946[x 1 0 0 0] 947[ ] 948[0 x 1 0 0] 949[ ] 950[0 0 x 1 0] 951[ ] 952[0 0 0 x 1] 953[ ] 954[0 0 0 0 x] 955 956 957 958make_identity(11); 959 960 961[1 0 0 0 0 0 0 0 0 0 0] 962[ ] 963[0 1 0 0 0 0 0 0 0 0 0] 964[ ] 965[0 0 1 0 0 0 0 0 0 0 0] 966[ ] 967[0 0 0 1 0 0 0 0 0 0 0] 968[ ] 969[0 0 0 0 1 0 0 0 0 0 0] 970[ ] 971[0 0 0 0 0 1 0 0 0 0 0] 972[ ] 973[0 0 0 0 0 0 1 0 0 0 0] 974[ ] 975[0 0 0 0 0 0 0 1 0 0 0] 976[ ] 977[0 0 0 0 0 0 0 0 1 0 0] 978[ ] 979[0 0 0 0 0 0 0 0 0 1 0] 980[ ] 981[0 0 0 0 0 0 0 0 0 0 1] 982 983 984 985on rounded; 986 987 % makes things a bit easier to read. 988random_matrix(3,3,100); 989 990 991[ - 8.11911717343 - 75.7167729277 30.62058083 ] 992[ ] 993[ - 50.0325962624 47.1655452861 35.8674263384 ] 994[ ] 995[ - 49.3715543826 - 97.5563670864 - 18.8861862756] 996 997 998on not_negative; 999 1000 1001random_matrix(3,3,100); 1002 1003 1004[43.8999853223 33.7140980286 33.75065406 ] 1005[ ] 1006[49.7333355117 98.9642944905 58.5331568816] 1007[ ] 1008[39.9146060895 67.7954727837 24.8684367642] 1009 1010 1011on only_integer; 1012 1013 1014random_matrix(3,3,100); 1015 1016 1017[16 77 49] 1018[ ] 1019[28 84 51] 1020[ ] 1021[84 56 57] 1022 1023 1024on symmetric; 1025 1026 1027random_matrix(3,3,100); 1028 1029 1030[89 74 91] 1031[ ] 1032[74 95 41] 1033[ ] 1034[91 41 87] 1035 1036 1037off symmetric; 1038 1039 1040on upper_matrix; 1041 1042 1043random_matrix(3,3,100); 1044 1045 1046[41 3 8 ] 1047[ ] 1048[0 31 80] 1049[ ] 1050[0 0 12] 1051 1052 1053off upper_matrix; 1054 1055 1056on lower_matrix; 1057 1058 1059random_matrix(3,3,100); 1060 1061 1062[69 0 0 ] 1063[ ] 1064[34 87 0 ] 1065[ ] 1066[78 72 13] 1067 1068 1069off lower_matrix; 1070 1071 1072on imaginary; 1073 1074 1075off not_negative; 1076 1077 1078random_matrix(3,3,100); 1079 1080 1081[ - 95*i - 72 - 57*i + 59 52*i + 46] 1082[ ] 1083[ - 40*i - 54 70*i 39*i + 28] 1084[ ] 1085[ - 40*i + 45 28*i - 81 9*i + 74 ] 1086 1087 1088off rounded; 1089 1090 1091 1092% toeplitz and vandermonde can take both a list of arguments or just 1093% the arguments. 1094toeplitz({1,2,3,4,5}); 1095 1096 1097[1 2 3 4 5] 1098[ ] 1099[2 1 2 3 4] 1100[ ] 1101[3 2 1 2 3] 1102[ ] 1103[4 3 2 1 2] 1104[ ] 1105[5 4 3 2 1] 1106 1107 1108toeplitz(x,y,z); 1109 1110 1111[x y z] 1112[ ] 1113[y x y] 1114[ ] 1115[z y x] 1116 1117 1118 1119vandermonde({1,2,3,4,5}); 1120 1121 1122[1 1 1 1 1 ] 1123[ ] 1124[1 2 4 8 16 ] 1125[ ] 1126[1 3 9 27 81 ] 1127[ ] 1128[1 4 16 64 256] 1129[ ] 1130[1 5 25 125 625] 1131 1132 1133vandermonde(x,y,z); 1134 1135 1136[ 2] 1137[1 x x ] 1138[ ] 1139[ 2] 1140[1 y y ] 1141[ ] 1142[ 2] 1143[1 z z ] 1144 1145 1146 1147% kronecker_product 1148 1149a1 := mat((1,2),(3,4),(5,6)); 1150 1151 1152 [1 2] 1153 [ ] 1154a1 := [3 4] 1155 [ ] 1156 [5 6] 1157 1158 1159a2 := mat((1,x,1),(2,2,2),(3,3,3)); 1160 1161 1162 [1 x 1] 1163 [ ] 1164a2 := [2 2 2] 1165 [ ] 1166 [3 3 3] 1167 1168 1169 1170kronecker_product(a1,a2); 1171 1172 1173[1 x 1 2 2*x 2 ] 1174[ ] 1175[2 2 2 4 4 4 ] 1176[ ] 1177[3 3 3 6 6 6 ] 1178[ ] 1179[3 3*x 3 4 4*x 4 ] 1180[ ] 1181[6 6 6 8 8 8 ] 1182[ ] 1183[9 9 9 12 12 12] 1184[ ] 1185[5 5*x 5 6 6*x 6 ] 1186[ ] 1187[10 10 10 12 12 12] 1188[ ] 1189[15 15 15 18 18 18] 1190 1191 1192 1193clear a1,a2; 1194 1195 1196 1197% High level algorithms. 1198 1199on rounded; 1200 1201 % makes output easier to read. 1202ch := cholesky(mat7); 1203 1204 1205ch := { 1206 1207 [1 0 0 ] 1208 [ ] 1209 [1 1.41421356237 0 ] 1210 [ ] 1211 [0 0.707106781187 0.707106781187] 1212 1213 , 1214 1215 1216 [1 1 0 ] 1217 [ ] 1218 [0 1.41421356237 0.707106781187] 1219 [ ] 1220 [0 0 0.707106781187] 1221 1222 } 1223 1224tp first ch - second ch; 1225 1226 1227[0 0 0] 1228[ ] 1229[0 0 0] 1230[ ] 1231[0 0 0] 1232 1233 1234tmp := first ch * second ch; 1235 1236 1237 [1 1 0] 1238 [ ] 1239tmp := [1 3.0 1] 1240 [ ] 1241 [0 1 1] 1242 1243 1244tmp - mat7; 1245 1246 1247[0 0 0] 1248[ ] 1249[0 0 0] 1250[ ] 1251[0 0 0] 1252 1253 1254off rounded; 1255 1256 1257 1258gram_schmidt({1,0,0},{1,1,0},{1,1,1}); 1259 1260 1261{{1,0,0},{0,1,0},{0,0,1}} 1262 1263gram_schmidt({1,2},{3,4}); 1264 1265 1266 1 2 2*sqrt(5) - sqrt(5) 1267{{---------,---------},{-----------,------------}} 1268 sqrt(5) sqrt(5) 5 5 1269 1270 1271on rounded; 1272 1273 % again, makes large quotients a bit more readable. 1274% The algorithm used for lu_decom sometimes swaps the rows of the input 1275% matrix so that (given matrix A, lu_decom(A) = {L,U,vec}), we find L*U 1276% does not equal A but a row equivalent of it. The call convert(A,vec) 1277% will return this row equivalent (ie: L*U = convert(A,vec)). 1278lu := lu_decom(mat5); 1279 1280 1281lu := { 1282 1283 [4 0 0 0 ] 1284 [ ] 1285 [1 0.75 0 0 ] 1286 [ ] 1287 [1 0.75 2.0 0 ] 1288 [ ] 1289 [3 0.25 4.0 4.33333333333] 1290 1291 , 1292 1293 1294 [1 1.25 0.25 0.5 ] 1295 [ ] 1296 [0 1 1 0.666666666667] 1297 [ ] 1298 [0 0 1 0 ] 1299 [ ] 1300 [0 0 0 1 ] 1301 1302 , 1303 1304 [3,3,3,4]} 1305 1306mat5; 1307 1308 1309[1 2 1 1] 1310[ ] 1311[1 2 3 1] 1312[ ] 1313[4 5 1 2] 1314[ ] 1315[3 4 5 6] 1316 1317 1318tmp := first lu * second lu; 1319 1320 1321 [4 5.0 1 2.0] 1322 [ ] 1323 [1 2.0 1 1 ] 1324tmp := [ ] 1325 [1 2.0 3.0 1 ] 1326 [ ] 1327 [3 4.0 5.0 6.0] 1328 1329 1330tmp1 := convert(mat5,third lu); 1331 1332 1333 [4 5 1 2] 1334 [ ] 1335 [1 2 1 1] 1336tmp1 := [ ] 1337 [1 2 3 1] 1338 [ ] 1339 [3 4 5 6] 1340 1341 1342tmp - tmp1; 1343 1344 1345[0 0 0 0] 1346[ ] 1347[0 0 0 0] 1348[ ] 1349[0 0 0 0] 1350[ ] 1351[0 0 0 0] 1352 1353 1354% and the complex case... 1355lu1 := lu_decom(mat6); 1356 1357 1358lu1 := { 1359 1360 [ 1 0 0 ] 1361 [ ] 1362 [ 4 - 4*i + 5 0 ] 1363 [ ] 1364 [i + 1 3 0.414634146341*i + 2.26829268293] 1365 1366 , 1367 1368 1369 [1 i 0 ] 1370 [ ] 1371 [0 1 0.19512195122*i + 0.243902439024] 1372 [ ] 1373 [0 0 1 ] 1374 1375 , 1376 1377 [3,2,3]} 1378 1379mat6; 1380 1381 1382[i + 1 i + 2 i + 3] 1383[ ] 1384[ 4 5 2 ] 1385[ ] 1386[ 1 i 0 ] 1387 1388 1389tmp := first lu1 * second lu1; 1390 1391 1392 [ 1 i 0 ] 1393 [ ] 1394tmp := [ 4 5 2.0 ] 1395 [ ] 1396 [i + 1 i + 2 i + 3.0] 1397 1398 1399tmp1 := convert(mat6,third lu1); 1400 1401 1402 [ 1 i 0 ] 1403 [ ] 1404tmp1 := [ 4 5 2 ] 1405 [ ] 1406 [i + 1 i + 2 i + 3] 1407 1408 1409tmp - tmp1; 1410 1411 1412[0 0 0] 1413[ ] 1414[0 0 0] 1415[ ] 1416[0 0 0] 1417 1418 1419 1420mat9inv := pseudo_inverse(mat9); 1421 1422 1423 [ - 0.199999999996 0.100000000013 ] 1424 [ ] 1425 [ - 0.0499999999988 0.0500000000037 ] 1426mat9inv := [ ] 1427 [ 0.0999999999982 - 5.57818374824e-12] 1428 [ ] 1429 [ 0.249999999995 - 0.0500000000148 ] 1430 1431 1432mat9 * mat9inv; 1433 1434 1435[ 0.999999999982 - 0.0000000000557817542157] 1436[ ] 1437[5.54267742814e-12 1.00000000002 ] 1438 1439 1440 1441simplex(min,2*x1+14*x2+36*x3,{-2*x1+x2+4*x3>=5,-x1-2*x2-3*x3<=2}); 1442 1443 1444{45.0,{x1 = 0,x2 = 0,x3 = 1.25}} 1445 1446 1447simplex(max,10000 x1 + 1000 x2 + 100 x3 + 10 x4 + x5,{ x1 <= 1, 20 x1 + 1448 x2 <= 100, 200 x1 + 20 x2 + x3 <= 10000, 2000 x1 + 200 x2 + 20 x3 + x4 1449 <= 1000000, 20000 x1 + 2000 x2 + 200 x3 + 20 x4 + x5 <= 100000000}); 1450 1451 1452{100000000,{x1 = 0,x2 = 0,x3 = 0,x4 = 0,x5 = 100000000}} 1453 1454 1455simplex(max, 5 x1 + 4 x2 + 3 x3, 1456 { 2 x1 + 3 x2 + x3 <= 5, 1457 4 x1 + x2 + 2 x3 <= 11, 1458 3 x1 + 4 x2 + 2 x3 <= 8 }); 1459 1460 1461{13.0,{x1 = 2.0,x2 = 0,x3 = 1}} 1462 1463 1464simplex(min,3 x1 + 5 x2,{ x1 + 2 x2 >= 2, 22 x1 + x2 >= 3}); 1465 1466 1467{5.04651162791,{x1 = 0.093023255814,x2 = 0.953488372093}} 1468 1469 1470simplex(max,10x+5y+5.5z,{5x+3z<=200,0.2x+0.1y+0.5z<=12,0.1x+0.2y+0.3z<=9, 1471 30x+10y+50z<=1500}); 1472 1473 1474{525.0,{x = 40.0,y = 25.0,z = 0}} 1475 1476 1477%example of extra variables (>=0) being added. 1478simplex(min,x-y,{x>=-3}); 1479 1480*** Warning: variable y not defined in input. Has been defined as >=0. 1481 1482***** Error in simplex: The problem is unbounded. 1483 1484 1485% unfeasible as simplex algorithm implies all x>=0. 1486simplex(min,x,{x<=-100}); 1487 1488 1489***** Error in simplex: Problem has no feasible solution. 1490 1491 1492% three error examples. 1493simplex(maxx,x,{x>=5}); 1494 1495 1496***** Error in simplex(first argument): must be either max or min 1497 1498simplex(max,x,x>=5); 1499 1500 1501***** Error in simplex(third argument): must be a list 1502 1503simplex(max,x,{x<=y}); 1504 1505 1506***** Error in simplex: Right hand side of inequalities must be numbers. 1507 1508 1509simplex(max, 346 X11 + 346 X12 + 248 X21 + 248 X22 + 399 X31 + 399 X32 + 1510 200 Y11 + 200 Y12 + 75 Y21 + 75 Y22 + 2.35 Z1 + 3.5 Z2, 1511{ 1512 4 X11 + 4 X12 + 2 X21 + 2 X22 + X31 + X32 + 250 Y11 + 250 Y12 + 125 Y21 + 1513 125 Y22 <= 25000, 1514 X11 + X12 + X21 + X22 + X31 + X32 + 2 Y11 + 2 Y12 + Y21 + Y22 <= 300, 1515 20 X11 + 15 X12 + 30 Y11 + 20 Y21 + Z1 <= 1500, 1516 40 X12 + 35 X22 + 50 X32 + 15 Y12 + 10 Y22 + Z2 = 5000, 1517 X31 = 0, 1518 Y11 + Y12 <= 50, 1519 Y21 + Y22 <= 100 1520}); 1521 1522 1523{99250.0, 1524 1525 {x11 = 75.0, 1526 1527 x12 = 0, 1528 1529 x21 = 225.0, 1530 1531 x22 = 0, 1532 1533 x31 = 0, 1534 1535 x32 = 0, 1536 1537 y11 = 0, 1538 1539 y12 = 0, 1540 1541 y21 = 0, 1542 1543 y22 = 0, 1544 1545 z1 = 0, 1546 1547 z2 = 5000.0}} 1548 1549 1550 1551% from Marc van Dongen. Finding the first feasible solution for the 1552% solution of systems of linear diophantine inequalities. 1553simplex(max,0,{ 1554 3*X259+4*X261+3*X262+2*X263+X269+2*X270+3*X271+4*X272+5*X273+X229=2, 1555 7*X259+11*X261+8*X262+5*X263+3*X269+6*X270+9*X271+12*X272+15*X273+X229=4, 1556 2*X259+5*X261+4*X262+3*X263+3*X268+4*X269+5*X270+6*X271+7*X272+8*X273=1, 1557 X262+2*X263+5*X268+4*X269+3*X270+2*X271+X272+2*X229=1, 1558 X259+X262+2*X263+4*X268+3*X269+2*X270+X271-X273+3*X229=2, 1559 X259+2*X261+2*X262+2*X263+3*X268+3*X269+3*X270+3*X271+3*X272+3*X273+X229=1, 1560 X259+X261+X262+X263+X268+X269+X270+X271+X272+X273+X229=1}); 1561 1562 1563{0, 1564 1565 {x229 = 0.5, 1566 1567 x259 = 0.5, 1568 1569 x261 = 0, 1570 1571 x262 = 0, 1572 1573 x263 = 0, 1574 1575 x268 = 0, 1576 1577 x269 = 0, 1578 1579 x270 = 0, 1580 1581 x271 = 0, 1582 1583 x272 = 0, 1584 1585 x273 = 0}} 1586 1587 1588 1589% An example with a bound section: 1590simplex(min, 0, 1591 {-n2 <= -1, -n1-n2 <= 0, 2*n1+n2 <= 0, 5*n1+2*n2 <= 0, 5*n1-n2 <= 0}, 1592 {-infinity <= n1, -infinity <= n2 <= 2}); 1593 1594 1595{0,{n1 = - 0.5,n2 = 1}} 1596 1597 1598 1599svd_ans := svd(mat8); 1600 1601 1602svd_ans := { 1603 1604 [ 0.289784137735 0.957092029805] 1605 [ ] 1606 [ - 0.957092029805 0.289784137735] 1607 1608 , 1609 1610 1611 [5.1491628629 0 ] 1612 [ ] 1613 [ 0 2.9130948854] 1614 1615 , 1616 1617 1618 [ - 0.687215403194 0.726453707825 ] 1619 [ ] 1620 [ - 0.726453707825 - 0.687215403194] 1621 1622 } 1623 1624tmp := tp first svd_ans * second svd_ans * third svd_ans; 1625 1626 1627 [ 0.99999998509 2.9999999859 ] 1628tmp := [ ] 1629 [ - 4.00000004924 2.99999995342] 1630 1631 1632tmp - mat8; 1633 1634 1635[ - 0.0000000149096008872 - 0.0000000141042799662] 1636[ ] 1637[ - 0.000000049243064737 - 0.000000046583274127 ] 1638 1639 1640 1641mat9inv := pseudo_inverse(mat9); 1642 1643 1644 [ - 0.199999999996 0.100000000013 ] 1645 [ ] 1646 [ - 0.0499999999988 0.0500000000037 ] 1647mat9inv := [ ] 1648 [ 0.0999999999982 - 5.57818374824e-12] 1649 [ ] 1650 [ 0.249999999995 - 0.0500000000148 ] 1651 1652 1653mat9 * mat9inv; 1654 1655 1656[ 0.999999999982 - 0.0000000000557817542157] 1657[ ] 1658[5.54267742814e-12 1.00000000002 ] 1659 1660 1661 1662% triang_adjoint(in_mat) calculates the 1663% triangularizing adjoint of in_mat 1664 1665triang_adjoint(mat1); 1666 1667 1668[1 0 0 0 0] 1669[ ] 1670[-2 1 0 0 0] 1671[ ] 1672[-1 2 -1 0 0] 1673[ ] 1674[0 0 0 0 0] 1675[ ] 1676[0 0 0 0 0] 1677 1678 1679triang_adjoint(mat2); 1680 1681 1682[1 0 0 0] 1683[ ] 1684[-2 1 0 0] 1685[ ] 1686[0 0 0 0] 1687[ ] 1688[0 0 0 0] 1689 1690 1691triang_adjoint(mat5); 1692 1693 1694[1 0 0 0] 1695[ ] 1696[-1 1 0 0] 1697[ ] 1698[-3 3 0 0] 1699[ ] 1700[10 -12 -4 6] 1701 1702 1703triang_adjoint(mat6); 1704 1705 1706[ 1 0 0 ] 1707[ ] 1708[ -4 i + 1 0 ] 1709[ ] 1710[4*i - 5 3 i - 3] 1711 1712 1713triang_adjoint(mat7); 1714 1715 1716[1 0 0] 1717[ ] 1718[-1 1 0] 1719[ ] 1720[1 - 1 2] 1721 1722 1723triang_adjoint(mat8); 1724 1725 1726[1 0] 1727[ ] 1728[4 1] 1729 1730 1731triang_adjoint(mat9); 1732 1733 1734***** Error in triang_adjoint: input matrix should be square. 1735 1736 1737% testing triang_adjoint with random matrices 1738 1739% the range of the integers is in one case from 1740% -1000 to 1000. in the other case it is from 1741% -1 to 1 so that the deteminant of the i-th 1742% submatrix equals very often to zero. 1743 1744% random matrix contains arbitrary real values 1745off only_integer; 1746 1747 1748tmp:=random_matrix(5,5,1000); 1749 1750 1751tmp := mat(( - 558.996086656*i + 1.71931812953,76.9987188735*i + 1.19004104683, 1752 1753 - 739.283479439*i - 1.32534106204,742.101952123*i + 1.35926854848, 1754 1755 680.515777254*i + 1.56403177895), 1756 1757 ( - 689.196170962*i + 1.49289170118, 1758 1759 - 232.584493916*i - 1.38227180395,280.109305836*i + 1.38865247861, 1760 1761 298.151479065*i - 1.19035182389, - 602.312143386*i - 1.82183796879), 1762 1763 ( - 739.195910955*i - 1.45944960213,859.293884826*i + 1.70488070379, 1764 1765 359.856032683*i - 1.2966991869, - 105.409833087*i - 1.9360631701, 1766 1767 234.350529301*i - 1.15598520849), 1768 1769 (155.629059348*i + 1.09264385739, - 16.1559469072*i - 1.9425176505, 1770 1771 725.11578405*i - 1.05760723025,783.020942195*i - 1.28625265346, 1772 1773 - 544.129360355*i + 1.74790906085), 1774 1775 (373.562370318*i - 1.95218354686, - 722.109349973*i + 1.56309793677, 1776 1777 - 746.664959169*i - 1.9915755693,186.154794517*i - 1.09842189916, 1778 1779 435.90998986*i - 1.46175649496)) 1780 1781 1782triang_adjoint tmp; 1783 1784 1785mat((1,0,0,0,0), 1786 1787 (689.196170962*i - 1.49289170118, - 558.996086656*i + 1.71931812953,0,0,0), 1788 1789 ( - 1253.37955588*i + 7.64148089854e+5, - 1516.42713845*i - 4.23429448803e+5 1790 1791 ,1078.01877642*i - 1.830851973e+5,0,0), 1792 1793 102791325687.0*i + 7.3752778526e+8 1794 (------------------------------------, 1795 i - 169.834887206 1796 1797 - 3.66748178757e+10*i - 6.62162769101e+6 1798 -------------------------------------------, 1799 i - 169.834887206 1800 1801 9.85957444629e+7*i - 1.01033337718e+6, 1802 1803 - 7.49414742893e+8*i - 2.25311577415e+6,0), 1804 1805 - 547052849318.0*i + 4.06181988045e+13 1806 (-----------------------------------------, 1807 i - 112.974983172 1808 1809 - 141265342333.0*i + 4.13350560163e+12 1810 -----------------------------------------, 1811 i - 112.974983172 1812 1813 845804392649.0*i - 9.62488227345e+13 1814 --------------------------------------, 1815 i - 112.974983172 1816 1817 876106032577.0*i - 2.66464796763e+13 1818 --------------------------------------, 1819 i - 112.974983172 1820 1821 1.47617976407e+12*i - 1.66771384004e+14 1822 -----------------------------------------)) 1823 i - 169.834887206 1824 1825 1826 1827tmp:=random_matrix(1,1,1000); 1828 1829 1830tmp := [ - 463.860434427*i + 1.35500571348] 1831 1832 1833triang_adjoint tmp; 1834 1835 1836[1] 1837 1838 1839 1840% random matrix contains complex real values 1841on imaginary; 1842 1843 1844tmp:=random_matrix(5,5,1000); 1845 1846 1847tmp := mat((107.345792105*i - 1.98704739339,188.868545358*i + 1.22417796742, 1848 1849 - 630.485915434*i + 1.32195292724, 1850 1851 - 542.510039297*i - 1.94318764036,359.860945563*i - 1.69174206177), 1852 1853 ( - 469.501213476*i - 1.17375946319, - 62.2197820375*i - 1.4051578261 1854 1855 , - 98.6604380996*i + 1.64691610034, 1856 1857 - 216.296595937*i + 1.56809020199,797.19877393*i - 1.31894550853), 1858 1859 (2.07054207792*i + 1.3891068942,393.038868455*i - 1.60894498437, 1860 1861 - 215.390393738*i - 1.00068640594, 1862 1863 - 195.674928032*i + 1.22123114986,211.921323796*i - 1.42499533273), 1864 1865 ( - 750.357435524*i - 1.19871674827, 1866 1867 - 792.333836712*i - 1.63151974094, - 494.87049225*i + 1.99554801527 1868 1869 ,638.173945822*i + 1.23793954612,111.418959978*i - 1.96273029328), 1870 1871 ( - 255.359922267*i + 1.99035939892, 1872 1873 - 575.376389757*i - 1.03533681609,463.961589382*i - 1.86476410547, 1874 1875 83.8856338571*i + 1.10369785887, - 129.597812786*i - 1.4917934624)) 1876 1877 1878triang_adjoint tmp; 1879 1880 1881mat((1,0,0,0,0), 1882 1883 (469.501213476*i + 1.17375946319,107.345792105*i - 1.98704739339,0,0,0), 1884 1885 (383.407897912*i + 1.84407237435e+5,1218.59364331*i + 41798.5118562, 1886 1887 769.235159465*i - 81990.7504399,0,0), 1888 1889 - 1.411092405e+10*i - 1.91497165215e+8 1890 (-----------------------------------------, 1891 i - 106.587367245 1892 1893 - 2.06157034475e+10*i + 1.09218575639e+8 1894 -------------------------------------------, 1895 i - 106.587367245 1896 1897 - 2.4008888901e+8*i + 13175.2571592, 1898 1899 - 1.02728261373e+8*i + 9.22309484944e+5,0), 1900 1901 - 203213290519.0*i - 3.07405185302e+12 1902 (-----------------------------------------, 1903 i - 184.764270765 1904 1905 311149245317.0*i + 2.05618234856e+13 1906 --------------------------------------, 1907 i - 184.764270765 1908 1909 212889617996.0*i - 4.13210409411e+13 1910 --------------------------------------, 1911 i - 184.764270765 1912 1913 - 7.79955805661e+10*i - 5.10418442965e+12 1914 --------------------------------------------, 1915 i - 184.764270765 1916 1917 7.62835257557e+10*i - 1.40944700076e+13 1918 -----------------------------------------)) 1919 i - 106.587367245 1920 1921 1922 1923tmp:=random_matrix(1,1,1000); 1924 1925 1926tmp := [276.278111177*i + 1.74724262616] 1927 1928 1929triang_adjoint tmp; 1930 1931 1932[1] 1933 1934 1935off imaginary; 1936 1937 1938 1939% random matrix contains rounded real values 1940on rounded; 1941 1942 1943tmp:=random_matrix(5,5,1000); 1944 1945 1946tmp := mat(( - 293.224093687, - 99.023221037, - 819.400851656,796.020234589, 1947 1948 593.862087611), 1949 1950 ( - 137.84203019,354.3234619, - 852.314261681, - 217.485901759, 1951 1952 256.139775139), 1953 1954 (324.37828726, - 56.5718498235, - 118.293003834,108.279501424, 1955 1956 23.2385400299), 1957 1958 ( - 976.556156754,684.207160793,146.328625386,502.848132905, 1959 1960 312.766816689), 1961 1962 (211.783458501,166.556239469,175.715904944,251.57997022,280.123720131 1963 1964 )) 1965 1966 1967triang_adjoint tmp; 1968 1969 1970mat((1,0,0,0,0), 1971 1972 (137.84203019, - 293.224093687,0,0,0), 1973 1974 ( - 1.07136859076e+5, - 48709.2122316, - 1.17545737812e+5,0,0), 1975 1976 (1.27980020917e+8, - 1.64635380167e+8,4.76863677307e+8,1.43208428244e+8,0), 1977 1978 (5.82963241185e+10,3.9383738062e+10, - 437637051137.0, - 111757830528.0, 1979 1980 261327212376.0)) 1981 1982 1983 1984tmp:=random_matrix(1,1,1000); 1985 1986 1987tmp := [406.584701921] 1988 1989 1990triang_adjoint tmp; 1991 1992 1993[1] 1994 1995 1996off rounded; 1997 1998 1999 2000% random matrix contains only integer values 2001on only_integer; 2002 2003 2004tmp:=random_matrix(7,7,1000); 2005 2006 2007 [969 210 8 244 -887 -39 -916] 2008 [ ] 2009 [-774 296 -475 -694 -909 560 89 ] 2010 [ ] 2011 [-390 -559 -551 -567 241 -306 -655] 2012 [ ] 2013tmp := [-478 809 181 -987 -144 929 -886] 2014 [ ] 2015 [188 267 -778 660 374 590 30 ] 2016 [ ] 2017 [ 73 971 -946 -43 -215 386 -365] 2018 [ ] 2019 [-792 -852 558 -797 343 219 110 ] 2020 2021 2022triang_adjoint tmp; 2023 2024 2025mat((1,0,0,0,0,0,0), 2026 2027 (774,969,0,0,0,0,0), 2028 2029 (548106,459771,449364,0,0,0,0), 2030 2031 (-108937808,399369604,-497500435,-461605941,0,0,0), 2032 2033 (-386678984240,-1001551613816,454549593485,637690866447,433944480084,0,0), 2034 2035 (-604165739229705,-320961967400919,-165015285307395,-1008712187269380, 2036 2037 -1670995725485274,1433408878792557,0), 2038 2039 (-181830640557070260,295390292387079435,816541226477288004, 2040 2041 850494616785589377,458176543109779557,-1709784109660828152, 2042 2043 -1475366833406131953)) 2044 2045 2046 2047tmp:=random_matrix(7,7,1); 2048 2049 2050 [0 0 0 0 0 0 0] 2051 [ ] 2052 [0 0 0 0 0 0 0] 2053 [ ] 2054 [0 0 0 0 0 0 0] 2055 [ ] 2056tmp := [0 0 0 0 0 0 0] 2057 [ ] 2058 [0 0 0 0 0 0 0] 2059 [ ] 2060 [0 0 0 0 0 0 0] 2061 [ ] 2062 [0 0 0 0 0 0 0] 2063 2064 2065triang_adjoint tmp; 2066 2067 2068[1 0 0 0 0 0 0] 2069[ ] 2070[0 0 0 0 0 0 0] 2071[ ] 2072[0 0 0 0 0 0 0] 2073[ ] 2074[0 0 0 0 0 0 0] 2075[ ] 2076[0 0 0 0 0 0 0] 2077[ ] 2078[0 0 0 0 0 0 0] 2079[ ] 2080[0 0 0 0 0 0 0] 2081 2082 2083 2084% random matrix contains only complex integer 2085% values 2086 2087on imaginary; 2088 2089 2090tmp:=random_matrix(5,5,1000); 2091 2092 2093tmp := mat((12*(38*i + 83),3*(153*i - 305),2*(141*i + 427), - 553*i + 617, 2094 2095 3*(83*i + 157)), 2096 2097 (164*i - 635, - 133*i + 991, - 373*i + 210,965*i - 608,2*(358*i - 55) 2098 2099 ), 2100 2101 ( - 230*i + 227,32*i + 339,2*(485*i - 219),707*i - 767, - 985*i - 51) 2102 2103 , 2104 2105 ( - 609*i + 647, - 441*i + 187,930*i - 349,250*i - 211,274*i - 451), 2106 2107 ( - 374*i - 135,793*i + 592, - 81*i - 1,89*i + 26,3*( - 40*i + 201))) 2108 2109 2110triang_adjoint tmp; 2111 2112 2113mat((1,0,0,0,0), 2114 2115 ( - 164*i + 635,12*(38*i + 83),0,0,0), 2116 2117 (293397*i - 414880,9*(14243*i - 47243),3*(253651*i + 180645),0,0), 2118 2119 - 253324472288717*i + 71265413812547 2120 (---------------------------------------, 2121 253651*i + 180645 2122 2123 2*( - 220885726602145*i - 1441709355714) 2124 ------------------------------------------, - 1436348339*i + 1393250309, 2125 253651*i + 180645 2126 2127 511458435*i - 1454012933,0), 2128 2129 13983048003979950612955437881*i - 71498490838832832842693585028 2130 (-----------------------------------------------------------------, 2131 65634686423804933*i - 9174596297286164 2132 2133 89295323223054915316808489269*i - 37624299403809895760446255007 2134 -----------------------------------------------------------------, 2135 65634686423804933*i - 9174596297286164 2136 2137 2*( - 71881165390656818494884812727*i - 25318671134083617432051412624) 2138 ------------------------------------------------------------------------, 2139 65634686423804933*i - 9174596297286164 2140 2141 134577377248105484011524135103*i + 3495516738012600790097438251 2142 -----------------------------------------------------------------, 2143 65634686423804933*i - 9174596297286164 2144 2145 6*(65634686423804933*i - 9174596297286164) 2146 --------------------------------------------)) 2147 253651*i + 180645 2148 2149 2150 2151tmp:=random_matrix(5,5,2); 2152 2153 2154 [i - 1 i i 0 - (i + 1)] 2155 [ ] 2156 [ 0 i -1 - i + 1 i + 1 ] 2157 [ ] 2158tmp := [ -1 0 0 - i + 1 -1 ] 2159 [ ] 2160 [ -1 - i - i - i i + 1 ] 2161 [ ] 2162 [i - 1 0 i + 1 -1 0 ] 2163 2164 2165triang_adjoint tmp; 2166 2167 2168[ 1 0 0 0 0 ] 2169[ ] 2170[ 0 i - 1 0 0 0 ] 2171[ ] 2172[ - (i + 1) i + 1 ] 2173[------------ ------- - (i + 1) 0 0 ] 2174[ i - 1 i - 1 ] 2175[ ] 2176[ - (i + 1) 2*(2*i + 1) - 2*i ] 2177[------------ 0 ------------- -------- 0 ] 2178[ i i - 1 i - 1 ] 2179[ ] 2180[ 2*(3*i - 4) 2*(i + 2) 5*(3*i + 1) - 7*i + 1 2*(i + 2) ] 2181[------------- ----------- ------------- ------------ -----------] 2182[ 4*i + 3 i - 1 4*i + 3 4*i + 3 i - 1 ] 2183 2184 2185 2186% Predicates. 2187 2188matrixp(mat1); 2189 2190 2191t 2192 2193matrixp(poly); 2194 2195 2196 2197squarep(mat2); 2198 2199 2200t 2201 2202squarep(mat3); 2203 2204 2205 2206symmetricp(mat1); 2207 2208 2209t 2210 2211symmetricp(mat3); 2212 2213 2214 2215if not rounded_was_on then off rounded; 2216 2217 2218 2219 2220END; 2221 2222Tested on x86_64-pc-windows CSL 2223Time (counter 1): 31 ms 2224 2225End of Lisp run after 0.03+0.06 seconds 2226real 0.23 2227user 0.04 2228sys 0.04 2229