Rational arithmetic algorithms for Quote notation

[luser droog]

As described in my new readme,
https://github.com/luser-dr00g/inca/tree/master/olmec
I'm trying to build an interpreter/calculator based on
Quote Notation https://en.wikipedia.org/wiki/Quote_notation .

Quote notation represents positive integers as a sequence of
digits, ideally in a prime base, with an additional zero figure
marked as repeating infinitely to the left.

Negative integers are represented in a radix-complement form
with a repeating figure of base-1. For decimal base, this
is "10's complement" ie. 9's comp +1. For binary base, this
is the regular 2's complement, with an infinitely repeating
1 bit.

The representation adapts to more general rational numbers,
but for simplicity I'm going to just build it up for integers
first.

I intend to represent this in C with a struct.

struct q {
int n,m;
int d[]; // d[0]..d[n-1] d[n]..d[n+m]
};

So, before basic functions like addition and subtraction
can be defined (easily), we need some helper functions for
making two number forms *conformable* for the algorithm.

Addition and subtraction proceed columnwise digit-by-digit
(see examples in the Wikipedia article), so conformable
means x->n == y->n && x->m == y->m. Otherwise, the smaller
one needs to be *unrolled* from the repeating m portion.

So a first draft for these:

void match_n(q**x,q**y){
if ((*x)->n == (*y)->n)
return;
if ((*x)->n < (*y)->n)
roll_out(x, (*y)->n);
else
roll_out(y, (*x)->n);
}

void roll_out(q**x,int n){
q*z=malloc(sizeof*z+(n+(*x)->m)*sizeof*z->d);
z->n=n;
z->m=(*x)->m;
memcpy(z->d,(*x)->d,((*x)->n+(*x)->m)*sizeof*z->d);
int i;
int *d,*s;
for (i=0,s=&z->d[(*x)->n],d=&z->d[(*x)->n+(*x)->m];
i < n-(*x)->n;
i++)
*d++=*s++;
*x=z;
}

void match_m(q**x,q**y){
if ((*x)->m == (*y)->m)
return;
if ((*x)->m < (*y)->m)
reshape_m(x,(*y)->m);
else
reshape_m(y,(*x)->m);
}

void reshape_m(q**x,int m){
q*z=malloc(sizeof*z+((*x)->n+m)*sizeof*z->d);
z->n=(*x)->n;
z->m=m;
memcpy(z->d,(*x)->d,((*x)->n+(*x)->m)*sizeof*z->d);
int i;
int *d,*s;
for (i=0,s=&z->d[(*x)->n],d=&z->d[(*x)->n+(*x)->m];
i < m-(*x)->m;
i++)
*d++=*s++;
*x=z;
}

Then for just integers, still,

q*qadd(q*x,q*y){
int carry,j;
match_n(&x,&y);
match_m(&x,&y)
for (carry=0,j=0; j < z->m+z->n; j++){
int r= x->d[j] + y->d[j] + carry;
z->d[j] = r % BASE;
carry = r / BASE;
}
return z;
}

Then we need a helper function to normalize a result
that's excessively "rolled-out". Still working on that one.
Comments appreciated.

[Jeff Higgins]

On 08/05/2015 07:57 PM, luser droog wrote:
> As described in my new readme,
> https://github.com/luser-dr00g/inca/tree/master/olmec

Exact arithmetic on the Stern–Brocot tree
<http://www.sciencedirect.com/science/article/pii/S1570866706000311>

Fractions and p-adic numbers
<http://cosmolearning.org/video-lectures/fractions-p-adic-numbers/>

[luser droog]

On Thursday, August 6, 2015 at 1:54:51 PM UTC-5, Jeff Higgins wrote:
> On 08/05/2015 07:57 PM, luser droog wrote:
> > As described in my new readme,
> > https://github.com/luser-dr00g/inca/tree/master/olmec
>

> Exact arithmetic on the Stern-Brocot tree

> <http://www.sciencedirect.com/science/article/pii/S1570866706000311>
>
> Fractions and p-adic numbers
> <http://cosmolearning.org/video-lectures/fractions-p-adic-numbers/>

Thanks for those. The video really helps.

I've ordered a pile of books on the subject so I'm kinda stalled
until they arrive. Here are two drafts that may not work,
and are definitely incomplete. But it's something to look at.

$ cat q.c
#include <math.h>
#include <stdlib.h>
#include <string.h>

typedef struct field {
int (*add)();
int (*mul)();
int (*le)();
} *field;

typedef struct B {
field f;
int b;
} B;

