On 2014-10-25, Jason Grout <jason-s...@creativetrax.com> wrote:

>   http://www.ams.org/notices/201410/rnoti-p1249.pdf

> P.S. It would be interesting to see if Sage can do the calculation they 
> identified as buggy in mathematica.  That would make for a cool 
> follow-up editorial.

I've reimplemnted the buggy determinant calculation in Maxima. 
Presumably from this it would be easy to redo it in any other system;
I don't know how Sage manages such calculations.

I am happy to report that Maxima, despite its many and varied bugs,
doesn't have this particular one:

    bfloat (determinant (big_matrix));
     => 1.951242191319868b9762

Reported value in paper is 1.95124219131987 * 10^9762.

Script is attached as a PS. The function foo(n) can be used to generate
random examples, as the authors did to find one which tickles the bug.

All in good fun,

Robert Dodier

PS.
basic_matrix : matrix
 ([-32, 69, 89, -60, -83, -22, -14, -58, 85, 56, -65, -30, -86, -9],
  [6, 99, 11, 57, 47, -42, -48, -65, 25, 50, -70, -3, -90, 31],
  [78, 38, 12, 64, -67, -4, -52, -65, 19, 71, 38, -17, 51, -3],
  [-93, 30, 89, 22, 13, 48, -73, 93, 11, -97, -49, 61, -25, -4],
  [54, -22, 54, -53, -52, 64, 19, 1, 81, -72, -11, 50, 0, -81],
  [65, -58, 3, 57, 19, 77, 76, -57, -80, 22, 93, -85, 67, 58],
  [29, -58, 47, 87, 3, -6, -81, 5, 98, 86, -98, 51, -62, -66],
  [93, -77, 16, -64, 48, 84, 97, 75, 89, 63, 34, -98, -94, 19],
  [45, -99, 3, -57, 32, 60, 74, 4, 69, 98, -40, -69, -28, -26],
  [-13, 51, -99, -2, 48, 71, -81, -32, 78, 27, -28, -22, 22, 94],
  [11, 72, -74, 86, 79, -58, -89, 80, 70, 55, -49, 51, -42, 66],
  [-72, 53, 49, -46, 17, -22, -48, -40, -28, -85, 88, -30, 74, 32],
  [-92, -22, -90, 67, -25, -28, -91, -8, 32, -41, 10, 6, 85, 21],
  [47, -73, -30, -60, 99, 9, -86, -70, 84, 55, 19, 69, 11, -84]) $

small_matrix : matrix
 ([528, 853, -547, -323, 393, -916, -11, -976, 279, -665, 906, -277, 103, -485],
  [878, 910, -306, -260, 575, -765, -32, 94, 254, 276, -156, 625, -8, -566],
  [-357, 451, -475, 327, -84, 237, 647, 505, -137, 363, -808, 332, 222, -998],
  [-76, 26, -778, 505, 942, -561, -350, 698, -532, -507, -78, -758, 346, -545],
  [-358, 18, -229, -880, -955, -346, 550, -958, 867, -541, -962, 646, 932, 168],
  [192, 233, 620, 955, -877, 281, 357, -226, -820, 513, -882, 536, -237, 877],
  [-234, -71, -831, 880, -135, -249, -427, 737, 664, 298, -552, -1, -712, -691],
  [80, 748, 684, 332, 730, -111, -643, 102, -242, -82, -28, 585, 207, -986],
  [967, 1, -494, 633, 891, -907, -586, 129, 688, 150, -501, -298, 704, -68],
  [406, -944, -533, -827, 615, 907, -443, -350, 700, -878, 706, 1, 800, 120],
  [33, -328, -543, 583, -443, -635, 904, -745, -398, -110, 751, 660, 474, 255],
  [-537, -311, 829, 28, 175, 182, -930, 258, -808, -399, -43, -68, -553, 421],
  [-373, -447, -252, -619, -418, 764, 994, -543, -37, -845, 30, -704, 147, 
-534],
  [638, -33, 932, -335, -75, -676, -934, 239, 210, 665, 414, -803, 564, -805]);

powers_list : [10^123,
  10^152, 10^185, 10^220, 10^397, 10^449,
  10^503, 10^563, 10^979, 10^1059, 10^1143,
  10^1229, 10^1319, 10^1412];

powers_matrix : apply (diag_matrix, powers_list);

big_matrix : basic_matrix . powers_matrix + small_matrix;

bfloat (determinant (big_matrix));

foo (n) := block ([basic_matrix, powers_matrix, small_matrix, big_matrix, a, b],
  basic_matrix : genmatrix (lambda ([i, j], random (2 * 99 + 1) - 99), n, n),
  powers_matrix : apply (diag_matrix, rest (powers_list, n - length 
(powers_list))),
  small_matrix : genmatrix (lambda ([i, j], random (2 * 999 + 1) - 999), n, n),
  big_matrix : basic_matrix . powers_matrix + small_matrix,
  a : determinant (big_matrix),
  b : determinant (big_matrix),
  is (a = b)) $

set_random_state (make_random_state (1729)) $

foo (4);

-- 
You received this message because you are subscribed to the Google Groups 
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-devel+unsubscr...@googlegroups.com.
To post to this group, send email to sage-devel@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.

Reply via email to