smallest |cos(x)| ?
10-07-2021, 12:14 PM
Post: #1
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
smallest |cos(x)| ?
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 ?
10-07-2021, 03:40 PM (This post was last modified: 10-08-2021 06:57 AM by Werner.)
Post: #2
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: smallest |cos(x)| ?
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
10-07-2021, 06:40 PM (This post was last modified: 10-07-2021 10:10 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
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
10-08-2021, 07:02 AM (This post was last modified: 10-08-2021 12:39 PM by Werner.)
Post: #4
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: smallest |cos(x)| ?
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
10-08-2021, 08:24 AM (This post was last modified: 10-08-2021 08:26 AM by EdS2.)
Post: #5
 EdS2 Senior Member Posts: 595 Joined: Apr 2014
RE: smallest |cos(x)| ?
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!)
10-08-2021, 02:16 PM
Post: #6
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: smallest |cos(x)| ?
While I've been able to speed it up a bit, a thought occurred to me:
(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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
10-08-2021, 07:54 PM
Post: #7
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
(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)
10-08-2021, 08:17 PM (This post was last modified: 10-12-2021 11:36 AM by Albert Chan.)
Post: #8
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
(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
10-09-2021, 01:15 PM
Post: #9
 ijabbott Senior Member Posts: 1,297 Joined: Jul 2015
RE: smallest |cos(x)| ?
(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?

— Ian Abbott
10-09-2021, 02:09 PM
Post: #10
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
(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...
10-09-2021, 05:01 PM (This post was last modified: 10-24-2021 12:49 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
(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
10-10-2021, 01:36 PM
Post: #12
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
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.
10-11-2021, 01:16 PM
Post: #13
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: smallest |cos(x)| ?
(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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
10-11-2021, 04:22 PM (This post was last modified: 10-13-2021 12:52 AM by Albert Chan.)
Post: #14
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
(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.
10-11-2021, 10:05 PM
Post: #15
 Albert Chan Senior Member Posts: 2,678 Joined: Jul 2018
RE: smallest |cos(x)| ?
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)