Function Repository Resource:

SimultaneousRationalize

Source Notebook

Approximately rationalize a set of numbers to give a common denominator

Contributed by: Daniel Lichtblau

ResourceFunction["SimultaneousRationalize"][{v1,v2,}]

finds rationals r1,r2, approximating v1,v2, such that the ri share a small common denominator.

ResourceFunction["SimultaneousRationalize"][{v1,v2,},weight]

uses weight to influence the size of the resulting denominator.

Details and Options

ResourceFunction["SimultaneousRationalize"] uses lattice reduction to obtain a common denominator.
Larger values of weight tend to give closer approximations, albeit with larger denominators.
The default value for weight is 20.
ResourceFunction["SimultaneousRationalize"] can take a set of rationals as input and produce an approximating set with a smaller common denominator.
ResourceFunction["SimultaneousRationalize"] is similar to the legacy package function NumberTheory`Rationalize`AffineRationalize.

Examples

Basic Examples (1) 

Find rational approximations to a set of approximate real values:

In[1]:=
data = {5.831666666666677`, 22.96202080402009`, 40.0120210909091`, 45.92557465116284`};
ResourceFunction["SimultaneousRationalize"][data]
Out[1]=

Scope (7) 

Start with a set of approximate real values:

In[2]:=
data = {5.831666666666677`, 22.96202080402009`, 40.0120210909091`, 45.92557465116284`};

Use the weight argument to find rational approximations with a small common denominator:

In[3]:=
ResourceFunction["SimultaneousRationalize"][data, 10]
Out[3]=

With a larger weight we obtain closer approximations:

In[4]:=
ResourceFunction["SimultaneousRationalize"][data, 100]
Out[4]=

With a sufficiently small weight we find integer approximations:

In[5]:=
ResourceFunction["SimultaneousRationalize"][data, 5]
Out[5]=

SimultaneousRationalize can work with rationals in the input:

In[6]:=
ratdata = Rationalize[data, 10^(-5)]
Out[6]=

Find rational approximations with a smaller common denominator:

In[7]:=
ResourceFunction["SimultaneousRationalize"][ratdata]
Out[7]=

The values in question are approximately commensurate, i.e. all are multiples of a common value. In this case the smallest element in the list:

In[8]:=
ResourceFunction["SimultaneousRationalize"][Rest@data/First@data, 10]
Out[8]=

Possible Issues (2) 

Values that are small in magnitude compared to unity will often be approximated by zero:

In[9]:=
data = {.003245, 5.831666666666677`, 22.96202080402009`, 40.0120210909091`, 45.92557465116284`};
ResourceFunction["SimultaneousRationalize"][data]
Out[9]=

Use a weight commensurate with the reciprocal of the smallest value to force the approximation to be nonzero:

In[10]:=
ResourceFunction["SimultaneousRationalize"][data, 200]
Out[10]=

Neat Examples (9) 

Simultaneous rationalization can be used on light curve data to estimate the period of a variable star. We will take for our example Kepler project data for the Cepheid star 2308. The reference period at the OGLE web site is given as 5.7373087 days. Data and more can also be found at the Wolfram Demonstration "CepheidVariableStarLightCurve" by Jeffrey Bryant. This example is also found in the resource function IrregularPeriodogram.

Import the data from OGLE:

In[11]:=
cep2308Idata = Import[
    "http://ogledb.astrouw.edu.pl/~ogle/CVS/data/I/08/OGLE-LMC-CEP-2308.dat"
    , "Table"][[All, 1 ;; 2]];

Separate into observation times and light intensities, and subtract the average from the set of values:

In[12]:=
{times2308, ovalues2308} = Transpose[cep2308Idata];
times2308 = times2308 - times2308[[1]];
mean2308 = Mean[ovalues2308];
values2308 = ovalues2308 - mean2308;
ListPlot[Transpose@{times2308, values2308}]
Out[12]=

We extract candidate period multiples as follows. The idea is that "common" time differences for approximately equal values correspond to times that are separated by a multiple of the basic period. (1) Compute an ϵ for proximity (a small fraction of the range of values). (2) Find all time differences between pairs that have their respective values within this ϵ. (3) Separate time differences into bins according to whether they are in proximity. This step uses an input threshold for proximity, with default of 0.1. (4) Cull out those time differences that do not appear often. This step uses an input minimum for bin sizes, with default of 15. (5) Take means of each retained bin. These are the candidate approximations for multiples of periodicities.

Utility function to find candidate periodicity multiples:

In[13]:=
periodMultiples[vals_, times_, minrunlen_ : 15, frac_ : 300, thresh_ : .1] := Module[
  {nfunc, eps, nbrs, nbrsSplit, xdiffs, splitdiffs},
  nfunc = Nearest[vals -> Range[Length[vals]]];
  eps = (Max[vals] - Min[vals])/frac;
  nbrs = Table[nfunc[valj, {Infinity, eps}], {valj, vals}];
  nbrsSplit = Map[{First[#], Rest[#]} &, nbrs];
  xdiffs = Sort[Flatten[
     Abs[Table[
       times[[nbhdj[[1]]]] - times[[nbrjk]], {nbhdj, nbrsSplit}, {nbrjk, nbhdj[[2]]}]]]];
  splitdiffs = Split[xdiffs, Abs[#2 - #1] < thresh &];
  Map[Mean, Select[splitdiffs, Length[#] >= minrunlen &]]
  ]

Compute candidate periodicity multiples for Cepheid 2308:

In[14]:=
periods = periodMultiples[values2308, times2308, 40]
Out[14]=

Remove the smallest value since it is too small to be a period and is most likely an artifact of the average sampling time gap:

In[15]:=
periods = Rest[periods]
Out[15]=

Use simultaneous Diophantine approximation to find integer ratios that approximate the period multiple ratios:

In[16]:=
multiples = ResourceFunction["SimultaneousRationalize"][
  Rest[periods]/First[periods], 5]
Out[16]=

Divide the approximated period multiples by the approximated multipliers (the initial multiplier is assumed to be unity):

In[17]:=
periodEstimates = periods/Prepend[multiples, 1]
Out[17]=

Take the average to get our period approximation:

In[18]:=
periodEstimate = Mean[periodEstimates]
Out[18]=

Fold the light curve by this estimated period:

In[19]:=
foldTimes[times_, period_] := period*FractionalPart[times/period]
In[20]:=
ListPlot[
 Transpose[{foldTimes[times2308, periodEstimate], values2308}]]
Out[20]=

The "fuzz" points are due to the estimate being off in the second decimal place, but still one can clearly discern the main wave pattern.

Version History

  • 1.0.0 – 01 May 2020

Related Resources

License Information