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 )