C Find determinant of n*n integer matrix with elements 1..n^2 C Find target values not expressible by such a determinant C Result is OEIS A088238 C C Hugo Pfoertner http://www.pfoertner.org/ C C Version history: C 03.10.03 Find all possible determinants in given range for n=5. C Generation of random permutation by sorting restored. C 28.09.03 Use Fisher, Yates algo from Knuth TAOCP vol. 2, p145 to C generate random permutations, Cramer's rule for n=5 C 25.09.03 Search for given target value C 23.09.03 Search by exchange steps C 19.09.03 5*5 matrix included C 17.09.03 Initial version C implicit integer (a-z) C C Whole range of representible determinants plus beginning of range C with not representible determinants (15 min CPU on Athlon 800Mhz) C Expected result: 102 not representible values >=6773999 parameter ( n=5, nn=n*n, nlo=0, nhi=6783000, conmax=50 ) C C Parameters for searching upper range (may take a few days CPU) C For people with a lot of CPU resources: Run with conmax=50 or 100 C This will allow to find OEIS A088217(5): C A088217(5)=2*(A085000(5)-l_not_rep(5))+1 C parameter ( n=5, nn=n*n, nlo=6773990, nhi=6839492, conmax=10 ) C Whole range for n=4 c parameter ( n=4, nn=n*n, nlo=0, nhi=40800, conmax=50 ) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. dimension a(nn) integer*1 s(nlo:nhi) real x(nn) integer ip(nn) C Name of uniform (0,1) RNG (replace by available RNG if necessary) real urand C C STATEMENT FUNCTIONS for 2*2, 3*3, 4*4, 5*5 determinants C D2(I1,I2,I3,I4)=I1*I4-I2*I3 C D3(I1,I2,I3,I4,I5,I6,I7,I8,I9)= & I1*D2(I5,I6,I8,I9) & -I2*D2(I4,I6,I7,I9) & +I3*D2(I4,I5,I7,I8) C D4(I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15,I16)= & I1*D3(I6,I7,I8,I10,I11,I12,I14,I15,I16) & -I2*D3(I5,I7,I8,I9,I11,I12,I13,I15,I16) & +I3*D3(I5,I6,I8,I9,I10,I12,I13,I14,I16) & -I4*D3(I5,I6,I7,I9,I10,I11,I13,I14,I15) C D5(I1,I2,I3,I4,I5,I6,I7,I8,I9,I10, & I11,I12,I13,I14,I15,I16,I17,I18,I19,I20, & I21,I22,I23,I24,I25)= & I1*D4(I7,I8,I9,I10,I12,I13,I14,I15, & I17,I18,I19,I20,I22,I23,I24,I25) & -I2*D4(I6,I8,I9,I10,I11,I13,I14,I15, & I16,I18,I19,I20,I21,I23,I24,I25) & +I3*D4(I6,I7,I9,I10,I11,I12,I14,I15, & I16,I17,I19,I20,I21,I22,I24,I25) & -I4*D4(I6,I7,I8,I10,I11,I12,I13,I15, & I16,I17,I18,I20,I21,I22,I23,I25) & +I5*D4(I6,I7,I8,I9,I11,I12,I13,I14, & I16,I17,I18,I19,I21,I22,I23,I24) C C Activate if program shall be used to find matrix for C a given single target C Read target value c write (*,*) ' Target:' c read (*,*) target c target = abs(target) C C Get starting time of computation from system and use C as initial value for seed generation C secnds () ... Time since beginning of day or year. IRAN = secnds () C IF ( MOD ( IRAN, 2 ) .EQ. 0 ) IRAN = IRAN + 1 C Preset array indicating if a certain determinant was found do 290 i = nlo, nhi s(i) = 0 290 continue C Alternatively read previous results from input file c open ( unit=20, file='det5.dat', form='unformatted', c & status='old', iostat=ios ) C if ( iostat .ne. 0 ) stop ' Missing input file.' c read (20) s c close ( unit=20 ) C C Preset counter of target values for which no determinant was C found in previous run over all targets lnotpr = -1 C C Starting point of next run through whole range to be searched C 600 continue C Count values for which no determinant has been found so far lnot = 0 do 610 i = nlo, nhi if ( s(i) .eq. 0 ) lnot = lnot + 1 610 continue C Compare with counts from previous run if ( lnot .ne. lnotpr ) then lnotpr = lnot confir = 0 write (*,*) ' lnot=', lnot else C No new determinants found in last run confir = confir + 1 write (*,*) ' Confirmed ', confir, ' times: ', lnot C Check stopping criterion if ( confir .gt. conmax ) goto 900 endif C Reset counter for unsuccessful searches in current run lnot = 0 C C Loop over whole target range C do 300 target = nlo, nhi C Omit already found determinant values if ( s(target) .gt. 0 ) goto 300 C C Best value found so far (only for single target diagnostic print) C dglob = 0 C C Generate random matrices as starting points for exchange sweeps C do 500 rantry = 1, 1000 C C Generate random permutation vector by Fisher/Yates method as C described in Knuth's TAOCP vol. 2 c do 510 i = 1, nn c a(i) = i c510 continue c do 520 j = nn, 2, -1 c k = 1 + int ( float(j) * urand(iran) ) c i = a(j) c a(j) = a(k) c a(k) = i c520 continue C C Generate random permutation by sorting a random array x do 10 i = 1, nn x(i) = urand(iran) 10 continue C Permutation vector by sorting (SLATEC) call spsort ( x, nn, ip, -1, ier ) C set matrix elements do 20 i = 1, nn a(ip(i)) = i 20 continue C Determinant of new random matrix C d = d4 (a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9), C & a(10),a(11),a(12),a(13),a(14),a(15),a(16)) d = d5 (a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9), & a(10),a(11),a(12),a(13),a(14),a(15),a(16),a(17), & a(18),a(19),a(20),a(21),a(22),a(23),a(24),a(25)) d = abs(d) C Mark determinant as already found if ( d .ge. nlo .and. d .le. nhi ) s(d) = 1 C Extremely unlikely hit by chance if ( d .eq. target ) goto 300 C l = 0 C Perform up to 100 sweeps of interchange attempts do 100 it = 1, 100 lic = 0 do 30 i = 1, nn-1 do 40 j = i+1, nn C Interchange matrix elements a(i) and a(j) C Remember values to be able to restore this state if interchange C is not successful ai = a(i) aj = a(j) a(j) = ai a(i) = aj C db = d4 (a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9), C & a(10),a(11),a(12),a(13),a(14),a(15),a(16)) db = d5 (a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9), & a(10),a(11),a(12),a(13),a(14),a(15),a(16),a(17), & a(18),a(19),a(20),a(21),a(22),a(23),a(24),a(25)) db = abs(db) C Mark found determinant value if ( db .ge. nlo .and. db .le. nhi ) s(db) = 1 C Check for exact hit if ( db .eq. target ) goto 300 C Check for reduction of difference to target value if ( abs(db-target) .lt. abs(d-target) ) then C Count improvements and store index of successful sweep C (activate for diagnostic output) C l = l + 1 C lit = it C Mark current sweep as successful lic = 1 d = db else C Restore matrix to state before trial interchange a(i) = ai a(j) = aj endif C End of interchange loops 40 continue 30 continue C C Check for improvement during last sweep; C if no improvement was found try new random matrix if ( lic .eq. 0 ) goto 200 C End of loop over repeated sweeps 100 continue 200 continue C C Update best value and print (activate to watch success C during search for single target) c if ( abs(d-target) .le. abs(dglob-target) ) then c dglob = d c write (*,*) ' d=', d, ' l, lit:', l, lit c print elements of matrix with matching determinant c write (*,1000) (a(k),k=1,nn) c1000 format ( (25i3),:,/,(25i3) ) c if ( abs(dglob) .eq. abs(target) ) goto 300 c endif C C End of loop over different random matrices used to start C search for one given target value of determinant C 500 continue C C If loop ran through, no matrix with target value of determinant C was found lnot = lnot + 1 write (*,*) lnot, target C C End of loop over target values, successful search causes an C immediate jump to this label 300 continue C C Optional write of "det found" array (may be used to restart when C scanning upper range) c open ( unit=20, file='det5.dat', form='unformatted', c & status='unknown' ) c write (20) s c close ( unit=20 ) c write (*,*) ' lnot=', lnot goto 600 900 continue write (*,*) ' final lnot=', lnot do 910 i = nlo, nhi if ( s(i) .eq. 0 ) write (21,1021) i 1021 format ( i8 ) 910 continue end