Phrases/Arith

From J Wiki
Jump to navigation Jump to search

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
)