Ok, I see from github that you're working on a Float125 and Float127 implementation. Why not Float128?, and why not use Julia BigFloats? Out of curiosity, I did a few tests on Julia's sin(BigFloat).
julia> p=Float64(pi) 3.141592653589793 julia> length(string(p)) 17 julia> sin(p) 1.2246467991473532e-16 Here the error is O(eps_machine), as it should be, since sin(x) is well-conditioned. But for BigFloat julia> p=BigFloat(pi) 3.141592653589793238462643383279502884197169399375105820974944592307816406286198 julia> length(string(p)) 80 julia> sin(p) 1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77 That's two orders of magnitude larger than machine epsilon 1e-79 (subtracting 1 from 80 because of the period). Is this what's troubling you? And a side note: unless I'm missing something, the Julia docs on BigFloat seem misleading or ambiguous. BigInt and BigFloat are discussed together as arbitrary precision arithmetic, though if BigInt just wraps GNU MFPR, it's large but finite precision. Right? John On Saturday, October 17, 2015 at 9:39:54 AM UTC-4, Jeffrey Sarnoff wrote: > > thanks, good pointer. > > I am working on routines for a double-double-like floating point type. > Straight double-double implementations have very nice arithmetic > properties; in my experimentation, most double-double trig routines resolve > fewer bits than I want. > > I want to take as much advantage of the built-in elementary functions as > possible -- the outer representation is a non-overlapping pair of Float64s > where the available value is their implicit extended precision sum: hipart > + lopart. It would have been spectacular to use sin(x + dx) = > sin(x)*cos(dx) + cos(x)*sin(dx) > and find sin(hipart + lopart) using only built-in trig and the module's > extended precision arithmetic: > sin(hipart)*cos(lopart)+cos(hipart)*sin(lopart). > Unfortunately, that requires much more than double precision to work well. > > > I have trig working at ~100 sigbits for angles in -2pi..2pi, but too > slowly. > > > On Saturday, October 17, 2015 at 9:12:38 AM UTC-4, John Gibson wrote: > >> Search for the comment that begins "OK kiddies, time for the pros...." in >> http://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions >> >> John >> >> On Saturday, October 17, 2015 at 9:00:35 AM UTC-4, John Gibson wrote: >>> >>> Why are you trying to roll your own sin(x) function? I think you will be >>> hard pressed to improve on the library sin(x) in either speed or accuracy. >>> >>> John >>> >>> On Saturday, October 17, 2015 at 3:38:17 AM UTC-4, Jeffrey Sarnoff wrote: >>>> >>>> I had tried to find a clean way to jump into the taylor series using >>>> the well approximated sin(x) or cos(x) and so additionally limit the terms >>>> used -- there may be / probably is a way in concert with an additional >>>> tabulation (which would be fine in this case). Taylor's theorem is not >>>> numerically crisp, so while I can identify the next term (using, perhaps >>>> eps(sin(x))/3) I don't know how to back out the delta between accumulation >>>> of the series to the prior term and the value sin(x::Float64). >>>> >>>> On Saturday, October 17, 2015 at 3:29:28 AM UTC-4, Jeffrey Sarnoff >>>> wrote: >>>>> >>>>> That has promise, Kristoffer. I did port something of that nature, >>>>> expecting it to work well -- but there was some numerical mush in more >>>>> than >>>>> a couple of trailing bits in some cases. >>>>> Using more terms did not help. Thinking about it just now, it might be >>>>> more robustly stable if I expand in one direction only upfrom or >>>>> downfrom >>>>> some pretabulated points. >>>>> >>>>> On Saturday, October 17, 2015 at 2:47:57 AM UTC-4, Kristoffer Carlsson >>>>> wrote: >>>>>> >>>>>> Use a truncated Taylor series around the point maybe? >>>>> >>>>>
