@ 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

Reply via email to