typedef struct I {
field f;
int b,d[1];
} I;

typedef struct Qpq {
field f;
I p,q;
} Qpq;

typedef struct Qnmd {
field f;
int n,m,b,d[1];
} Qnmd;

typedef struct Qenmd {
field f;
int e,n,m,b,d[1];
} Qenmd;

Qnmd *qneg(Qnmd *q){
Qnmd *z;
z=malloc(sizeof*z + (q->n+q->m)*sizeof*z->d);
*z=*q;

for (int j=0;j<=q->n;j++){
z->d[j]= q->b-1 - q->d[j];
}
return z;
}

Qnmd *Qnmd_(int i){
Qnmd *t=&(Qnmd){.b=10};
if (i < 0) return qneg(Qnmd_(-i));
t->m=0;
t->n=ceil(log(i==0?1:i)/log(t->b));
if (t->n==0) t->n=1;

Qnmd *z;
z=malloc(sizeof*z + t->n*sizeof*z->d);
*z=*t;

for (int j=0; j<=z->n; j++){
z->d[j]= i % z->b;
i /= z->b;
}
return z;
}

#include <stdio.h>

void prnq(Qnmd *q){
printf("(%d, %d, (", q->n,q->m);
for (int j=0; j< q->n+q->m; j++)
printf("%s%d", j?", ":"", q->d[j]);
printf(")\n");
}

int main() {
for (int j=0; j<10; j++){
Qnmd *q=Qnmd_(j);
prnq(q);
}

for (int j=0; j<10; j++){
Qnmd *q=Qnmd_(-(j+1));
prnq(q);
}
}

$ cat qn.c
enum { base = 10 };

/*
d[n+m-1] ... d[n+1] d[n] ' d[n-1] ... d[2] d[1] d[0]

Sum(i=0..n d[i]*b^i) + Sum(i=n+1..n+m d[i]*b^i/(b^m-1))
*/
typedef struct q {
int b,n,m,d[];
} *q;

q qneg(q x){
return qsub(qlong(0),x);
}

q qlong(long i){
if (i<0) return qneg(qlong(-i));
int m=0;
int n=ceiling(log(i)/log(base));
q z=malloc(sizeof*z + (n+m)*sizeof*z->d);
z->b=base; z->n=n; z->m=m;
for(int j=0; j<n; j++){
z->d[j] = i % base;
i /= base;
}
return z;
}

int qsign(q x){
switch(((x->n==0))<<1+(x->m==0)){
case (1<<1)+1: return 0;
case 1<<1: return -1;
case 1: return 1;
case 0: return x->d[x->n-1] - x->d[x->n+x->m-1];
}
}

q qadd(q x,q y){
int m=lcm(x->m,y->m);
int n=max(x->n,y->n)+m+2;
int t,c=0,tt=0;
q z=malloc(sizeof*z + (n+m)*sizeof*z->d);
z->b=base; z->n=n; z->m=m;
for(int k=0; k<n+m; k++){
t = x->d[k] + y->d[k] + c;
z->d[k] = t % base;
c = t / base;
tt=t;
}
}

q qsub(q x,q y){
int m=lcm(x->m,y->m);
int n=max(x->n,y->n)+m+2;
q z=malloc(sizeof*z + (n+m)*sizeof*z->d);
z->b=base; z->n=n; z->m=m;
}

q qmul(q x,q y){
int m=lcm(x->m,y->m);
int n=x->n+y->n+m+1;
q z=malloc(sizeof*z + (n+m)*sizeof*z->d);
z->b=base; z->n=n; z->m=m;
}

q qdiv(q x,q y){
int m=lcm(x->m,euler_phi(y));
int n= y->n<=y->m ? x->n + y->n + m + 1 :
y->m<y->n && y->n<=x->n ? x->n - y->n + m + 1:
m+1;
q z=malloc(sizeof*z + (n+m)*sizeof*z->d);
z->b=base; z->n=n; z->m=m;
}


fwiw

[Bill Davy]

I should get a life really, but ...

q qlong(long i){
if (i<0) return qneg(qlong(-i));
int m=0;
int n=ceiling(log(i)/log(base));

Seems likely to have probelms for MIN_INT (for example, 0x8000 which
is -32768, and -0x8000 is -32768 as well) and for zero (log(0) feels large).

Assuming 2's complement of course.

And if it fails for two simple edge cases, then I would suspect there is a
lot more to find.


"luser droog" <luser...@gmail.com> wrote in message
news:103d1a31-db8f-46f2...@googlegroups.com...