smallest |cos(x)| ?

[Albert Chan]

Last month, I had proposed a better complex tan(z) for Free42

We considered HP71B math-rom implementation.

Quote:tan((x,y)) = ( sin(x)*cos(x)/d , sinh(y)*cosh(y)/d )
d is sinh(y)^2+cos(x)^2

Does d ever goes to zero ?
What is smallest cos(x)^2 ? (numerically)

To make cos(x)^2 small, x must be close to (2k+1) * pi/2
This is Python script for searching smallest |cos(x)|, for IEEE double

Code:

from gmpy2 import *
def eps(x, f=next_above):
    'find smallest |cos(x)|, x should be close to (2k+1)*pi/2'
    x2 = f(x); y = cos(x); y2 = cos(x2)
    if abs(y) < abs(y2): x,y, x2,y2, f = x2,y2, x,y, next_below
    while y*y2 > 0:
        x, x2 = x2, f(x2)
        y, y2 = y2, cos(x2)
    if abs(y) < abs(y2): return x, y
    return x2, y2

Looping (2k+1)*pi/2, for k = 1 to 1000000 step 2, k = 29 give tiniest |cos(x)|

>>> eps(const_pi()/2 * 29)
(mpfr('45.553093477052002'), mpfr('-6.1898063658835771e-19'))

Others had search for this too. From Accurate trigonometric functions for large arguments

>>> eps(mpfr('0x1.6ac5b262ca1ffp+849'))
(mpfr('5.3193726483265414e+255'), mpfr('-4.6871659242546277e-19'))

---

How to setup a Free42 program to search for this ?

[Werner]

Quick hack (corrected):

Code:

00 { 145-Byte Prgm }
01▸LBL "EPS"
02 RAD
03 LSTO "X"
04 COS
05 LSTO "Y"
06 1ᴇ-33
07 LSTO "E"
08 LASTX
09 XEQ 14
10 LSTO "X2"
11 COS
12 LSTO "Y2"
13 ABS
14 RCL "Y"
15 ABS
16 X>Y?
17 GTO 10
18 RCL "X"
19 X<> "X2"
20 STO "X"
21 RCL "Y"
22 X<> "Y2"
23 STO "Y"
24 -1
25 STO× "E"
26▸LBL 10
27 RCL "Y"
28 RCL× "Y2"
29 X≤0?
30 GTO 00
31 RCL "X2"
32 STO "X"
33 XEQ 14
34 STO "X2"
35 COS
36 X<> "Y2"
37 STO "Y"
38 GTO 10
39▸LBL 00
40 LASTX
41 ABS
42 RCL "Y2"
43 ABS
44 X<Y?
45 GTO 00
46 RCL "X"
47 RCL "Y"
48 RTN
49▸LBL 00
50 RCL "X2"
51 LASTX
52 RTN
53▸LBL 14
54 ENTER
55 XEQ "MANT"
56 STO÷ ST Y
57 RCL+ "E"
58 ×
59 END

Of course, MANT has been coded long ago, thanks to Dieter:

Code:

00 { 25-Byte Prgm }
01*LBL "MANT"
02 ABS
03 1    @ avoid upper exponent limit. ONLY needed for 9.99999998854e499 < x < 1e500
04 X<Y?
05 10^X
06 /
07 ENTER
08 X#0?
09 LOG
10 INT
11 10^X
12 /
13 1    @ multiply by 10 for x < 1
14 X>Y?
15 10^X
16 *
17 END

Then, the loop, for k=1..10000:

Code:

00 { 84-Byte Prgm }
01▸LBL "ACH"
02 1
03 LSTO "Y"
04 CLX
05 LSTO "X"
06 LSTO "K"
07 LSTO "M"
08▸LBL 02
09 0.5
10 RCL+ "K"
11 PI
12 ×
13 XEQ "EPS"
14 ABS
15 RCL "Y"
16 X≤Y?
17 GTO 00
18 R↓
19 STO "Y"
20 R↓
21 STO "X"
22 RCL "K"
23 STO "M"
24▸LBL 00
25 ISG "K"
26 X<>Y
27 10000
28 RCL "K"
29 X<Y?
30 GTO 02
31 RCL "X"
32 RCL "Y"
33 RCL "M"
34 END

Which yields the smallest value for k=953, being 8.2E-35, for x=2995.508595197867852874130465957006
If I haven't made an error, that is ;-)

Cheers, Werner

[Albert Chan]

Thanks, Werner !

