[GAP Forum] Comparing real numbers

Leonard Soicher L.H.Soicher at qmul.ac.uk
Thu Sep 20 11:37:48 BST 2012


Dear GAP-Forum,

The following (not user-documented) DESIGN package service function, or 
parts of its code, may be of some use. Sturm sequences of polynomials
are used. 

   DESIGN_IntervalForLeastRealZero( f, a, b, eps)

Suppose that  f  is a univariate polynomial over the rationals, 
a,b are rational numbers with a<=b, and eps is a positive rational. 
If  f  has no real zero in the closed interval  [a,b],  then this
function returns the empty list  [].  Otherwise, let  alpha  be
the least real zero of  f  such that  a<=alpha<=b.  Then this function
returns a list  [c,d]  of rational numbers, with  c<=d,  d-c<=eps, 
and  c<=alpha<=d.  Moreover,  c=d  if and only if  alpha  is rational
(in which case  alpha=c=d).

Regards,
Leonard

On Thu, Sep 20, 2012 at 11:21:04AM +0200, Marek Mitros wrote:
> Hi again,
> 
> Since nobody answers I have created following function for comparing
> real numbers from field CF(5). It works for numbers tau,
> tau+1/10^4000. It does not for tau, tau + 1/10^5000 (tau is defined
> below).
> 
> Regards,
> Marek
> 
> 
> # Function for comparing real numbers from CF(5) field
> # The idea is to use fibonacci sequence approximating irrational
> number tau= (Sqrt(5)-1)/2 = 0,61803398874989484820458683436564
> # Define basis [1,tau] of the Field(tau). Any number x0 from this
> field can be presented as a*1+b*tau
> # for rational numbers a,b. Using approximation f(n)/f(n+1) for tau we
> obtain approximation
> # x(n)=a+b*f(n)/f(n+1) for x0.
> # Case 1: compare irrational number x0 with rational number p0.
> # In this case we look for n such that both x(n) and x(n+1) are
> greater than p0 or both smaller than p0.
> # Since all the time irrational x is between x(n) and x(n+1) we know
> that x is smaller or greater than p0.
> # Case 2: compare irrational numbers x0 with y0.
> # Look for n such that we have both [x(n), x(n+1)] are smaller or
> bigger than both [y(n), y(n+1)];
> 
> # function returns
> tau:=(Sqrt(5)-1)/2;
> F:=Field(tau);
> bas:=Basis(F, [1,tau]);
> 
> # At first try do not look at efficiency; just obtain good result;
> returns true when x<y
> comp_real5:=function(x,y)
>   local cc,a,b,c,d,n, maxn, nstep;
>   maxn:=10000;   nstep:=100;
>   if not ([x,y] in F^2) then Print ("Numbers outside the tau field,
> cannot compare\n"); return fail; fi;
>   if x=y then return false;
>   elif IsRat(x) and IsRat(y) then return x<y; # both rationals
>   elif IsRat(y) then
>      cc:=Coefficients(bas, x); a:=cc[1]; b:=cc[2];
>      n:=First(nstep*[1..maxn/nstep],
> n->((a+b*Fibonacci(n)/Fibonacci(n+1)<y) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<y))
>          or ((a+b*Fibonacci(n)/Fibonacci(n+1)>y) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)>y)) );
>      if n= fail then Print("Cannot compare using maxn
> approximation\n"); return fail; fi;
>      return (a+b*Fibonacci(n)/Fibonacci(n+1)<y) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<y);
>   elif IsRat(x) then
>      cc:=Coefficients(bas, y); a:=cc[1]; b:=cc[2];
>      n:=First(nstep*[1..maxn/nstep],
> n->((a+b*Fibonacci(n)/Fibonacci(n+1)>x) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)>x))
>          or ((a+b*Fibonacci(n)/Fibonacci(n+1)<x) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<x)) );
>      if n= fail then Print("Cannot compare using maxn
> approximation\n"); return fail; fi;
>      return (a+b*Fibonacci(n)/Fibonacci(n+1)>x) and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)>x);
>   else # both irrational
>      cc:=Coefficients(bas, x); a:=cc[1]; b:=cc[2];
>      cc:=Coefficients(bas, y); c:=cc[1]; d:=cc[2];
>      n:=First(nstep*[1..maxn/nstep],
> n->((a+b*Fibonacci(n)/Fibonacci(n+1)<c+d*Fibonacci(n)/Fibonacci(n+1))
>                          and
> (a+b*Fibonacci(n)/Fibonacci(n+1)<c+d*Fibonacci(n+1)/Fibonacci(n+2))
>                          and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<c+d*Fibonacci(n)/Fibonacci(n+1))
>                          and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<c+d*Fibonacci(n+1)/Fibonacci(n+2)))
>          or    ((a+b*Fibonacci(n)/Fibonacci(n+1)>c+d*Fibonacci(n)/Fibonacci(n+1))
>             and
> (a+b*Fibonacci(n)/Fibonacci(n+1)>c+d*Fibonacci(n+1)/Fibonacci(n+2))
>             and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)>c+d*Fibonacci(n)/Fibonacci(n+1))
>             and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)>c+d*Fibonacci(n+1)/Fibonacci(n+2)))
> );\
>      if n= fail then Print("Cannot compare using maxn
> approximation\n"); return fail; fi;
>      return ((a+b*Fibonacci(n)/Fibonacci(n+1)<c+d*Fibonacci(n)/Fibonacci(n+1))
>                          and
> (a+b*Fibonacci(n)/Fibonacci(n+1)<c+d*Fibonacci(n+1)/Fibonacci(n+2))
>                          and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<c+d*Fibonacci(n)/Fibonacci(n+1))
>                          and
> (a+b*Fibonacci(n+1)/Fibonacci(n+2)<c+d*Fibonacci(n+1)/Fibonacci(n+2)));
>    fi;
> end;
> 
> 
> 
> 
> On 9/19/12, Marek Mitros <marek at mitros.org> wrote:
> > Hi,
> >
> > I wonder how I could compare real numbers in GAP i.e. cyclotomics
> > which are real. I try Float but it doesn't work for cyclotomics - any
> > advice ? Say my numbers are limited to CF(5) field. Here is example:
> > sigma:=(E(5)+E(5)^4)/2; # = RealPart(E(5)) = (Sqrt(5)-1)/4 = 0.309...
> >  is real number. How I can compare whether it is bigger or smaller
> > that 1/2 ? It is known that this number is equal to Another real
> > number in CF(5) is
> > tau:=-RealPart(E(5)^2); # =-(E(5)^2+E(5)^3)/2=(Sqrt(5)+1)/4 =
> > 0.809...=sigma+1/2
> >
> > How could I develop function comparing elements of the field generated
> > by sigma and tau ? I mean comparison as real numbers. I can obtain
> > coefficients for given number x:
> > Coefficients(Basis(Field(sigma)), x)
> > Using this can I find the position of the number on real axis ? I.e.
> > compare with another field element or rational number ?
> >
> > I would really be nice to have possibility of comparing real subset of
> > cyclotomics are real numbers in GAP.
> >
> > Regards,
> > Marek
> >
> 
> _______________________________________________
> Forum mailing list
> Forum at mail.gap-system.org
> http://mail.gap-system.org/mailman/listinfo/forum



More information about the Forum mailing list