[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

616.0. "Solving n = a*k + b" by CLT::GILBERT (eager like a child) Tue Nov 25 1986 02:40

I'd like to find near-solutions to the equation:  n = a*k + b,
where a and b are (given) real values, and the integers n and k
are desired.

Let {x} indicate the fractional part of x; then a 'near-solution'
has { a*k + b } < { a*j + b }, for all 0 <= j < k.  Surely there
is a better way to find near-solutions than to try successive k;
can someone provide an algorithm?
T.RTitleUserPersonal
Name
DateLines
616.1Something to tryMODEL::YARBROUGHTue Nov 25 1986 12:3614
>I'd like to find near-solutions to the equation:  n = a*k + b,
>where a and b are (given) real values, and the integers n and k
>are desired.

You might try approximating a and b with rationals, a1/a2, b1/b2; so

	n = a1/a2*k + b1/b2, or

	n*a2*b2 = a1*b1*k + a2*b1, or

	n*c1   =  k*c2 + c3

where the coefficients are now integers. This might be solvable with 
diophantine methods. If not, try different approximations for a, b.
616.2OoopsMODEL::YARBROUGHTue Nov 25 1986 12:385
Sorry, that should be

	n*a2*b2 = a1*b2*k + a2*b1

- Lynn
616.3BEING::POSTPISCHILAlways mount a scratch monkey.Tue Nov 25 1986 13:124
    This is somewhat related to topic 473.
    
    
    				-- edp
616.4CLT::GILBERTeager like a childFri Nov 28 1986 18:12128
    The following Pascal program contains a routine (solvex) that should
    do what's wanted.  See also note 597 for the application that this
    program solves.

program sum_of_sines (input, output);
const
	max_n = 100;
type
	r1 = single;
	r2 = double;
	r4 = quadruple;

function f4(x:integer): r4; begin f4 := quad(x); end;
function tan4(x:r4): r4; begin tan4 := sin(x)/cos(x); end;
function cot4(x:r4): r4; begin cot4 := cos(x)/sin(x); end;

[external(mth$hint)]
function trunc (x : r4) : r4; external;
 
[external(mth$hacos)]
function arccos (x : r4) : r4; external;
 
[external(mth$hnint)]
function round (x : r4) : r4; external;
 
function t4(x:r4): r4;
    begin
    if x < 0 then x := x - trunc(x) + f4(1);
    t4 := x - trunc (x);
    end;

procedure solvex (xa,xb: r4);
    { }
    {	Find an approximate integer solution to: n = a*k + b
    { }
    var
	za,zb:	r4;
	one:	r4;
	z0,z1,z2:	r4;
	k0,k1,k2:	r4;
	x,t0,t1:	r4;
	i:	integer;
	z4,kcos4,ksin4,khalf4:	r4;
	n4:	r4;
    begin
    one := 1;
    khalf4 := f4(1)/f4(2);
    kcos4 := cos (khalf4);
    ksin4 := sin (khalf4);

    za := t4(xa);
    zb := t4(xb);

    z1 := t4(za);		k1 := 1;

    z0 := t4(zb);		k0 := 0;
    while t4(z0+z1) > z0 do
	begin
	z0 := t4(z0+z1);	k0 := k0 + k1;
	end;
    z0 := t4(z0+z1);		k0 := k0 + k1;
{    writeln ('z0', k0:14:0, z0);}

    t0 := trunc ( one/z1 );
    z2 := t4(t0*za) - one;	k2 := t0 * k1;
{    writeln ('z2', k2:14:0, z2);}

    for i := 0 to 200 do
	begin

	t0 := trunc (z0 / (-z2));
	z0 := t4(z0+t0*z2);	k0 := k0 + t0 * k2;
{	writeln ('z0', k0:14:0, z0);}

	t1 := trunc (z1 / (-z2));
	z1 := t4(z1+t1*z2);	k1 := k1 + t1 * k2;
{	writeln ('z1', k1:14:0, z1);}

	z2 := z1 + z2;		k2 := k1 + k2;
{	writeln ('z2', k2:14:0, z2);}

	if t0 <> 0
	then
	    begin
	    x := xa * k0 + xb;
	    writeln ('**      ', x, x - round(x) :14);

	    n4 := trunc (x);
	    z4 := (kcos4 - cos(n4+khalf4)) / 2 / ksin4;
	    writeln ('**', n4:20:0, z4, z4-round(z4):14);
	    end;

	end;

    end;

procedure init;
    var
	x4,y4,z4,kcos4,ksin4,khalf4:	r4;
	pi:	r4;
	i,j:	integer;
    	n,k:	integer;
    begin
{ }
{   We want to find n such that cos(n+1/2) = cos(1/2) - 2*sin(1/2)
{   In other words, n = arccos(cos(1/2)-2*sin(1/2)) - 1/2 + k*pi
{ }
    pi := 4*arctan(f4(1));
    khalf4 := f4(1)/f4(2);
    kcos4 := cos (khalf4);
    ksin4 := sin (khalf4);
    z4 := arccos (kcos4-2*ksin4) - khalf4;

    writeln (pi, z4);
    writeln (cos( z4+khalf4), kcos4-2*ksin4);
    writeln (cos(-z4-khalf4), kcos4-2*ksin4);

    solvex (2*pi, z4);
    solvex (2*pi, -1-z4);
{
    z4 := (kcos4 - cos(f4(n)+khalf4)) / 2 / ksin4;
    writeln (n,z4);
}
    end;

begin
init;
end.