In the context of floating point operations, I am usually comfortable with 9!:18'' as epsilon.
For rational arithmetic, 2^_1024x isn't actually all that much of an issue. Consider: #":*:*:2^_1024x 1236 But when floating point operations are in the mix, they become the controlling issue (unless, of course, we introduce an adequate rational approximation...). Thanks, -- Raul On Wed, Jan 25, 2023 at 11:49 AM Ian Clark <[email protected]> wrote: > > @ Raul > While you were thinking, so was I, but along different lines. However I > think we came to the same conclusion. > > When I first saw your 2^_1024x -- my knee-jerk was: best not make it too > small. > > I haven't studied math/calculus, but when the "derivative" primitives were > removed from JE (d. D. D:) I went back to the finite-difference methods I'd > used in earlier APL implementations of math/cal. These needed a small > step-size, but not too small. I messed around (let's be honest) until I > found something that worked. Usually. > > My fiercely pragmatic reasoning went like this. If a novice user tries > playing with finite difference methods using math/cal, they ought to be > given a safe "epsilon", i.e. one that cancels out if used (once) in both > the numerator and denominator of some scientific calculation. > > Now only a fool will deliberately go and write: > epsilon * y % epsilon > But it might happen out-of-sight, in the sort of cascade of calculations > which turns up when playing with math/cal -- or TABULA, which uses it. > > Here's my test-rig, which I've just written (to preserve sanity, in a head > with too much going on in it): > > NB.21 - epsilon in math/calculus > > 0 : 0 > > Wed 25 Jan 2023 15:52:28 > > - > > Note these articles: > > https://en.wikipedia.org/wiki/Finite_difference > > https://en.wikipedia.org/wiki/Calculus_of_variations > > ) > > clear'' > > > case=. 1 > > > 3 : 0 case > > select. y > > case. 0 do. epsilon=: 2^_1020x > > case. 1 do. epsilon=: 2^_1021x > > case. 2 do. epsilon=: 2^_1022x > > case. 3 do. epsilon=: 2^_1023x > > case. 4 do. epsilon=: 2^_1024x > > end. > > ) > > > float=: 1.0 * ] > > > h=: hiddencalc=: {{ (%&epsilon) (*&epsilon) y }} > > h2=: hiddencalc2=: {{ (*&epsilon (%&epsilon y)) }} > > h3=: hiddencalc3=: {{ (epsilon * y % epsilon) }} > > h4=: hiddencalc4=: {{ (epsilon %~ y * epsilon) }} > > > smoutput 'epsilon' ; epsilon > > smoutput 'float epsilon' ; float epsilon > > > smoutput |.(float hiddencalc z) ; z=: 1e_2 > > smoutput |.(float hiddencalc2 z) ; z > > smoutput |.(float hiddencalc3 z) ; z > > smoutput |.(float hiddencalc4 z) ; z > > > The user is invited to play with choices of nouns: case and: z . > > > Now case=. 1 --is the highest I can go, and get z back again after applying > hiddencalc to it. > Which, I fancy, is your: 2^_1021x in place of: 2^_1024x > > Ian > > On Wed, 25 Jan 2023 at 15:01, Raul Miller <[email protected]> wrote: > > > Actually, after thinking about this, an epsilon at or near 2^_1021x > > (or 4e_308) would probably be a better default. Having an epsilon > > which can be safely represented as a floating point value is probably > > for the best, especially initially. > > > > Technically, 1e_318 might be used, but anything much smaller than > > 4e_308 can only be represented as a denormalized floating point > > number, which means casting to rational and then back to float turns > > it into a zero, and I haven't quite convinced myself that this is a > > bug. > > > > Thanks, > > > > -- > > Raul > > > > On Wed, Jan 25, 2023 at 8:45 AM Raul Miller <[email protected]> wrote: > > > > > > I haven't looked at math/calculus recently. > > > > > > However, it sounds like it could use a few epsilons... > > > > > > Perhaps a default of 2^_1024x would be appropriate? > > > > > > Thanks, > > > > > > -- > > > Raul > > > > > > On Wed, Jan 25, 2023 at 6:17 AM Ian Clark <[email protected]> wrote: > > > > > > > > Thanks Raul. > > > > > > > > This saves me locating the source of math/tabula and friends, to > > update it. > > > > It's several years now since I touched the code, and I've forgotten > > how to > > > > do it -- and it's possibly changed. > > > > > > > > But I ought to be grateful for the current j904 breaking the code, > > because > > > > it has alerted me to latent bugs regarding the use of rational numbers > > by > > > > math/cal and math/uu. For me this is the tip of a murky iceberg. > > > > > > > > It's not so much that I want to retain a sensible looking constant for > > > > rational-infinity, as the fact that infinities (or exploding large > > finite > > > > number representations in general) can arise in so many ways, and > > there's > > > > no guarantee they will equate with whatever I settle on as a > > "reference" > > > > rational-infinity. In particular, math/cal's use of Newton's method > > with > > > > rational numbers is particularly fraught, with an emerging host of > > spooky > > > > reasons for non-convergence I don't have words to describe (some play > > on > > > > the terms Moiré, resonance, Nyquist … might be needed). This is the > > > > flagship feature of math/cal we're talking about. > > > > > > > > Hitherto math/cal has played whack-a-mole with issues as they arose. I > > fear > > > > that the move to new arithmetic routines will deliver a fresh load of > > moles > > > > to whack. And just when my attention is diverted elsewhere. > > > > > > > > Ian > > > > > > > > On Wed, 25 Jan 2023 at 04:36, Raul Miller <[email protected]> > > wrote: > > > > > > > > > I've noticed an odd quirk here. > > > > > > > > > > 1.2r3.4 > > > > > 0.352941 > > > > > _r3.4 > > > > > |ill-formed number > > > > > > > > > > This issue is present in j903. > > > > > > > > > > I have opted to retain this quirk for j904, because it doesn't seem > > to > > > > > be important and it makes the implementation a bit simpler. > > > > > > > > > > (Also, other than this quirk, _r1 and friends will work in the next > > > > > update to j904.) > > > > > > > > > > Thanks, > > > > > > > > > > -- > > > > > Raul > > > > > > > > > > On Tue, Jan 24, 2023 at 10:38 PM Raul Miller <[email protected]> > > > > > wrote: > > > > > > > > > > > > Right... Aside from adding libgmp support, a change from j903 is > > that > > > > > > j903 had an extended precision infinity which was used in parsing > > > > > > numeric constants, but j904 does not. > > > > > > > > > > > > And, when I was rewriting the bit that handles rational constants, > > I > > > > > > overlooked some of the ways of representing rational infinity. > > > > > > > > > > > > I'm testing a fix for this problem right now. It should be ready > > soon. > > > > > > > > > > > > Thanks, > > > > > > > > > > > > -- > > > > > > Raul > > > > > > > > > > > > On Tue, Jan 24, 2023 at 10:18 PM Henry Rich <[email protected]> > > > > > wrote: > > > > > > > > > > > > > > Decision, decisions. How /should/ you specify an extended > > infinity? > > > > > I say > > > > > > > > > > > > > > _x > > > > > > > |ill-formed number > > > > > > > > > > > > > > There could be alternatives. _r(any finite) and (any non0)r0 are > > both > > > > > > > reasonable. > > > > > > > > > > > > > > NOTE that the GMP library that we have moved to has no way to > > represent > > > > > > > extended infinity. Raul has chosen 1r0 as our internal > > representation > > > > > > > of extended infinity, so infinity will always have rational > > precision, > > > > > > > not extended integer. > > > > > > > > > > > > > > For display, we get it right: > > > > > > > > > > > > > > 1r0 > > > > > > > _ > > > > > > > > > > > > > > It seems that _r(any) should be converted to infinity - and _x > > and __x > > > > > > > too I think. This is in Raul's area. > > > > > > > > > > > > > > If you have a dependency on the internal representation of > > infinity it > > > > > > > will be on you to update it. > > > > > > > > > > > > > > Henry Rich > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > On 1/24/2023 9:00 PM, Ian Clark wrote: > > > > > > > > j903 accepts -- but j904 rejects -- this way of defining > > rational > > > > > infinity: > > > > > > > > > > > > > > > > _r1 > > > > > > > > > > > > > > > > |ill-formed number > > > > > > > > > > > > > > > > | _r1 > > > > > > > > > > > > > > > > | ^ > > > > > > > > > > > > > > > > JVERSION > > > > > > > > > > > > > > > > Engine: j904/j64arm/darwin > > > > > > > > > > > > > > > > Beta-k: commercial/2023-01-24T04:42:28 > > > > > > > > > > > > > > > > Library: 9.04.10 > > > > > > > > > > > > > > > > Qt IDE: 2.0.3/6.2.4(6.2.4) > > > > > > > > > > > > > > > > Platform: Darwin 64 > > > > > > > > > > > > > > > > Installer: J904 install > > > > > > > > > > > > > > > > InstallPath: /applications/j904 > > > > > > > > > > > > > > > > Contact: www.jsoftware.com > > > > > > > > > > > > > > > > A workaround is to use 1r0 instead: > > > > > > > > > > > > > > > > > > > > > > > > 1r0 > > > > > > > > > > > > > > > > _ > > > > > > > > > > > > > > > > datatype 1r0 > > > > > > > > > > > > > > > > rational > > > > > > > > > > > > > > > > > > > > > > > > Not a lot of j-ers willl have a use for rational [minus] > > infinity, > > > > > but IMO > > > > > > > > a beginner might find it more intuitive to define it as _r1 > > rather > > > > > than 1r0 > > > > > > > > . Maybe it's no big deal in itself, but it breaks 3 addons, > > viz: > > > > > math/cal, > > > > > > > > math/uu -- and in consequence math/tabula: > > > > > > > > > > > > > > > > > > > > > > > > load'math/uu' NB. Launch UU only > > > > > > > > |ill-formed number in script, executing monad 0!: 0 > > > > > > > > |any word beginning with a digit or _ must be a valid number > > > > > > > > | BADRAT=: __r1 > > > > > > > > | ^ > > > > > > > > |[-33] /applications/j904/addons/math/uu/uu.ijs > > > > > > > > > > > > > > > > load'math/cal' > > > > > > > > |ill-formed number in script, executing monad 0!: 0 > > > > > > > > |any word beginning with a digit or _ must be a valid number > > > > > > > > | BAD_EXE_VALUE=: __r1 > > > > > > > > | ^ > > > > > > > > |[-92] /applications/j904/addons/math/cal/cal.ijs > > > > > > > > > > > > > > > > > > > > > > > > So I'd call it a bug. > > > > > > > > > > > > > > > > > > > > > > > > There's been plenty of time to discover this. Sorry it's taken > > me so > > > > > long. > > > > > > > > > > > > > > > > > > > > > > > > Ian Clark > > > > > > > > > > > > > > > ---------------------------------------------------------------------- > > > > > > > > For information about J forums see > > > > > http://www.jsoftware.com/forums.htm > > > > > > > > > > > > > > > > ---------------------------------------------------------------------- > > > > > > > For information about J forums see > > http://www.jsoftware.com/forums.htm > > > > > > > ---------------------------------------------------------------------- > > > > > For information about J forums see > > http://www.jsoftware.com/forums.htm > > > > > > > > > ---------------------------------------------------------------------- > > > > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
