# DigBench / NLA

Programs consisting of nonlinear polynomial invariants.

## A division algorithm (cohendiv)

```int cohendiv(int x, int y){
assert(x>0 && y>0);

int q=0; int r=x;

while (r>=y) {
int a=1; int b=y;
while (r >= 2*b){
//x==q*y+r, b==y*a, , r>=0, r>=2*y*a
a = 2*a; b = 2*b;
}
r=r-b; q=q+a;
}

//x == q*y+r, q == x % y, r == x % y
return q;
}
```

## Another division algorithm (divbin)

```int divbin(int A, int B){
/*
a division algorithm, by Kaldewaij
returns A//B
*/
assert(A > 0 && B > 0);
int q = 0; int r = A; int b = B;
while(1){
if (!(r>=b)) break;
assert(q == 0 && r == A  &&  A >= 0  &&  B > 0);
b=2*b;
}

while(1){
if (!(b!=B)) break;
//(A == q*b + r && r >= 0);
q = 2*q; b = b/2;
if (r >= b) {
q = q + 1; r = r - b;
}
}
//q==A/B, r==A%B
return q;
}
```

## Another division algorithm (hard)

```int hard(int A, int B) {
//hardware integer division program, by Manna
//returns q==A//B

assert(A >= 0);
assert(B >= 1);

int r,d,p,q;

r=A;
d=B;
p=1;
q=0;

while(1){
if (!(r >= d)) break;
// A >= 0 && B > 0 && q == 0 && r == A && d == B*p
d = 2 * d;
p  = 2 * p;
}

while(1){
if (!(p!=1)) break;
// A == q*B+r
// d == B*p
d=d/2;
p=p/2;

if(r>=d){
r=r-d;
q=q+p;
}
}

// r == A % B
// q == A / B
return q;
}
```

## Another division algorithm (mannadiv)

```int mannadiv (int x1, int x2){
//Division algorithm from
//"Z. Manna, Mathematical Theory of Computation, McGraw-Hill, 1974"
//return x1 // x2

int y1,y2,y3;
y1 = 0; y2 = 0; y3 = x1;
while(1){
if(!(y3 != 0)) break;
// y1* x2 + y2 + y3 == x1
if (y2 + 1 == x2) {
y1 = y1 + 1;
y2 = 0;
y3 = y3 - 1;
}
else {
y2 = y2 + 1;
y3 = y3 - 1;
}
}

// y1 == x1 / x2
return y1;
}
```

## Square root of natural number (sqrt)

```int sqrt1(int n){
/* computing the floor of the square root of a natural number */
assert(n >= 0);

int a,s,t;
a=0;
s=1;
t=1;

while(1){

if(!(s <= n)) break;

assert(a*a <= n);
assert(t == 2*a + 1);
assert(s == (a + 1)*(a + 1));

a=a+1;
t=t+2;
s=s+t;
}

// assert(a==(int)floor(sqrt(n)));
return a;
}
```

## Another square root algorithm (dijkstra)

```int dijkstra(int n){
/* program computing the floor of the square root, by Dijkstra */
assert(n>=0);

int p,q,r,h;

p=0;
q=1;
r=n;
h= 0;

while(1){
if (!(q<=n)) break;

q=4*q;
assert(n >= 0 && p == 0 && r==n );
}
//q == 4^n

while(1){

if (!(q!=1)) break;

assert(r < 2*p + q);
assert(p*p + r*q == n*q);

//h^2*p - 4*h*n*q + 4*h*q*r + 4*n*p*q - p*q^2 - 4*p*q*r == 0
//h^2*n - h^2*r - 4*h*n*p + 4*h*p*r + 4*n^2*q - n*q^2 - 8*n*q*r + q^2*r + 4*q*r^2 == 0
//h^3 - 12*h*n*q - h*q^2 + 12*h*q*r + 16*n*p*q - 4*p*q^2 - 16*p*q*r == 0

q=q/4;
h=p+q;
p=p/2;
if (r>=h){
p=p+q;
r=r-h;
}
}

assert(p == (int)floor(sqrt(n)));
return p;
}
```

## Another square root algorithm (freire1)

```int freire1 (float a){
//algorithm for finding the closest integer to square root,
//from  www.pedrofreire.com/crea2_en.htm?

float x;
int r;

x=a/2.0;
r=0;

while(1){
if (!(x>r)) break;
x=x-r;
r=r+1;
}

assert(r==(int)round(sqrt(a)));

return r;
}
```

