48 lines
1.3 KiB
Text
48 lines
1.3 KiB
Text
sternBrocot: procedure options(main);
|
|
%replace MAX by 1200;
|
|
declare S(1:MAX) fixed;
|
|
|
|
/* find the first occurrence of N in S */
|
|
findFirst: procedure(n) returns(fixed);
|
|
declare (n, i) fixed;
|
|
do i=1 to MAX;
|
|
if S(i)=n then return(i);
|
|
end;
|
|
end findFirst;
|
|
|
|
/* find the greatest common divisor of A and B */
|
|
gcd: procedure(a, b) returns(fixed) recursive;
|
|
declare (a, b) fixed;
|
|
if b = 0 then return(a);
|
|
return(gcd(b, mod(a, b)));
|
|
end gcd;
|
|
|
|
/* calculate S(i) up to MAX */
|
|
declare i fixed;
|
|
S(1) = 1; S(2) = 1;
|
|
do i=2 to MAX/2;
|
|
S(i*2-1) = S(i) + S(i-1);
|
|
S(i*2) = S(i);
|
|
end;
|
|
|
|
/* print first 15 elements */
|
|
put skip list('First 15 elements: ');
|
|
do i=1 to 15;
|
|
put edit(S(i)) (F(2));
|
|
end;
|
|
|
|
/* find first occurrences of 1..10 and 100 */
|
|
do i=1 to 10;
|
|
put skip list('First',i,'at',findFirst(i));
|
|
end;
|
|
put skip list('First ',100,'at',findFirst(100));
|
|
|
|
/* check GCDs of adjacent pairs up to 1000th element */
|
|
do i=2 to 1000;
|
|
if gcd(S(i-1),S(i)) ^= 1 then do;
|
|
put skip list('GCD of adjacent pair not 1 at i=',i);
|
|
stop;
|
|
end;
|
|
end;
|
|
put skip list('All GCDs of adjacent pairs are 1.');
|
|
end sternBrocot;
|