mini challenge: find the smallest cosine of an integer

[EdS2]

.
Naturally, we should use radians for this, and we're interested in the absolute value of the cosine. Some calculators will have a maximum integer that they are happy to take the cosine of, which will of course affect your result. But there are no prizes, so start your search - whether manually, programatically, or mathematically. Show your working!

Spoiler padding follows...

.
..
.
...
.
..
.
....
.
..
.
...
.
..
.
.....
.
..
.
...
.
..
.
....
.
..
.
...
.
..
.

cos(0) is the first, and the worst...
cos(1) is better, and could hardly be otherwise...
cos(2) is better yet
cos(5) is a new record
cos(8) is another winner
cos(11) is pretty good! (And surely no coincidence)

[Albert Chan]

For getting smallest |cos(x)|, we are only using angle-reduction code of cosine function.
Tiny |cos(x)| implied x ≈ (2k+1)*(pi/2)

|cos(x)| = |sin(x - (2k+1)*(pi/2))| ≈ |(x - (2k+1)*(pi/2))|

Example, |cos(11)| ≈ |11 - 7*(pi/2)| ≈ 0.0044257

So, it depends how good angle-reduction code is. Most calculators are not good.

Free42-Decimal happened to have very good angle-reduction code (at the cost of *huge* table)
This is the smallest |cos(x)|, for x < 1E50

9.341730789500356812974471132146251e46 COS ABS → 1.028415848209791685669808767152452e-35

Above x (47 digits integer) is derived from convergents of y = (pi/2) / 10^(47 - 34)
see https://www.hpmuseum.org/forum/thread-17...#pid153032

[PeterP]

Brute force, HP 41

Code:


LBL 'SC'
STO 01...... that is the delta to a multiple of pi/2 or 3pi/2
pi
2
/
STO 02...... current value to check
LBL 00
  pi
  ST+ 02
  RCL 02
  FRC
  CHS
  1
  +
  ABS
  RCL 01
  X<Y? ...... is the current test value further away from an integer than the chosen delta
GTO 00
  BEEP
  VIEW 02
  STOP

0.01 XEQ 'SC' yields 10.994 hence the integer 11
0.001 XEQ 'SC' yields 15,253.9997 hence the integer 15,254
0.0001 XEQ 'SC' yields 24,461.99999 hence the integer 24,462

Upon which the accuracy limit of the HP 41 is reached.

Cheers

PeterP

[Albert Chan]

Limiting x to integer ≈ (2k+1)*(pi/2)

pi/2 ≈ x / (2k+1)

Pick pi/2 convergents with odd denominator, then numerator = x

c:\> spigot -d22 -C pi/2
1/1
2/1
3/2
11/7
344/219
355/226
51819/32989
52174/33215
260515/165849
573204/364913
4846147/3085153
5419351/3450066
37362253/23785549
42781604/27235615
122925461/78256779
411557987/262005952
534483448/340262731
2549491779/1623056876
3083975227/1963319607
17969367914/11439654911
21053343141/13402974518
881156436695/560961610149
902209779836/574364584667

HP50g:
902209779836 COS       → -5.28216132855E-13 (red part should be 3683)

[Albert Chan]

If calculator had 12-digits precision, and we would like to search x of 14 digits:

pi/2/100 ≈ (x/100) / (2k+1)

Pick (pi/2/100) convergents with odd denominator, then numerator = x/100
Note: To get x of 14 digits, numerator must be 12 digits.

c:\> spigot -d20 -C pi/2/100 | tail -5
1324880093/84344486322
8539351431/543631996417
9864231524/627976482739
205823981911/13103161651197
833160159168/53040623087527

HP50g:
833160159168E2 COS       → -2.94461252684E-14 (red part should be 2017582)

[EdS2]

Thanks both, that's two very different approaches!

Peter, I like the simplicity and directness of your program. It's a good trick, I think, to look for n close to odd multiples of pi/2, and leave the cosines out of the loop.

But there's something curious about your result. You get 24462 as the first n to have a cos less than 0.0001. And yet
cos(24462) ≈ 0.0112.

