smallest cos(x) ?

10072021, 12:14 PM
Post: #1




smallest cos(x) ?
Last month, I had proposed a better complex tan(z) for Free42
We considered HP71B mathrom implementation. Quote:tan((x,y)) = ( sin(x)*cos(x)/d , sinh(y)*cosh(y)/d ) 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 * 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.1898063658835771e19')) Others had search for this too. From Accurate trigonometric functions for large arguments >>> eps(mpfr('0x1.6ac5b262ca1ffp+849')) (mpfr('5.3193726483265414e+255'), mpfr('4.6871659242546277e19'))  How to setup a Free42 program to search for this ? 

10072021, 03:40 PM
(This post was last modified: 10082021 06:57 AM by Werner.)
Post: #2




RE: smallest cos(x) ?
Quick hack (corrected):
Code: 00 { 145Byte Prgm } Of course, MANT has been coded long ago, thanks to Dieter: Code: 00 { 25Byte Prgm } Then, the loop, for k=1..10000: Code: 00 { 84Byte Prgm } Which yields the smallest value for k=953, being 8.2E35, for x=2995.508595197867852874130465957006 If I haven't made an error, that is ;) Cheers, Werner 

10072021, 06:40 PM
(This post was last modified: 10072021 10:10 PM by Albert Chan.)
Post: #3




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 = 100digitsofpi %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.200102230416340029960966876550294E35 c:\> eps 82 1.532138639232766081410107492878833E35 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.200102230416340029960966876550294E35 pi 82425623.5x ?68 d =34 ? 2.5894773325515822091636756232496249999999998467861360767233918589893E+8 2.589477332551582209163675623249625E+8  abs ?34 1.532138639232766081410107492878833E35 cos(x) = sin(ε) = ε  ε^3/3! + ... ≈ ε Free42 Decimal: 2995.508595197867852874130465957006 COS → 8.200102230416340029960966876550293e35 258947733.2551582209163675623249625 COS → 1.532138639232766081410107492878833e35 Intel Decimal128 angle reduction code is very good 

10082021, 07:02 AM
(This post was last modified: 10082021 12:39 PM by Werner.)
Post: #4




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 

10082021, 08:24 AM
(This post was last modified: 10082021 08:26 AM by EdS2.)
Post: #5




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!) 

10082021, 02:16 PM
Post: #6




RE: smallest cos(x) ?
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 

10082021, 07:54 PM
Post: #7




RE: smallest cos(x) ?
(10082021 02:16 PM)Werner Wrote: instead of doing (k+0.5)*pi, which has 2 rounding errors (one for pi, one for *), do: "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.686E35 If x has cos(x) minimized, x scaled to 34 digits, is very close to integer. Say, combined relative error = 3.7E35 (as long as error < 5E35, we are safe) ulp_error(x) < 10^34 ULP * 3.7E35 = 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) 

10082021, 08:17 PM
(This post was last modified: 10122021 11:36 AM by Albert Chan.)
Post: #8




RE: smallest cos(x) ?
(10082021 07:02 AM)Werner Wrote: A bit counterintuitive that the smallest COS happens for such a large argument value; This give me an idea ! For IP(x) of 9 digits, calculate (pi/2*10^(349)), 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.532138639232766081410107492878833e35 247276870.5 PI * COS → 4.596415917698298244230322478636498e35 Best semiconvergent is not as good: 2k+1 = 18404082 + 164851247*3 = 512957823 256478911.5 PI * COS → 5.589328359211609676578079200248261e34 

10092021, 01:15 PM
Post: #9




RE: smallest cos(x) ?
(10082021 07:02 AM)Werner Wrote: Running my programs to 100 million would take about 45 minutes ;) There should be infinitely many such values for x producing a smaller ε, shouldn't there? — Ian Abbott 

10092021, 02:09 PM
Post: #10




RE: smallest cos(x) ?
(10092021 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 ...) (10072021 06:40 PM)Albert Chan Wrote: Free42 Decimal: Using convergents / semiconvergents idea, for x < 1E50, these two have smaller ε: 9.510911580848780648392560778912209e40 COS → 1.328353856527533918552271223903534e35 9.341730789500356812974471132146251e46 COS → 1.028415848209791685669808767152452e35 To show how hard to get similar sized ε, this is true x, before rounded. 30274171828040719778093696476072003369775.5 * pi = 95109115808487806483925607789122090000000.000000000000000000000000000000000013283538... 29735652643654715492245370225672550765119900844.5 * pi = 93417307895003568129744711321462509999999999999.999999999999999999999999999999999989... 

10092021, 05:01 PM
(This post was last modified: 10242021 12:49 PM by Albert Chan.)
Post: #11




RE: smallest cos(x) ?
(10092021 02:09 PM)Albert Chan Wrote: Both effect tends to cancel out, giving min(ε) not much smaller than machine epsilon. 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.028415848209791685669808767152452e35 Update: Brute force proof, for Free42Decimal, cos(x) > 10^37 see https://www.hpmuseum.org/forum/thread17...#pid153539 

10102021, 01:36 PM
Post: #12




RE: smallest cos(x) ?
Trivia: 2 convergents cannot have both even denominator (or both even numerator)
Code: def genConvergents(coefs): 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. 

10112021, 01:16 PM
Post: #13




RE: smallest cos(x) ?
(10082021 07:54 PM)Albert Chan Wrote:(10082021 02:16 PM)Werner Wrote: instead of doing (k+0.5)*pi, which has 2 rounding errors (one for pi, one for *), do: 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 

10112021, 04:22 PM
(This post was last modified: 10132021 12:52 AM by Albert Chan.)
Post: #14




RE: smallest cos(x) ?
(10112021 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: 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 ≡ 1E34) 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.469009289781287037341230607307734e35  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 ≈ (10.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 overestimated, relerr(PI) = −0.369 ULP 180/PI underestimated, 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. 

10112021, 10:05 PM
Post: #15




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 34dgits, relerr ≈ (10.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 34dgits, relerr ≈ 0.0004941/0.7874 < 0.00063 ULP 

« Next Oldest  Next Newest »

User(s) browsing this thread: