| >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.
|
| 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.
|