I suppose that adding pi some twenty thousand times might have caused a precision problem, but it seems a bit early for that to have happened. (In fact, it's about right: five digits to the left of the point and looking for small values in five digits to the right.)

I think your program might do somewhat better by adding (pi-3) instead of pi.

The Online Encylopedia of Integer Sequences is quite handy. For this investigation, these two entries are of some use:
A004112 Numbers k where |cos(k)| (or |cosec(k)| or |cot(k)|) decreases monotonically to 0; also |tan(k)|, |sec(k)|, |sin(k)| increases.
A096456 Numerators of convergents to Pi/2.

From the first, we find a link to this list of record-holders for small cos, and from that we see:
cos(24129) ≈ 0.002375894172
cos(24484) ≈ -0.002345749902
which I think must only be coincidentally close to your 24462.

But perhaps more notably we see
cos(40459) ≈ 0.000989256
which ought to be the first one less than 0.001

Albert, I am not at all surprised that continued fractions and convergents would come into a mathematical solution! But your method for dealing with trailing zeros to find answers for calculators with limited precision is very nice.

I wonder if the OEIS would like to add a series for the numerators of pi/2 convergents where the denominators are odd.

I note also that Free 42 shows its strengths again, quite wonderfully.

(I note also that wolfram alpha online needs to be given integers with lots of trailing zeros, to compute cosine with adequate precision: giving it those large inputs in floating form returns zero for cosine.)

[Albert Chan]

(10-21-2021 09:44 AM)EdS2 Wrote:  A004112 Numbers k where |cos(k)| (or |cosec(k)| or |cot(k)|) decreases monotonically to 0

Since it is impossible to input x to infinite precision, |cos(x)| does not decrease to zero.
With limited precision, |cos(x)| is not much smaller than machine-epsilon.

Quote:cos(24129) ≈ 0.002375894172
cos(24484) ≈ -0.002345749902
cos(40459) ≈ 0.000989256

FYI, all 3 x's comes from semi-convergents of pi/2:

344/219
355/226
51819/32989

24129 = 344 + 355*67
24484 = 344 + 355*68
40459 = 344 + 355*113

Quote:I wonder if the OEIS would like to add a series for the numerators of pi/2 convergents where the denominators are odd.

No need. Just use convergents of pi instead

pi/2 ≈ x / (2k+1)
pi ≈ (2x) / (2k+1)

With convergents fully reduced, even numerator guaranteed odd denominator.
So, just pick convergents/semi-convergents of pi, with even numerator = 2x

Note: convergents of pi/2 may have shifted to semi-convergents of pi.
Example, convergents of pi: 333/106, 355/113

(333+355) / (106+113) = 2*(344/219)       // (344/219) is a convergent of pi/2

[Albert Chan]

(10-21-2021 12:01 PM)Albert Chan Wrote:  With convergents fully reduced, even numerator guaranteed odd denominator.

2 neighboring convergents, p1/q1, p2/q2, we have

p1*q2 - p2*q1 = ± 1       // or, gap = |p1/q1 - p2/q2| = 1/|q1*q2|

If g = gcd(p1,q1) ≠ 1, RHS must be multiples of g, contradict above identity.
Thus, convergents always fully reduced.

[EdS2]

Peter, I notice something else about your program: it seeks numbers which are just smaller than an integer, but in our case numbers which are just larger are also of interest. I think it should only take a line or two to address this.

[EdS2]

Thinking a little more about Peter's program, if we separately accumulate the 3 and the pi-3, and if we carry a 0.5 from the accumulation of fractions whenever it exceeds 0.5, into the accumulation of integers (now an accumulation of half-integers), we'll find that we lose very little accuracy as proceed, which might extend the useful range of the search.

Indeed, I think we can do a little better, by also accumulating the even smaller fraction which holds the difference between the limited precision value of pi we're using and the real value of pi. We need to add this accumulation of tiny pieces to the accumulation of fractions when comparing to see how close we are to an integer. I don't think we'll need to split this tiny accumulation into big and small parts.

Before I decided to trim the fractional accumulation by 0.5, I was going to trim it whenever it exceeded 1.0. With this approach, whenever we go over 1.0, we can as a next step accumulate 6 times our fraction, because we know the fraction is less than 1/7. So that could be a speedup. If that works, it could I think be modified for the more accurate approach of trimming at 0.5.

The more mathematically inclined reader might be interested in this related paper I found:
Rosenholtz, I. (1999). Tangent Sequences, World Records, π, and the Meaning of Life: Some Applications of Number Theory to Calculus. Mathematics Magazine, 72(5), 367. doi:10.2307/2690793

Quote:the author knows of no other instance in which knowing several million digits of pi is actually useful

Quote:To summarize: If we desire to find |tan(n)|/n records, we need only look at numerators of those convergents of the continued fraction expansion of pi/2 having odd denominators

[EdS2]

Albert, thanks for your further thoughts. On this comment

Quote:Since it is impossible to input x to infinite precision, |cos(x)| does not decrease to zero

I think we are in danger of mixing mathematical results with the results of calculation. The statement

Quote:A004112 Numbers k where |cos(k)| (or |cosec(k)| or |cot(k)|) decreases monotonically to 0

surely stands, in mathematics, even if it fails in calculation. (It fails because our calculations have finite resources.)

Is it certain, though, that there is no calculator input of a positive integer which causes the calculator to return a cos of zero? It's not clear to me, presently.

[Albert Chan]

(10-23-2021 06:53 AM)EdS2 Wrote:  Is it certain, though, that there is no calculator input of a positive integer which causes the calculator to
return a cos of zero? It's not clear to me, presently.

Let's explore what condition needed to have calculator return |cos(x)| of 0.0
Assumed correct angle-reduction code, and no small-differences-incorrectly-set-to-zero “feature"

Again, using 12-digits calculator, and integer x of 14 digits as example.

Instead of convergents of pi/200, we do 200/pi
Valid convergent/semi-convergent must have odd numerator, and denominator of 12 digits.

c:> spigot -d20 -c -l 200/pi
63;1,1,1,23,36,1,1,1,8,2,1,52,13,4,2,6,1,20,4,63

These are convergents of last 3 CF coefs:

20: p0/q0 = 13103161651197/205823981911
04: p1/q1 = 53040623087527/833160159168
63: p2/q2 = 3354662416165398/52694914009495

gap = |p1/q1 - p2/q2| = |p1*q2-p2*q1| / (q1*q2) = 1/(q1*q2)

We expected next convergent much better than current convergent.
Rough estimate, |200/pi - p1/q1| < gap = 1/(q1*q2)

But, gap error is really shared between 2 convergents.
Perhaps we increase q2 a bit: q2 = q0 + q1*CF < q1 * (CF+1)

|200/pi - 53040623087527/833160159168| ≈ 2.250E-26

1/(833160159168^2 * (63+1)) ≈ 2.251E-26

|cos(x)|                   // x = 100*q1, with odd p1
≈ |x - p1*(pi/2)|
= (100*q1) * |1 - (p1/q1)/(200/pi)|
≈ (100*q1) * 1/(q1^2*(CF+1)) / (200/pi)
= (pi/2) / (q1*(CF+1))
≈ 2.946E-14            // |cos(833160159168E2)| ≈ 2.945E-14, estimate not too bad

Note that factor 100 get cancelled. We can generalize, for big x, with valid convergent p1/q1 (*)

|cos(x)| ≈ (pi/2) / (q1*(CF+1)) ≈ (pi/2) / q2

To have |cos(x)| underflow to 0.0, q1*(CF+1) overflowed
For a 12-digits calculator that overflow if decimal exponent reached 500

(CF+1) > 10^(500-12) = 10^488

---

12-digits calculator, x all the way to 500 digits, x of 395 digits has smallest |cos(x)|

Free-42 Decimal:
0.700035629422e395 COS ABS → 2.056399907114786850992546180602331e-15

For this number, CF = 1090, estimate formula is pretty close

(pi/2) / (700035629422 * (1090+1)) ≈ 2.057e-15

(*) if p1/q1, p2/q2 are not neighboring convergents, we need another estimate formula.
Let pm/qm and p2/q2 be neighboring convergents.

p1*q2 - p2*q1 = (p2-k*pm)*q2 - p2*(q2-k*qm) = k*(p2*qm - pm*q2) = ±k
gap = |p1/q1 - p2/q2| = k / (q1*q2)

→ |cos(x)| ≈ (pi/2) / (q2/k)

Example, x of 137 digits had a higher CF = 2573, yet |cos(x)| is not small
Note: we reduced numerator of (2e125/pi), without disturbing its parity.

c:\> spigot -d28 -C "2e125/pi mod 2" | tail -3
17206041685/23300343711
202744600696/274555819883
521679063632493/706455424902670

Convergent with 12-digits denominator has even numerator.
Best we could use is semi-convergent, with k = 2570

706455424902670 - 274555819883 * 2570 = 846967803360

Estimated |cos(x)| = (pi/2) / (706455424902670/2570) ≈ 5.714e-12

Free42-Decimal:
846967803360e125 COS ABS → 5.714370521060907484280466201567459e-12

[PeterP]

Ever since your first very insightful post I had been playing around with both of these ideas as indeed, the error accumulates much quicker than I had originally thought.

Alas, things get really slow when you add the test for 0.5. And just doing Pi-3 actually only gives you about 1-2 extra digits and the error still accumulates quickly.

I’m on my phone now but can post the improved code thanks to your insights later but would not advise using it - even in iOS i41CX it runs pretty slowly….

Thank you for the nice mini-challenge. Thank you Albert for super impressive and instructive math, I love reading your posts and learn from you. Might have been nice to get someone from the 71b genius group share a code as well.

Thank you again Ed for posting it!

[EdS2]

Thanks Peter, that's very good to hear!

I would be interested to see your code. I have my 15C, 35S, and 30b ready, and all of those in principle could be programmed up for a search. (But instead I find myself playing around with Basic. And not even on a 71B.)

And thanks Albert, again: I now have a glimmering of understanding that we can't get an actual zero.

[ijabbott]

(10-24-2021 08:27 AM)EdS2 Wrote:  And thanks Albert, again: I now have a glimmering of understanding that we can't get an actual zero.

I think it's partly a quality of implementation issue. I wouldn't be surprised if there is some calculator out there that would return 0 as the cosine of some integral number of radians, particularly when operating outside the domain of exact, consecutive integer representations. For example, when a decimal, floating point number represents a continuous span of at least 10, some points in that continuous span would have a cosine of 0, so a poor quality calculator might give up and return 0 (or some other value in the range [-1,1]) for all such numbers!

[Albert Chan]

eps.bat Wrote:@spigot -C "2e%1/pi mod 2" | awk -vFS="/" -ve="e%1" "length($2)>34 {print 1.5708/$2, q1 e; exit} {q1=$2}"

Running this 1 liner overnight, x all the way to 6145 digits, we have a brute force proof

Free42-Decimal: |cos(x)| > 10^-37, any real x

if p1/q1 is valid convergent (odd p1, neighboring p2/q2), then |cos(x)| < ε = (pi/2)/q2
The requirement of q1 being 34 digits is not necessary.

Example, both below suggested the same x (note that 2nd mantissa has 33 digits, not 34)

c:\> eps 5981      → 1.60573e-37 4272221537462525975023616239425960e5981
c:\> eps 5982      → 1.60509e-37 427222153746252597502361623942596e5982

427222153746252597502361623942596e5982 COS ABS → 1.605082938363676519340189975883263e-37

min(|cos(x)|) for integer x (also, for real x)

c:\> eps 1381      → 1.03065e-37 2344813655066356855719930664718056e1381

2344813655066356855719930664718056e1381 COS ABS → 1.030557387629248882465543827418861e-37

eps.bat even work with negative exponent (below is min(|cos(x)|) for non-integer x)

c:\> eps -25       → 1.55097e-35 2589477332551582209163675623249625e-25

2589477332551582209163675623249625e-25 COS ABS → 1.532138639232766081410107492878833e-35

---

eps.bat only get right a little over half the time.
Neighboring convergents cannot both have p1,p2 even → p1 is more likely odd.
see https://www.hpmuseum.org/forum/thread-17...#pid153072

We need to check actual p1/q1, for odd p1 (or, just try cos(x))
For even p1, we may have to settle for semi-convergents, or p0/q0 (p0 guaranteed odd)
(it is highly unlikely to produce small |cos(x)| though, see my previous post)

pq.bat Wrote:@spigot -C "2e%1/pi mod 2" | awk -vS="/" "length-index($0,S)>34 {exit} {print}" | tail -2

Example: this get wrong eps (being conservative is OK; we don't miss anything)

c:\> eps 4639       → 2.06902e-38 9308532438209917461067659354862169e4639
c:\> pq 4639
1528624876453580931164937665986407/8775804805425001620222943906879558 // p0/q0
1621418726123738114106180500221048/9308532438209917461067659354862169 // p1/q1

p1 is even, we have to settle for p0/q0: cos(x = q0*1E4639) ≈ 1.68729e-34
However, this suggested sin(x = q1*1E4639) is very close to 0.

9308532438209917461067659354862169e4639 SIN → -2.069013989544476107762446187843978e-38

[PeterP]

Ever since your first very insightful post I had been playing around with both of these ideas as indeed, the error accumulates much quicker than I had originally thought.

Alas, things get really slow when you add the test for 0.5. And just doing Pi-3 actually only gives you about 1-2 extra digits and the error still accumulates quickly.

I’m on my phone now but can post the improved code thanks to your insights later but would not advise using it - even in iOS i41CX it runs pretty slowly….

Thank you for the nice mini-challenge. Thank you Albert for super impressive and instructive math, I love reading your posts and learn from you. Might have been nice to get someone from the 71b genius group share a code as well.

Thank you again Ed for posting it!

[Albert Chan]

(10-24-2021 06:56 AM)PeterP Wrote:  Alas, things get really slow when you add the test for 0.5. And just doing Pi-3 actually only gives you about 1-2 extra digits and the error still accumulates quickly.

A novel way is use SIN(X)

>PI, SIN(PI) @ REM sum = 3.141592653589793238462643 (25 digits of pi)
3.14159265359 -2.06761537357E-13

This program get numerator of convergents of pi
With internal 25 digits of stored pi, all convergents correct

Building of convergents is fast, using formula, P = P0 + k*P1
(P's are all integers, no error is introduced building P's)

10 P0=0 @ P1=1 @ P=P0 @ S=2
20 P=P+P1 @ S2=ABS(SIN(P))
30 IF S>S2 THEN S=S2 @ DISP P @ P0=P1 @ P1=P @ P=P0
40 IF P<1.E+12 THEN 20

>RUN
1
3
22
333
355
103993
104348
208341
312689
833719
1146408
4272943
5419351
80143857
165707065
245850922
411557987
1068966896
2549491779
6167950454
14885392687
21053343141

sin(2x) = 2*sin(x)*cos(x). If |cos(x)| is small, so does |sin(2x)|
We pick convergents/semi-convergents with even numerator.

2x = 1+3, 1+3*3, 1+3*5, 22, 333+355, 333+355*3, 333+355*5 ... 333+355*291, 104348, ...

x = 2, 5, 8, 11, 344, 699, 1054 ... 51819, 52174, ...

[PeterP]

Great stuff Albert! Thank you for sharing.

[Albert Chan]

To complete symmetry, this is min(|cos(x)|) for Free-42 binary, if angle reduction code is correct.

eps.bat Wrote:@spigot -C "2*2^%1/pi mod 2" | awk -vFS="/" -ve="*0x1p%1" "$2>=2^53 {print 1.5708/$2, q1 e; exit} {q1=$2}"

Just a change of constants needed for eps.bat
This take less than 1 minute, x upto 971+53 = 1024 bits integer.

c:\> (for /L %I in (-52,1,971) do @(eps %I)) | gsort -g | head
4.69203e-19 6381956970095103*0x1p797
6.19292e-19 3205513981387887*0x1p-46
6.20075e-19 6411027962775774*0x1p-47
9.38406e-19 6381956970095103*0x1p798
1.23858e-18 3205513981387887*0x1p-45
...

The list is suggestion only, confirmation necessary. First does give minimum |cos(x)|
see Accurate trigonometric functions for large arguments

>>> from gmpy2 import *
>>> x = mul_2exp(6381956970095103,797) # = 0x1.6ac5b262ca1ffp+849
>>> x
mpfr('5.3193726483265414e+255')
>>> cos(x)
mpfr('-4.6871659242546277e-19')