Phrases/Arith
Basic operations
In J504 a root of degree higher than 2 (or, rather <.@%:) was not implemented correctly for extended argument (corrected in J601). The workaround was to use
introot=: 4 : '(] [`]@.(y&>:@^&x@]) +)/ 1x,|.*/\.2x#~x<.@%~2x<.@^.y'
(see also Essays/IntRoot)
prp=: [: *./ ([: *./ 1: = {: +. }:)\
checks if integers from a given list are pairwise relatively prime
Interpolation
Verbs for linear and piecewise-linear interpolation
NB. x is 2 $ interval y is value(s) (any rank)
NB. Result is fractional position(s) within interval
invlerp =: ( ({.@[ - ]) % -/@[ )"1 _ :. lerp
NB. x is 2 $ interval y is fractional position(s) within interval (any rank)
NB. Result is value(s)
lerp =: (+/@:* (,~ -.))"1 0"1 _ :. invlerp
NB. x is a list of points, in ascending order, y is value
NB. result is interval . fractional position for the interval bracketing y
NB. If y is out of bounds, use the endpoint
NB. After the test for too low/too high, we keep the vector of items of x less than y, and
NB. the first item of x ge y. The count of these, -2, is the integer part, and we invlerp
NB. within the last 2 to find the fraction
piecewiseinvlerp =: (( ( (] (-&2@#@] + (invlerp~ _2&{.)) >:@(i.&0)@:< {. [) ` (<:@#@[)) @. ((>: {:)~)) ` 0:) @. ((< {.)~) "1 0
NB. x is a list of vectors, representing data values, each vector at one point
NB. y is interval . fractional position within interval
NB. Result is value
piecewiselerp =: 13 : '((,~ -.) 1 | y) +/@:(*"_1) ((#x) | (, >:) <. y) { x'"_ 0
Different way that allows for higher order interpolation is described in this forum post
NB. Resampling
NB. y is (x values),:(y values)
NB. x is new x values
NB. Result is new y values
resamp =: 4 : 0
NB. Intervals are numbered by the index of the sample that
NB. ends the interval. So, interval 0 is before the first sample
NB. and interval {:$y is after the last. We calculate the
NB. interval number for each x and then, if it is one of those
NB. off-the-end intervals, adjust to the nearest interior interval.
NB. This means we extrapolate out-of-range values using the slope
NB. of the outermost intervals.
ix =. 1 >. (<:{:$y) <. (0{y) I. x
NB. Calculate the interpolating polynomial for each interval.
NB. Here we use linear interpolation, so the polynomial is (y value),(dy/dx)
NB. Create a polynomial for the first interval (off the beginning),
NB. using the slope of the first internal interval
intpoly =. (1 { y) ,. (,~ {.) %~/ 2 -/\"1 y
NB. The value to return is the interpolating polynomial, evaluated
NB. given the distance between the desired value and the origin point
NB. (i. e. right endpoint) of the interval
(ix { intpoly) p. x-(<0;ix) { y
)
Interpolation on a circle
When values needed to be interpolated are angular, say, degrees, then using these routines may cause incorrect result. For example, if we measure some angle 355° at 11:00, it increases with time and becomes 5° at 12 and we want to know what was the angle at 11:24, then using 'lerp' directly gives
355 5 lerp 0.4 215
What we need to do is to adjust values that crossed 360° boundary, perform interpolation and coerce result back to 0°..360° range.
NB. Assumes that y is vector of increasing values modulo x. For example NB. if x=360, then y may mean angular values expressed in degrees. NB. Whenever next number is smaller than the previous, this must mean that it NB. crossed x boundary modadjust=:] + [: +/\ [ * 0 , 2 >/\ ] circadjust=:360&modadjust
(360 circadjust 355 5) lerp 0.4 359
Fractions
Continued fractions
Continued fractions are fractions of a form a+1/(b+1/(c+...+1/z), where a,b,c are integers >0 and z is integer>1. Simplest way to represent a finite continued fraction in J is a list of integers.
NB.* ecf v evaluate continued fraction ecf=: (+%)/ NB.* cf v convert a number to continued fraction cf=: ]`(<. , $:@%@(1&|))@.(~: <.)
Continued fraction representation of an irrational number is an infinite list. Chopped at some point to a finite list it gives "good" (in certain sense "best") approximation of a number. The following phrase give rational approximations of π. First in this list is mentioned in the bible, second discovered by Archimedes and fourth is sometimes attributed to Tsu Ch'ung-Chi (Zu Chongzhi).
cf x:o.1 3 7 15 1 292 1 1 1 2 1 3 1 1 5 1 3 11 1 1 4 1 2 14 ecf\ cf x:o.1 3 22r7 333r106 355r113 103993r33102 104348r33215 ...
Another useful fact about continued fractions is that quadratic irrationalities have periodic continued fraction representation.
cf x:%:2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 45 2 1 2 11 1 3 6 cf x:%:3 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 5 7 2 6 1 2 1 4 2 NB. non-periodic tails here are due to inexact floating point representation
Converting a number (possibly obtained by numeric methods) to continued fraction may hint it is being a quadratic irrationality.
See also this essay.
Decimal and other positional fractions
require 'strings'
DIGITS=:'0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'
NB.* d2f v convert decimal representation of fraction into
NB. normal fraction. Accepts 'x.yyy(zzzz)' where (zzzz) -- period
d2f=: (10x & $:) : (4 : 0)
r=. '(' taketo y NB. regular part
p=. (#~ e.&DIGITS) '(' takeafter y NB. periodic part
p10=. # '.'takeafter r
s=. *0.5 - '-' e. r
r=. r-.'-'
i=. x #. DIGITS i. '.' taketo r
f=. x #. DIGITS i. '.' takeafter r
a=. x #. DIGITS i. p
s*i+(x^-p10) * f+a%<:x^#p
NB. p1=.(".@(}:&.((-p10)&|.)) '-' taketo&.|. '-',r)% x ^ p10
NB. p2=.(".@(,&'x') %&x: <:@(x.&^)@#) p
NB. s*p1 + p2 % x ^ p10
)
NB. convert rational into periodic fraction on given base
NB. optional left argument is base and (optional) maximum number of
NB. digits in fractionalpart to look for
NB. For example:
NB. f2d 1r7 NB. same as (10x;100) f2d 1r7
NB. (2;15) f2d 1r29
NB. (2;26) f2d 1r29
f2d=: formatrational=: (10x & $:) : (4 : 0)
'base maxdigits'=. 2{.(boxxopen x),<100
base=. x:base
assert. (base=<.base)*.(1<base)
r=. ''
if. 0>y do. r=. '-' [ y=. -y end.
r=. r,DIGITS{~base&#.^:_1<.y
y=. 1x|y
if. 0<y do. r=. r,'.' end.
rm=. ,y NB. list of remainders
while. 0<y do.
v=. y*base
y=. 1x|v
i=: rm i. y
if. i>:#rm do.
r=.r,DIGITS{~<.v
rm=. rm,y
else.
r=. r,DIGITS{~<.v
NB. cycle detected in periodic fraction
r=. r (}.~ , '('"_ , {.~ , ')'"_) i-#rm
break.
end.
if. 0>maxdigits=. maxdigits-1 do. r=. r,'...' break. end.
end.
r
)