## Cubic root (freire2)

```int freire2(float a){
//algorithm for finding the closest integer to the cubic root,
//from www.pedrofreire.com/crea2_en.htm? */
//a is a float

float x,s;
int r;

x=a;
r=1;
s=3.25;

while(1){

if (!(x-s > 0.0)) break;

//(4*r*r*r - 6*r*r + 3*r) + (4*x - 4*a) == 1
//4*s - 12*r*r == 1

x = x - s;
s = s + 6 * r + 3;
r = r + 1;
}

//r == (int)round(pow(a,(1.0/3.0)))

return r;
}
```

## Consecutive cube (cohencu)

```int cohencu(int a){
/* printing consecutive cube, by Cohen */

int n,x,y,z;
n=0; x=0; y=1; z=6;

while(1){

if (!(n<=a)) break;
//z == 6*n + 6);
//y == 3*n*n + 3*n + 1);
//x == n*n*n);

n=n+1;
x=x+y;
y=y+z;
z=z+6;
}
return x;
}
```

## Extended GCD (euclidex1)

```int euclidex1 (int x, int y){

int a,b,p,q,r,s;

a=x;
b=y;
p=1;
q=0;
r=0;
s=1;

while(1){
if (!(a!=b)) break;
//1 == p*s - r*q, a == y*r + x*p, b == x*q + y*s
if (a>b){
a = a-b;
p = p-q;
r = r-s;
}
else {
b = b-a;
q = q-p;
s = s-r;}
}

return a;
}
```

## Another extended GCD (euclidex2)

```int euclidex2 (int x, int y){
/* extended Euclid's algorithm */

int a,b,p,q,r,s;

a=x; b=y; p=1; q=0; r=0; s=1;

while(1){

if (!(b!=0)) break;

int c=a,k=0;

while(1){

if (!(c>=b)) break;

//a == k*b+c);
//a == y*r+x*p);
//b == x*q+y*s);

c=c-b;
k=k+1;
}

a=b;
b=c;

int temp;
temp=p;
p=q;
q=temp-q*k;
temp=r;
r=s;
s=temp-s*k;
}

//mygcd(x,y)==y*r+x*p );
return a;
}
```

## Another Extended GCD (xgcd) UNTESTED !

```/* def xgcd(a,b): */
/* if b == 0: */
/* return [1,0,a] */
/* else: */
/* x,y,d = xgcd(b, a%b) */
/* return [y, x - (a//b)*y, d] # Note that a//b is floor(a/b) in Python. */
```

/* Todo: lots of good euclid like algorithms here

## Another Extended GCD (euclidex3)

```int euclidex3(int x, int y){
/* extended Euclid's algorithm */

int a,b,p,q,r,s;

a=x; b=y;  p=1;  q=0;  r=0;   s=1;

//a==y*r+x*p
//b==x*q+y*s

while(1){

if(!(b!=0)) break;

int c=a,k=0;

while(1){
if (!(c>=b)) break;
int d=1,v=b;

while(1){

if (!(c>= 2*v)) break;
//a == y*r+x*p, b == x*q+y*s, a == k*b+c, v == b*d

d = 2*d;
v = 2*v;

}
c=c-v;
k=k+d;

//a==y*r+x*p, b==x*q+y*s, a==k*b+c

}

a=b;
b=c;
int temp;
temp=p;
p=q;
q=temp-q*k;
temp=r;
r=s;
s=temp-s*k;
}

//mygcd(x,y)==y*r+x*p );

return a;
}
```

## GCD and LCM (lcm1)

```int lcm1(int a, int b){
/* algorithm for computing simultaneously the GCD and the LCM,
by Sankaranarayanan */

int x,y,u,v;

x=a;
y=b;
u=b;
v=0;

while(1){
if (!(x!=y)) break;

//x*u + y*v == a*b

while(1){
if (!(x>y)) break;

x=x-y;
v=v+u;
}

while(1){
if (!(x<y)) break;

y=y-x;
u=u+v;
}
}
//x==gcd(a,b)

int r = u+v;
return r; //lcm
}
```

## Another GCD and LCM (lcm2)

```int lcm2 (int a, int b){
/* algorithm for computing simultaneously the GCD and the LCM, by Dijkstra */

int x,y,u,v;

x=a;
y=b;
u=b;
v=a;

while(1){
if (!(x!=y)) break;

//x*u + y*v == 2*a*b
if (x>y){
x=x-y;
v=v+u;
}
else {
y=y-x;
u=u+v;
}

}

//x==gcd(a,b)
int r = (u+v)/2;
return r; //lcm
}
```

## Product of two natural numbers (prodbin)

```int prodbin (int a, int b){
/* algorithm for computing the product of two natural numbers */

assert(b>=0);

int x,y,z;

x = a;
y = b;
z = 0;

while(1){
if (!(y!=0)) break;
// z+x*y==a*b

if (myodd(y)){
z = z + x;
y = y - 1;
}
x = 2 * x;
y = y / 2;

}

// z == a*b
return z;
}

