ちょっとだけ補足します。この三つのプログラムいずれにも、
>> p[0]= - s1*a -2*s2*b + 2*c;
>> p[1]= -2*s1*a - s2*b + 2*c;
>> p[2]= -2*s1*a -2*s2*b + 3*c;
のようなものが登場します。初見だと、何のことだか全く判らないかもしれないので、説明を加えます。

原理は、単純です。次の問題を考えてみてください。
「直角三角形があります。斜辺とある一辺の長さの差は１、斜辺と残りの一辺の差は２です。三辺は何か？」
答えはもちろん　3,4,5　です。では、次の問題は？
「直角三角形があり、斜辺とある一辺の長さの差はｘ、斜辺と残りの一辺の差はｙ。直角三角形の三辺を、ｘとｙで表せ。」
あるいは、「a^2+b^2=c^2、x=c-a、y=c-b　の関係があるとき、a,b,cをx,yで表せ」

この問題を解くことにより　(y±√(2xy))^2 + (x±√(2xy))^2 　=　(x+y±√(2xy))^2　；　(複合同順)
という恒等式を見いだすことができ、2xy=(a+b-c)^2に注意して、書き換えれば、
(c-b±(a+b-c))^2 + (c-a±(a+b-c))^2　=(2c-a-b±(a+b-c))^2　が得られます。
プラス側からは、a^2+b^2=c^2という、前提にしていたつまらない式が現れますが、マイナス側からは、
(-a-2b+2c)^2 + (-2a-b+2c)^2 = (-2a-2b+3c)^2　・・・　（★）
を得ます。これがプログラムに登場する式の骨格です。a^2+b^2=c^2　が成立するときに（★）が成立するなら、
aを-aに置き換えた、( a-2b+2c)^2 + ( 2a-b+2c)^2 = ( 2a-2b+3c)^2 も成立し、．．．というのが、
ある原始ピタゴラス数から、別の原始ピタゴラス数を生成するアルゴリズムとなっています。