pseudobase-prime and the squid operator, revisited

[os]

required reading
gnu-bc implementation for the squid() operator follows:



/* example usage: */
/* j=num2p(200,f[]) --> 3 0 2 in f[] */
/* n=p2num(f[],j) returns 200 */

define num2p(n, *f[]) {
auto i, j, l, s[], z;
z = scale;
scale = 0;
l = n;
for(i = 2; i <= l; ++i) {
s[i - 2] = 1;
}
for(i = 2; i <= l; ++i) {
for(j = 2; i*j <= l; ++j) {
s[i*j - 2] = 0;
}
}
for(i = 2; i <= l; ++i) {
print i, ": ", s[i - 2], "\n";
}
j = 0;
for(i = 2; i <= l; ++i) {
if(s[i - 2]) {
f[j] = 0;
print n, " % ", i, " == ", n % i;
while(n % i == 0) {
++f[j];
n = n / i;
print "for ", i, " (p#", j, "): n -> ", n, "; f -> ", f[j], "\n";
}
++j;
}
}
for(i = 0; i < j; ++i) {
print f[i], " "
}
print "\n"
scale = z;
return j;
}

define p2num(p[], r) {
auto i, j, l, s[], z, o, k;
z = scale;
scale = 0;
l = r;
k = 0;
while(k < r) {
l = l * 2;
for(i = 2; i <= l; ++i) {
s[i - 2] = 1;
}
for(i = 2; i <= l; ++i) {
for(j = 2; i*j <= l; ++j) {
s[i*j - 2] = 0;
}
}
k = 0;
for(i = 2; i <= l; ++i) {
if(s[i - 2]) {
++k;
}
}
}
print "l=", l, " k=", k, "\n";
for(i = 2; i <= l; ++i) {
print i, ": ", s[i - 2], "\n";
}
o = 1;
j = 0;
for(i = 2; i <= l; ++i) {
if(s[i - 2]) {
if(p[j]) {
o = o * i^p[j];
print "after ", i, " (p#", j, "): o -> ", o, " from p:", p[j], "\n";
}
++j;
}
}
scale = z;
return o;
}

define squid(a,b) {
auto m, n, p[], q[], r[], j, k, l;
j = num2p(a,p[]);
k = num2p(b,q[]);
l = j + k;
print "j:", j, " k:", k, " -> l:", l, "\n";
for(m = 0; m < l; ++m) {
r[m] = 0;
}
for(m = 0; m < j; ++m) {
for(n = 0; n < k; ++n) {
r[m + n] = r[m + n] + p[m]*q[n];
print "@", m, ",", n, " r->", r[m+n], " from ", p[m], ",", q[n], "\n";
}
}
return p2num(r[],l);
}

Yes, it really does construct a sieve of eratosthenes each time, and you do have to manually have to track the upper bound of how many primes to use in the table (returned from num2p, passed into p2num). It's just a quick hack, but I think it's still faster than the shell script version (see "required reading" link).