## Another product of two natural numbers (prod4br)

int prod4br(int x, int y){
/* algorithm for computing the product of two natural numbers */

int a,b,p,q;

a = x;
b = y;
p = 1;
q = 0;

while(1){
if (!(a!=0 && b!=0)) break;

//q+a*b*p==x*y)

if (a % 2 ==0 && b % 2 ==0 ){
a =a/2;
b = b/2;
p = 4*p;
}
else if (a % 2 ==1 && b % 2 ==0 ){
a =a-1;
q = q+b*p;
}
else if (a % 2 ==0 && b % 2 ==1 ){
b =b-1;
q = q+a*p;
}
else {
a =a-1;
b =b-1;
q = q+(a+b+1) * p;
}

}

//q == x*y
return q;

}
```

## Divisor for factorization (fermat)

```int fermat1(int A, int R){
// program computing a divisor for factorisation, by Knuth 4.5.4 Alg C ?

// TODO: check for infinite loop

//A >= 1
//(R-1)*(R-1) < A
//A <= R*R
//A%2 ==1

int u,v,r;
u=2*R+1;
v=1;
r=R*R-A;

// 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1
while(1){
if (!(r!=0)) break;

//4*(A+r) == u*u-v*v-2*u+2*v

while(1){
if (!(r>0)) break;
r=r-v;
v=v+2;
}

while(1){
if (!(r<0)) break;
r=r+u;
u=u+2;
}

}

//u != v
int o = (u-v)/2;
return o;
}
```

## Another divisor for factorization (fermat2)

```int fermat2(int A, int R){
/* program computing a divisor for factorisation, by Bressoud */

//A >= 1
//(R-1)*(R-1) < A
//A <= R*R
//A%2 ==1

int u,v,r;

u=2*R+1;
v=1;
r=R*R-A;

//// 4*(A+r)==u*u-v*v-2*u+2*v && v%2==1 && u%2==1 && A>=1
while(1){
if (!(r!=0)) break;

//4*(A + r) == u*u - v*v - 2*u + 2*v

if (r>0){
r=r-v;
v=v+2;
}
else{
r=r+u;
u=u+2;
}
}

//u!=v
int o = (u-v)/2;
return o;
}
```

## Another divisor for factorization (knuth)

```int knuth(int n, int a){
//algorithm searching for a divisor for factorization, by Knuth

int r,k,q,d,s,t;

d=a;
r= n % d;
t = 0;

if (d <= 2){
printf("#E: d - 2 <= 0.  Use a diff val for d\n");
return 0;
}
k=n % (d-2);
q=4*(n/(d-2) - n/d);
s=sqrt(n);

while(1){
//assert(d*d*q - 2*q*d - 4*r*d + 4*k*d  + 8*r == 8*n);
//assert(k*t == t*t);
//assert(d*d*q - 2*d*q - 4*d*r + 4*d*t + 4*a*k - 4*a*t - 8*n + 8*r == 0);
//assert(d*k - d*t - a*k + a*t == 0);
if (!((s>=d)&&(r!=0))) break;

if (2*r-k+q<0){
t=r;
r=2*r-k+q+d+2;
k=t;
q=q+4;
d=d+2;
}
else if ((2*r-k+q>=0)&&(2*r-k+q<d+2)){
t=r;
r=2*r-k+q;
k=t;
d=d+2;
}
else if ((2*r-k+q>=0)&&(2*r-k+q>=d+2)&&(2*r-k+q<2*d+4)){
t=r;
r=2*r-k+q-d-2;
k=t;
q=q-4;
d=d+2;
}
else {/* ((2*r-k+q>=0)&&(2*r-k+q>=2*d+4)) */
t=r;
r=2*r-k+q-2*d-4;
k=t;
q=q-8;
d=d+2;
}

}
return d;
}
```
```void petter(k){
//geo
int x = 0;
int y = 0;

while(1){
if (!(y != k)) break;

x = x + pow(y,5);
y = y + 1;
}
}
```
```/*from Petter"s master thesis
http://www2.cs.tum.edu/~petter/da/da.pdf

IMPORTANT: don't try with very large power (k)
which makes the number *very* large and becomes problematic for computers.
e.g. 15*8 = -1732076671 on my Macbook Pro
This is a _computer_dependent_problem, not mathematics problem
*/

