Kawazoe(@riverplus)さんのCodeIQの問題です。
問題はこちらの御世話になりました。解答例もあります。
->https://blog.goo.ne.jp/r-de-r/e/45659ad975496af9b1b7de04a8d2305d
->https://blog.goo.ne.jp/r-de-r/e/44e77fc0784eef7457ffdd88c8771609
こちらにも情報があります。
->https://posfie.com/@riverplus/p/SeiCU5C

求める外接三角形の辺の長さは原始ピタゴラス数の整数倍です。
原始ピタゴラス数は3つの数が互いに素（ある2数が互いに素)なピタゴラス数です。
さらに整数辺の直角三角形の公式(a,b,c)=(I^2-J^2,2*I*J,I^2+J^2)(aが奇数)と、
三辺が(a1,b1,c1)の外接直角三角形の斜辺が、内接円の半径をrとしますと、
c1=a1-r+b1-r=a1+b1-c　(2*r=c)をとなるのを利用します。
この時I,Jは互いに素かつI>Jとなります。

Prologです。
c^2=I^2+J^2<=2*L^2よりI<=√(L/√(2))
先ずloop1/6で1からNまで互いに素でない2数、a,b>Lとなる数、偶数のaを除き、
loop2/6で原始ピタゴラス数をk倍した数が題意を満たすかどうか調べます。
すなわち外接三角形の面積の倍 (a1+b1+(a1+b1-2*r))*r=a1*b1 を変形して
a1=r*(a+b+r)/(2*b),b1=r*(a+b+r)/(2*a)が整数なら返り値を1増やします。

%SWI-Prolog ver 7.42,7.64,9.25
%:-initialization(start).       %Online Prolog Compiler,ideone.    
%start.

loop2(0,_,_,_,R,R).
loop2(K,X,Y,Z,R1,R):-
   X1 is X*K,Y1 is Y*K,Z1 is Z*K,V is (X1 + Y1 + Z1) * Z1,
   ((mod(V,X1*2)=:=0,mod(V,Y1*2)=:=0)->R2 is R1+1;R2 = R1),
   K1 is K-1,loop2(K1,X,Y,Z,R2,R).
   
loop1(N0,N,I,_,R,R):-
   N<I,!.
loop1(N0,N,I,J,R1,R):-
   J =:= I,!,I1 is I+1,loop1(N0,N,I1,1,R1,R).
loop1(N0,N,I,J,R1,R):-
   gcd(I,J)=\=1,!,J1 is J+1,loop1(N0,N,I,J1,R1,R).
loop1(N0,N,I,J,R1,R):-
   X0 is I^2-J^2,J1 is J+1,((0=:=mod(X0,2);X0>N0)->loop1(N0,N,I,J1,R1,R);
  (Y0 is 2*I*J,(Y0>N0->(I1 is I+1,loop1(N0,N,I1,1,R1,R));
  (Z0 is I^2+J^2,K is min(div(N0,X0),div(N0,Y0)),loop2(K,X0,Y0,Z0,R1,R2),
   loop1(N0,N,I,J1,R2,R))))).

solve(N0,R):-N is floor(sqrt(N0/sqrt(2))),loop1(N0,N,2,1,0,R).

start:-
    f(N0,A),(solve(N0,R)->(R==A->Str=" ok  ";Str=" no  ");
    write(fail)),write(" "),write(Str),writeln(N0->R),fail.
start.

f(10,1).
f(50,7).
f(100,15).
f(456,77).
f(1000, 173).      %抜けていたので追加（2025.4.5）
f(7654,1405).
f(75319,14133).
f(90346,16972).
f(100000,18797).

以前juliaで解いていたのですが、投稿しようとして見直したときに、
自分でも何をしているのかわかりませんでした。
そこで新たにプログラムを書きましたが、
元のほうがやや速いようですのでもう一度考えてみました。

求める外接三角形の辺の長さは原始ピタゴラス数の整数倍ですが、
原始ピタゴラス数とその奇数倍のピタゴラス数では
a,cが奇数、bが偶数ですので、a*cは奇数,a+b-cは偶数となり
a*c/(a+b-c)が整数となりませんので、題意に会う外接三角形はありません。

a<bの場合を除外するとa<bとなる原始ピタゴラス数が除かれますが、
原始ピタゴラス数の二倍のピタゴラス数（2*a,2*b,2*c）があるとすれば
2*(i^2-j^2)=(h＊i)^2-(2*j/h)^2からhを求めるとh=√2。これはあり得ません。
一方
2*b=2*2*i*j=(i+j)^2-(i-j)^2,2*a=2*(i^2-j^2)=2*(i+j)(i-j),
2*b>2*aなので、ピタゴラス数(2*b,2*a,2*c)が求められます。
ただしI<=√(L/√(2))のときi1=i+jとするとi<=Iなのでi1<=2*Iです。

またa1/a=b1/b=(a1+b1-c)/cからa1=a*c/(a+b-c),b1=b*c/(a+b-c)となります。

よってi<=n,j<iとなるi,jで、
i,jが互いに素でない場合,a>Lの場合,a<bの場合を除外し、
(a,b,c)のk倍の(a1,b1,c1)に関してa1*c1/(a1+b1-c1),b1*c1/(a1+b1-c1)
が整数なら返り値に1を足します。

#julia ver 1.41,1.53,1.10.3

function solve(n0,v)
    n = Int(floor(sqrt(sqrt(2)*n0))) 
    sum = 0
    for i in 2 : n
     for j in 1 : i-1
       if gcd(i,j) !=  1 
          continue　　　
       end
       a = i^2 - j^2　　
       if a>n0 　　　　　
          continue
       end
       b = 2*i*j
       if b>a  
         break
       end
       c = i^2 + j^2
       for k in 1 : div(n0,a)
         a1 = a*k
         b1  =b*k
         c1 = c*k
         am = rem(a1*c1, a1+b1-c1)  
         bm = rem(b1*c1, a1+b1-c1)　
         if am == 0 && bm == 0
           sum+=1
         end
        end
      end
    end
     s = v == sum ? "ok " : "no "
     println(s,n0,"->",sum)
end

solve(10, 1)
solve(50, 7)
solve(100, 15)
solve(456, 77)
solve(1000, 173)
solve(7654, 1405)
solve(75319, 14133)
solve(90346, 16972)
solve(100000, 18797)