-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathconsts.f
299 lines (299 loc) · 9.67 KB
/
consts.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
C Alternative search function for constants
C This is an alternative search function that is used to supplement,
C but not replace, the query implemented in the database.
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C ----------------------------------------------------------------------
C Copyright 2024 Hugo Pfoertner
C
C Licensed under the Apache License, Version 2.0 (the "License");
C you may not use this file except in compliance with the License.
C You may obtain a copy of the License at
C
C http://www.apache.org/licenses/LICENSE-2.0
C
C Unless required by applicable law or agreed to in writing, software
C distributed under the License is distributed on an "AS IS" BASIS,
C WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
C See the License for the specific language governing permissions and
C limitations under the License.
C
C Version history:
C 2024-05-21 A-number exclusion selectable by shortcut
C 2024-02-20 Initial version derived from oeiseach.f
C
C Local size parameters
C ndmax ... maximum digits in constant
C mf ... 1 + defined conversion formulas
parameter (ndmax=21, mf=13, lcline=160, laline=500)
real*16 q, one, two, three, half, othird, pi, e
character c*50, eform*10, eout*35, T(mf)*7
C Short names for conversion formulas
data t /'1*c', '2*c ', 'c/2 ', '3*c ', 'c/3 ',
& 'c^2 ', 'c^(1/2)',
& '1/c ', 'Pi*c ', 'e*c ', 'log(c)', 'exp(c)','Cl(n)-c'/
integer w(ndmax,mf)
C 16 byte floats to read data from file "stripped"
real*16 r(500), rvalue
integer*8 timprv, timnow
C Lines in file "stripped" are assumed not to exceed 500 characters
character cline*(lcline), aline*(laline), prefil*23
logical matchr, pexist, axclud
external matchr
C Name of file with information from the preceding call
prefil = 'prev_call_oeisearch.txt'
C
C Floating point constants computed with length 128 bits
zero = real(0,16)
one = real(1,16)
two = real(2,16)
three = real(3,16)
half = one/two
othird = one/three
pi = (two+two)*atan(one)
e = exp(one)
C
C Exclusion indicator for A-numbers
naex = 0
10 continue
C Get full command line
call get_command (cline, lcmd )
C
C count of significant nonzero digits
C Position of decimal point (DP).
C Use backward search to exclude occurrences of "." in program name
ndp = index(cline,'.',.true.)
C first blank right of DP
nb1 = index(cline(ndp:),' ') + ndp - 1
C
C write termination character to rear end of line to
C enable reading by Fortran's list directed input
cline (lcline:lcline) = '/'
C extract optional excluded A-number
C check for character starting A number exclusion
i = index(cline, '!')
nex = i
naex = 0
axclud = .false.
if ( i .ne. 0 ) then
C check whether ! is the last character in input
if (len_trim(cline(1:lcline-1)) .eq. i ) then
C try to read A-number found in previous call
timnow = time8()
inquire (file=prefil, exist=pexist)
if (pexist) then
open (unit=9, file=prefil, form='formatted', status='old')
read (9,*) timprv, nax
close (unit=9,status='delete')
if (timnow-timprv .lt. 120) then
naex = nax
axclud = .true.
ios = 0
endif
endif
else
j = max(1,index(cline(i:), 'A' ))
read ( cline(i+j:),*,iostat=ios) naex
endif
C replace "!" by another termination character to stop reading
C of search item list at this point
cline(i:i) = '/'
if ( ios .ne. 0 ) then
write (*,*)
& 'unreadable A-number, continue without exclusion'
else
C silently limit A-number to plausible range
if ( naex .lt. 1 .or. naex .gt. 999999 ) then
naex = 0
else
write (*,1003) naex
1003 format ( 'Hitting A', I6.6,
& ' will not terminate the search')
endif
endif
endif
C Constant may also be terminated by "!" instead of a blank
if (nex .gt. 0) nb1 = min(nb1, nex)
C
C Position of last blank before DP
nb0 = index(cline(1:ndp-1),' ',.true.)
if (nb0 .eq. 0) stop 'Missing separator before number'
C first nonzero digit left of decimal point
nd0 = -1
do 13 k = nb0+1, ndp-1
if ( cline(k:k) .ne. '0' ) then
nd0 = k
goto 14
endif
13 continue
14 continue
if (nd0 .ne. -1) then
nsig = nb1 - nd0 - 1
else
C no significant digits left of DP
C first nonzero digit after DP
nz1 = nb1
do 15 k = ndp+1, nb1-1
if (cline(k:k) .ne. '0') then
nz1 = k
goto 16
endif
15 continue
16 continue
nsig = nb1 - nz1
if (nsig .eq. 0) stop 'Only digits zero found'
endif
C
C read floating point number
read (cline(nb0:nb1-1),*) q
C
C truncate input to 20 significant digits
nsig = min (nsig, ndmax)
C Variable format from which only the mantissa part is used
C The E-format always produces a mantisse of the form 0.xxxxx...
eform='(e30.xxe4)'
write(eform(6:7),'(i2.2)') nsig
write (eout,fmt=eform) q
call wrivec(eout,nv,w(1,1))
write (eout,fmt=eform) two*q
call wrivec(eout,nv,w(1,2))
write (eout,fmt=eform) half*q
call wrivec(eout,nv,w(1,3))
write (eout,fmt=eform) three*q
call wrivec(eout,nv,w(1,4))
write (eout,fmt=eform) othird*q
call wrivec(eout,nv,w(1,5))
write (eout,fmt=eform) q*q
call wrivec(eout,nv,w(1,6))
write (eout,fmt=eform) sqrt(q)
call wrivec(eout,nv,w(1,7))
write (eout,fmt=eform) one/q
call wrivec(eout,nv,w(1,8))
write (eout,fmt=eform) pi*q
call wrivec(eout,nv,w(1,9))
write (eout,fmt=eform) e*q
call wrivec(eout,nv,w(1,10))
write (eout,fmt=eform) exp(q)
call wrivec(eout,nv,w(1,11))
write (eout,fmt=eform) log(q)
call wrivec(eout,nv,w(1,12))
write (eout,fmt=eform) ceiling(q)-q
call wrivec(eout,nv,w(1,13))
C
C Start of search in file "stripped"
open (unit=10,file='stripped',form='formatted',
& status='old',iostat=ios)
if ( ios .ne. 0 ) stop 'file stripped not found'
C skip header lines in file stripped
1000 format (A)
read (10,1000) c
read (10,1000) c
read (10,1000) c
read (10,1000) c
C get largest representable number to preset terms
C > 10^4000 for 16 byte floats, so no risk of interfering
C with any DATA item
rvalue = 1.0D0
rvalue = huge(rvalue)
C This dummy write inibits segfault in Windows executable
C not understood
write (cline,*) 'Huge:', rvalue
write (*,*) 'Searching in file stripped ...'
100 continue
read (10,1000,end=200) aline
read (aline,'(1X,i6.6)') numa
C If shortcut exclusion is active,
C skip all A-numbers <= A-number of previous success
if (axclud .and. numa .le. naex) goto 100
lastc = index ( aline, ',', .true. )
aline(lastc:lastc) = '/'
C preset terms before list directed read
r = rvalue
read (aline(10:lastc),*) r
C
C Loop over transformed inputs
do 101 nof = 1, mf
C call comparison function
if ( matchr(r, nv, w(1,nof), rvalue, nrel) ) then
C write (*,1006) nof, offset(nof), gcd(nof), w(1:n,nof)
write (eout(1:nv),1101) w(1:nv,nof)
1101 format(21 (i1,:))
write (*,1100) T(nof), eout(1:nv), numa, nrel
1100 format ( 'Input ', A,': match ', A,' found in A', i6.6,
& ' at position ', i0 )
C Termination after first match
if (numa .ne. naex) goto 900
endif
101 continue
goto 100
200 continue
C
call cpu_time (etime)
write (*,1099) etime
1099 format ('No success; search time', F6.2, ' s')
goto 999
900 continue
if ( numa .gt. 0 .and. numa .ne. naex ) then
call cpu_time (etime)
write (*,1900) numa, etime
1900 format ( 'URL: https://oeis.org/A', i6.6,
& /, 'End of search; time', F6.2, ' s' )
C Write information on success to file
open (unit=9, file=prefil, form='formatted', status='unknown')
write (9,*) time8(), numa
endif
999 continue
close (unit=10)
end
C
subroutine wrivec(eout, nw, d)
integer nw
character*(*) eout
integer d(*)
idp = index(eout,'.')
iex = index(eout,'E')-2
C write (*,*) idp, iex, eout
k = 0
do 10 i = idp+1, iex
k = k + 1
read(eout(i:i),'(i1)') d(k)
10 continue
nw = k
C write (*,1000) d(1:k)
C1000 format (i0, 19(',',i0,:))
end
C
logical function matchr ( r, n, s, rvalue, nrel )
C Find contiguous match of search vector s in original DATA
C vector R.
C Modified version: Exit when sequence terms < 0 or > 9 are found
real*16 r(*), rvalue, dr, diglim, zero
integer*4 s(n)
C Difference between 8 byte integer and 16 byte float will
C be exact for the possible range of search values.
parameter (dr=5.0D-1, diglim=9.5D0, zero=0.0D0)
nrel = 0
C Empirical limit: not more than 120 data items in single line of
C file "stripped"
do 110 i = 1, 120
if ( r(i) .eq. rvalue
& .or. r(i) .gt. diglim .or. r(i) .lt. zero ) then
C signal no success if loop reaches empty part of vector
matchr = .false.
return
endif
C Check if first searched item is present in list
if ( abs(r(i) - s(1)) .lt. dr ) then
C Check if remaining searched items are at the
C subsequent positions of the list
do 120 j = 2, n
C first non-match exits from inner loop
if (abs(r(i+j-1) - s(j)) .gt. dr ) goto 110
120 continue
nrel = i
matchr = .true.
return
endif
110 continue
C End of function matchr
end