int geo_pw(const int a, const int z, const int k){
//computes (z-1) * sum(z^i)[k=0..k-1]

int i = 0; int s = 0;
for(; i<k ; ++i) s = s + a*pow(z,i);
return s;
}
```

## Geometric Series (geo1)

```int geo1(const int z, const int k){
/* computes x=(z-1)* sum(z^k)[k=0..k-1] , y = z^k
returns 1+x-y == 0

The loop is same as geo2 but the result is different.
*/

assert(k>0);

int x = 1; int y = z; int c = 1;

while(1){
if (!(c < k)) break;

c = c + 1;
x = x*z + 1;
y = y*z;

}

x = x *(z - 1);

//y == pow(z,k)
//geo_pw(1,z,k)*(z-1) == x
//1+x-y == 0 //documented inv

return x;
}
```

## Another geometric series (geo2)

```int geo2(const int z, const int k){
//computes x = sum(z^k)[k=0..k-1], y = z^(k-1)
assert (k>0);

int x = 1; int y = 1; int c = 1;

while(1){
if (!(c < k)) break;

c = c + 1;
x = x*z + 1;
y = y*z;
}

//y == pow(z,(k-1))
//geo_pw(1,z,k) == x);

return x;
}
```

## Another geometric series (geo3)

```int geo3(const int z, const int a, const int k){
//computes x = sum(a*z^k)[k=0..k-1], y = z^(k-1)
assert (k>0);

int x = a; int y = 1;  int c = 1;
while(1){
if (!(c < k)) break;

c = c + 1;
x = x*z + a;
y = y*z;

}

//y == pow(z,(k-1))
//geo_pw(a,z,k) == x

return x;
}
```
```/*
# // yielding x=y
# //todo: too easy, not worth doing
# def ps1 ( String []
# y = 0
# x = 0
# while( args !=null){
# y=y +1
# x=x +1
*/

int ps_pw(const int k, const int z){
// computes sum(i^z)[k=1..k]
int i = 1; int s = 0;
for(; i <= k; ++i) s = s + pow(i,z);
return s;
}
```

## Power Sum (ps1)

```int ps1(const int k){
int y = 0; int x = 0;int c = 0;
while(1){
if (!(c < k)) break;

c = c + 1;
y = y + 1;
x = x + 1;
}

//y == 1*k
//x == ps_pw(k,0)
return x;
}
```

## Power Sum (ps2)

```int ps2(int k){

int y = 0; int x = 0; int c = 0;

while(1){
//assert((y*y) - 2*x + y == 0);
if (!(c < k)) break;

c = c + 1;
y=y +1;
x=x+y;

}

//y == 1*k
//x == ps_pw(k,1)
return x;
}
```

## Power Sum (ps3)

```int ps3(int k){
int y = 0;
int x = 0;
int c = 0;

while(1){
if (!(c < k)) break;

c = c + 1;
y=y +1;
x=y*y+x;

}

//y == 1*k
//x == ps_pw(k,2)
return x;
}
```

## Power Sum (ps4)

```int ps4 (int k){
int y = 0;
int x = 0;
int c = 0;

while(1){
//assert(4*x-(y*y*y*y)-2*(y*y*y)-(y*y) == 0);

if (!(c < k)) break;

c = c +1 ;
y=y +1;
x=y*y*y+x;

}

//y == 1*k
//x == ps_pw(k,3)
return x;
}
```

## Power Sum (ps5)

```int ps5 (int k){
int y = 0;
int x = 0;
int c = 0;

while(1){
if (!(c < k)) break;

c = c +1 ;
y=y +1;
x=y*y*y*y+x;
}

//y == 1*k
//x == ps_pw(k,4)
return x;
}
```

## Power Sum (ps6)

```int ps6 (int k){
int y = 0;
int x = 0;
int c = 0;

while(1){
if (!(c < k)) break;

c = c +1 ;
y=y +1;
x=y*y*y*y*y+x;
}

//y == 1*k
//x == ps_pw(k,5)
return x;
}
```

Updated