I found another approach.
Instead of search min(|cos(x)|), I search min(|x - round(x,34)|), where x = (2k+1)*pi/2 = (k+0.5)*pi
Searching k all the way to 100 million, we found a winner: x = 82,425,623.5 pi

This required arbitrary precision package, with more digits of pi
We use 100 digits of pi, search 1 million of k at a time.

The code is a 1 liner

eps.bat Wrote:@gawk "BEGIN{n=1e6; while(n--) print}" | rpn =100 pi sp %1 1e6x .5+ x s d =34 - abs [ rp r + s d = - abs $s k1

rpn is a sample program, for MAPM. (+ - * are exact, rounding must be called explicitly)
It is not programmable, thus the need for gawk to produce 1 million blank lines.

=100 pi sp # p = 100-digits-of-pi
%1 1e6x .5+ x s # start = (%1*1e6 + 0.5) * p
d =34 - abs # eps1 = abs(start - round(start,34))

[ # start a pipe, for each line of input, perform below commands
rp r + s # start += p
d = - abs # eps2 = abs(start - round(start,34))
$s k1 # sort (eps1, eps2) in ascending order, then drop the big one

c:\> eps 0
8.200102230416340029960966876550294E-35

c:\> eps 82
1.532138639232766081410107492878833E-35

Then, locate x within the million k, and confirm:

c:\> rpn
=100
pi 953.5x ?68 d =-34 ?
2995.5085951978678528741304659570060000820010223041634002996096687655
2995.508595197867852874130465957006
- abs ?34
8.200102230416340029960966876550294E-35

pi 82425623.5x ?68 d =-34 ?
2.5894773325515822091636756232496249999999998467861360767233918589893E+8
2.589477332551582209163675623249625E+8
- abs ?34
1.532138639232766081410107492878833E-35

|cos(x)| = sin(ε) = ε - ε^3/3! + ... ≈ ε

Free42 Decimal:

2995.508595197867852874130465957006 COS       → -8.200102230416340029960966876550293e-35
258947733.2551582209163675623249625 COS       → 1.532138639232766081410107492878833e-35

Intel Decimal128 angle reduction code is very good

[Werner]

Running my programs to 100 million would take about 45 minutes ;-)
BTW corrected an error in EPS - it will now find that x = 82,425,623.5 pi.
A bit counterintuitive that the smallest COS happens for such a large argument value; I thought it would lose about 8 digits of precision after the argument reduction. In fact, how do you know there isn't a larger value that has an even smaller ε?
Cheers, Werner

[EdS2]

What an excellent investigation!

Does this approach to tan(z) work out as suitably accurate in this case, when the divisor is so small? (Hoping the question makes sense!)

[Werner]

While I've been able to speed it up a bit, a thought occurred to me:
instead of doing
(k+0.5)*pi, which has 2 rounding errors (one for pi, one for *), do:
->RAD(180*k + 90), which has only one, and thus would not need to examine the next numbers before or after. Depends how ->RAD is implemented of course.
Cheers, Werner

[Albert Chan]

(10-08-2021 02:16 PM)Werner Wrote:  instead of doing (k+0.5)*pi, which has 2 rounding errors (one for pi, one for *), do:
->RAD(180*k + 90), which has only one ...

"Bad" PI does not matter. (Besides, ->RAD might use the same PI)
Consider 34th digit of stored PI = 3, but should be 2.884...

x = (k+0.5)*PI
rel_error(x) = rel_error(PI) ≈ (0.116 ULP) / (PI*10^33 ULP) ≈ 3.686E-35

If x has |cos(x)| minimized, x scaled to 34 digits, is very close to integer.
Say, combined relative error = 3.7E-35 (as long as error < 5E-35, we are safe)

ulp_error(x) < 10^34 ULP * 3.7E-35 = 0.37 ULP

If we assume multiply has correct rounding, 0.37 ULP will have removed.
Searching for min(|cos(x)|, my OP eps() code can be replaced with plain cos(x)

[Albert Chan]

(10-08-2021 07:02 AM)Werner Wrote:  A bit counterintuitive that the smallest COS happens for such a large argument value;
I thought it would lose about 8 digits of precision after the argument reduction.

This give me an idea !

For IP(x) of 9 digits, calculate (pi/2*10^(34-9)), then throwaway rounded integer part

Absolute of fractional part is how close expression to integer.
No fractional part implied x - rounded(x,34) = 0

c:\> spigot -d15 -C abs(frac(pi/2*1e25+.5)-.5)
0/1
1/11
1/12
25/299
26/311
207/2476
233/2787
440/5263
673/8050
9189/109913
28240/337789
65669/785491
1472958/17618591
1538627/18404082
13781974/164851247
539035613/6447602715

x = (2k+1) * pi/2 = [1e8, 1e9)
(2k+1) = (2/pi) * [1e8, 1e9) ≈ 0.6366197724 * [1e8, 1e9)

From the convergents, 2k+1 = 164851247 satisfied.
3*(2k+1) = 494553741 also valid, but produced error of 3ε

82425623.5 PI * COS         → 1.532138639232766081410107492878833e-35
247276870.5 PI * COS       → -4.596415917698298244230322478636498e-35

Best semi-convergent is not as good: 2k+1 = 18404082 + 164851247*3 = 512957823

256478911.5 PI * COS       → -5.589328359211609676578079200248261e-34

[ijabbott]

(10-08-2021 07:02 AM)Werner Wrote:  Running my programs to 100 million would take about 45 minutes ;-)
BTW corrected an error in EPS - it will now find that x = 82,425,623.5 pi.
A bit counterintuitive that the smallest COS happens for such a large argument value; I thought it would lose about 8 digits of precision after the argument reduction. In fact, how do you know there isn't a larger value that has an even smaller ε?
Cheers, Werner

There should be infinitely many such values for x producing a smaller ε, shouldn't there?

[Albert Chan]

(10-09-2021 01:15 PM)ijabbott Wrote:  There should be infinitely many such values for x producing a smaller ε, shouldn't there?

The question is not if there are infinite smaller ε's, but whether ε^2 will underflow.

Going from IP(x) of 9 diigits to 10, it is 10 times harder to get smaller ε.
On the other hand, we have 10 times many more candidates to try.

Both effect tends to cancel out, giving min(ε) not much smaller than machine epsilon.
(this is a conjecture, but experiments tends to support it ...)

(10-07-2021 06:40 PM)Albert Chan Wrote:  Free42 Decimal:

2995.508595197867852874130465957006 COS      → -8.200102230416340029960966876550293e-35
258947733.2551582209163675623249625 COS      → 1.532138639232766081410107492878833e-35

Using convergents / semi-convergents idea, for x < 1E50, these two have smaller ε:

9.510911580848780648392560778912209e40 COS  → -1.328353856527533918552271223903534e-35
9.341730789500356812974471132146251e46 COS  → -1.028415848209791685669808767152452e-35

To show how hard to get similar sized ε, this is true x, before rounded.

30274171828040719778093696476072003369775.5 * pi
= 95109115808487806483925607789122090000000.00000000000000000000000000000000001328​3538...

29735652643654715492245370225672550765119900844.5 * pi
= 93417307895003568129744711321462509999999999999.99999999999999999999999999999999​9989...

[Albert Chan]

(10-09-2021 02:09 PM)Albert Chan Wrote:  Both effect tends to cancel out, giving min(ε) not much smaller than machine epsilon.
(this is a conjecture, but experiments tends to support it ...)

Perhaps Loch's theorem explain this.

Quote:this can be interpreted as saying that each additional term in the continued fraction representation of a "typical" real number increases the accuracy of the representation by approximately one decimal place.

When I did the convergents, steps required to get q = 2k+1 is roughly the same.
If valid 2k+1 is big, convergents also start with big denominator.

Relative error of final convergent tends towards the same ballpark.

Example, IP(x) of 47 digits. Convergents of y = (pi/2) / 10^(47 - 34):

0/1
1/6366197723675
1/6366197723676
5/31830988618379
...
933595789363114285779172744945397/5943455389076782359664376669525821535544153034
9341730789500356812974471132146251/59471305287309430984490740451345101530239801689

y = 1.570796326794896619231321691639751 ... * 10^-13
p = 9341730789500356812974471132146251
q = 59471305287309430984490740451345101530239801689

x = p * 10^13                   ≈ 9.341730789500356812974471132146251e46
ε = abs(q*y - p) * 10^13 ≈ 1.028415848209791685669808767152452e-35

Update:
Brute force proof, for Free42-Decimal, |cos(x)| > 10^-37
see https://www.hpmuseum.org/forum/thread-17...#pid153539

[Albert Chan]

Trivia: 2 convergents cannot have both even denominator (or both even numerator)

Code:

def genConvergents(coefs):
    'Return generator for convergents of n/d'
    n0, d0, n1, d1 = 0, 1, 1, 0
    for coef in coefs:
        n0, n1 = n1, n0 + coef * n1
        d0, d1 = d1, d0 + coef * d1
        yield (n1, d1)

From the loops, assume d0 is odd, d1 is even (same reasoning for numerator)

d1 ← d0 + coef*d1 = odd + even = odd

With this trivia, our search for "best" denominator = 2k+1 guaranteed exist.

[Werner]

(10-08-2021 07:54 PM)Albert Chan Wrote:  

(10-08-2021 02:16 PM)Werner Wrote:  instead of doing (k+0.5)*pi, which has 2 rounding errors (one for pi, one for *), do:
->RAD(180*k + 90), which has only one ...

"Bad" PI does not matter. (Besides, ->RAD might use the same PI)
Consider 34th digit of stored PI = 3, but should be 2.884...

x = (k+0.5)*PI
rel_error(x) = rel_error(PI) ≈ (0.116 ULP) / (PI*10^33 ULP) ≈ 3.686E-35

If x has |cos(x)| minimized, x scaled to 34 digits, is very close to integer.
Say, combined relative error = 3.7E-35 (as long as error < 5E-35, we are safe)

ulp_error(x) < 10^34 ULP * 3.7E-35 = 0.37 ULP

If we assume multiply has correct rounding, 0.37 ULP will have removed.
Searching for min(|cos(x)|, my OP eps() code can be replaced with plain cos(x)

Somehow, this is not true.. the division by 2 (or mult. by 0.5 ..) may introduce an error of one ULP, due to unbiased rounding, as is the case for k=0.
For the following values of k , (k+0.5)*pi is not the same as k*180+90 ->RAD, and the latter is correct:
k=0, 2, 8, 15, 18, 25, ..
consequently, (k+0.5)*pi does not yield the lowest |cos(x)|, but the ->RAD version does.

Cheers, Werner

[Albert Chan]

(10-11-2021 01:16 PM)Werner Wrote:  For the following values of k , (k+0.5)*pi is not the same as k*180+90 ->RAD, and the latter is correct:
k=0, 2, 8, 15, 18, 25, ..
consequently, (k+0.5)*pi does not yield the lowest |cos(x)|, but the ->RAD version does.

Assuming we already tried k=0,1,2 (i.e, IP(x) of single digit), difference does not matter.

We are now searching for exact(x) *close* to round(x,34).
For this, only error of concern is inaccuracy of PI, (rel_err < 0.37 ULP)
(Here, ULP ≡ 1E-34)

If |cos(x)| is minimized, both ways of calculating x will match.
If they are different, it is guaranteed |cos(x)| is not minimized.

Example, for IP(x) of 2 digits:

c:\> spigot -d3 -C "abs(pi/2*1e32 rem 1)"
0/1
1/6
1/7
15/104

7 2 / PI *       → 10.99557428756427633461925184147826
7 90 * ->RAD → 10.99557428756427633461925184147826
COS               → -9.469009289781287037341230607307734e-35

---

From Free42 core_helpers.cc, deg_to_rad(x) = x / (180 / PI)

180/pi = 57.29577951308232087679815481410517 0332 ...

Rounded to 34 digits, this has relerr ≈ 0.0332 / 0.573 ≈ 0.058 ULP

This is much more accurate than stored PI (relerr ≈ 0.369 ULP)
If it had used flipped constant:

pi/180 = 0.01745329251994329576923690768488612 7134 ...

relerr ≈ (1-0.7134) / 0.1745 ≈ 1.64 ULP ... very bad

Comment:
Here, relerr calculations ignored sign, but sign may be important.
I use this convention for sign of error: exact = approx + error

PI is over-estimated, relerr(PI) = −0.369 ULP
180/PI under-estimated, relerr(180/PI) = +0.058 ULP

But, ->RAD() use division of constant, we subtract relerr instead of add.
Or, we can think of multiply by inverse of constant, with relerr of −0.058 ULP

Since both formula have relerr of same sign, ->RAD() always better.

[Albert Chan]

Just for fun, this formula is more accurate than ->RAD()

(k+0.5)*PI ⇒ (k+0.5)*29/(29/PI)

29/pi = 9.230986699329929474595258275605832 997998659 ...

Rounded to 34-dgits, relerr ≈ (1-0.99800)/0.92310 < 0.00217 ULP

---

With flipped constant, this is even better.

(k+0.5)*PI ⇒ (k+0.5)*399*(PI/399)

pi/399 = 0.007873665798470659745520409481903516 000494158 ...

Rounded to 34-dgits, relerr ≈ 0.0004941/0.7874 < 0.00063 ULP