@s mod and \let\Xmod=\bmod % this is CWEB magic for using "mod" instead of "%" @*Intro. This program computes a function suggested by Roland Bacher's article in {\sl American Mathematical Monthly\/ \bf130} (2023), 824--836, namely the number of solutions to $n=ab+cd$ where $\max(a,b)<\min(c,d)$. Bacher proved the astonishing fact that there are exactly $(p+1)/2$ nonnegative solutions when $n=p$ is an odd prime. I'm curious about the number of solutions for other values of~$n$. The goal of this program is to compute a ``b-file'' for this sequence, whieh has been assigned {\sl OEIS\/} number A368207. @d mod % /* used for percent signs denoting remainder in \CEE/ */ @d maxn 10001 @c #include #include int th[maxn]; int xs[maxn]; int sum; @; int main(void) { register int d,dd,i,j,k,kk,m,n,nn; printf("# computed by Don Knuth's CWEB program on 18 December 2023\n"); @; for (n=1;n; printf("\n"); /* a b-file ends with a blank line */ } @ The function |th[n]| is the largest divisor $d$ of $n$ whose square is less than or equal to~$n$. (Thus $n=dd'$ where $d\le d'$. It's sequence A0033676 in the {\sl OEIS}.) The function |xs[n]| is twice the sum of $4d-2$ over all such divisors $d$, except that we use only $2d-1$ when $d^2=n$. In the following loop, $m^2\le n\le|nn|=(m+1)^2-1$. @= for (m=n=1,nn=3;n= { sum=0; for (k=1,kk=n-1;ki;j--) if (!(kk mod j)) solution(k/i,i,j,kk/j); } } } printf("%d %d\n", n,sum+xs[n]); } @ @= void solution(int a,int b,int c,int d) { if (a==b) { if (c==d) sum++; else sum+=2; } else if (c==d) sum+=2; else sum+=4; } @*Index.