Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 92 additions & 0 deletions lagoon/vere/noun/jets/i/lagoon.c
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "noun.h"
#include "softfloat.h"
#include "softblas.h"
#include "jets/i/twoc.h" // shared two's-complement kernels (%int2 array ops)

#include <math.h> // for pow()
#include <stdio.h>
Expand Down Expand Up @@ -2199,6 +2200,67 @@
return u3nc(u3nq(u3nt(M_, P_, u3_nul), u3k(bloq), c3__i754, u3_nul), r_data);
}

/* %int2 element-wise binary op over a ray: native two's-complement per lane
** (lane width = 2^bloq bits, always a machine width 8..128, so all native --
** c3_b/s/w/d/__int128). Mirrors the i754 byte-marshalling (the data atom's
** leading 0x1 sentinel sits just past the last lane and is preserved). add/
** sub/mul are direct typed wraps; div/mod use the shared signed twoc kernels
** and DECLINE (u3_none -> Hoon, which crashes) on a zero divisor. Result is
** the new data atom; the caller re-wraps it with the (unchanged) meta.
*/
typedef enum { _LA_ADD, _LA_SUB, _LA_MUL, _LA_DIV, _LA_REM } _la_iop;

static u3_noun
_la_int2_binop(u3_noun x_data,
u3_noun y_data,
u3_noun shape,
u3_noun bloq,
_la_iop op)
{
c3_d bl = u3x_atom(bloq);
if ( bl < 3 || bl > 7 ) {
return u3_none;
}
c3_d len = _get_length(shape);
c3_d lb = (c3_d)1 << (bl - 3); // bytes per lane
c3_d syz = len * lb;
c3_d w = (c3_d)1 << bl; // lane width in bits

c3_y* xb = (c3_y*)u3a_malloc(syz * sizeof(c3_y));
c3_y* yb = (c3_y*)u3a_malloc((syz + 1) * sizeof(c3_y));
u3r_bytes(0, syz, xb, x_data);
u3r_bytes(0, syz + 1, yb, y_data);

#define _LA_EW(TY, MSK, SB, DIVFN, REMFN) \
{ TY* X = (TY*)xb; TY* Y = (TY*)yb; \
for ( c3_d i = 0; i < len; i++ ) { \
switch ( op ) { \
case _LA_ADD: Y[i] = (TY)(X[i] + Y[i]); break; \
case _LA_SUB: Y[i] = (TY)(X[i] - Y[i]); break; \
case _LA_MUL: Y[i] = (TY)(X[i] * Y[i]); break; \
case _LA_DIV: \
if ( 0 == Y[i] ) { u3a_free(xb); u3a_free(yb); return u3_none; } \
Y[i] = (TY)DIVFN(X[i], Y[i], (MSK), (SB)); break; \
case _LA_REM: \
if ( 0 == Y[i] ) { u3a_free(xb); u3a_free(yb); return u3_none; } \
Y[i] = (TY)REMFN(X[i], Y[i], (MSK), (SB)); break; \
} } }

switch ( bl ) {
case 3: _LA_EW(c3_b, _tn_msk64(w), w - 1, _tn_div64, _tn_rem64); break;
case 4: _LA_EW(c3_s, _tn_msk64(w), w - 1, _tn_div64, _tn_rem64); break;
case 5: _LA_EW(c3_w, _tn_msk64(w), w - 1, _tn_div64, _tn_rem64); break;
case 6: _LA_EW(c3_d, _tn_msk64(w), w - 1, _tn_div64, _tn_rem64); break;
case 7: _LA_EW(_twoc_u128, _tn_msk128(w), w - 1, _tn_div128, _tn_rem128); break;
}
#undef _LA_EW

u3_noun r_data = u3i_bytes((syz + 1) * sizeof(c3_y), yb);
u3a_free(xb);
u3a_free(yb);
return r_data;
}

u3_noun
u3wi_la_add(u3_noun cor)
{
Expand Down Expand Up @@ -2239,6 +2301,12 @@
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);

case c3__int2: {
u3_noun r_data = _la_int2_binop(x_data, y_data, x_shape, x_bloq, _LA_ADD);
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);
}

default:
return u3_none;
}
Expand Down Expand Up @@ -2286,6 +2354,12 @@
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);

case c3__int2: {
u3_noun r_data = _la_int2_binop(x_data, y_data, x_shape, x_bloq, _LA_SUB);
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);
}

default:
return u3_none;
}
Expand Down Expand Up @@ -2333,6 +2407,12 @@
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);

case c3__int2: {
u3_noun r_data = _la_int2_binop(x_data, y_data, x_shape, x_bloq, _LA_MUL);
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);
}

default:
return u3_none;
}
Expand Down Expand Up @@ -2380,6 +2460,12 @@
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);

case c3__int2: {
u3_noun r_data = _la_int2_binop(x_data, y_data, x_shape, x_bloq, _LA_DIV);
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);
}

default:
return u3_none;
}
Expand Down Expand Up @@ -2427,6 +2513,12 @@
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);

case c3__int2: {
u3_noun r_data = _la_int2_binop(x_data, y_data, x_shape, x_bloq, _LA_REM);
if (r_data == u3_none) { return u3_none; }
return u3nc(u3nq(u3k(x_shape), u3k(x_bloq), u3k(x_kind), u3k(x_tail)), r_data);
}

default:
return u3_none;
}
Expand Down
22 changes: 14 additions & 8 deletions libmath/desk/lib/twoc.hoon
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
~% %non ..part ~ :: jet registration; nest non in hex (cf /lib/unum, /lib/math)
|%
:: Two's-complement signed integers, MODULAR (operations wrap mod 2^width
:: rather than crashing on overflow -- a deliberate divergence from a
Expand All @@ -11,6 +12,7 @@
:: width (bex bloq), so the logic is defined once in +twid.
::
++ twid
~/ %twid
|_ wid=@
:: +len: bit width. +len-mod: the modulus 2^width.
++ len wid
Expand Down Expand Up @@ -42,14 +44,14 @@
:: Arithmetic. Because the low `len` bits of a sum/product are independent
:: of sign, the same bit-level add and multiply serve signed and unsigned.
::
++ add |=([a=@ b=@] ^-(@ (mod (^add a b) len-mod)))
++ add ~/ %add |=([a=@ b=@] ^-(@ (mod (^add a b) len-mod)))
:: +neg: two's-complement negation, ~a + 1 (flip the len bits, add one).
:: neg(min) = min (wraps).
++ neg |=(a=@ ^-(@ (mod +((not 0 len (mod a len-mod))) len-mod)))
++ sub |=([a=@ b=@] ^-(@ (add a (neg b))))
++ mul |=([a=@ b=@] ^-(@ (mod (^mul (mod a len-mod) (mod b len-mod)) len-mod)))
++ neg ~/ %neg |=(a=@ ^-(@ (mod +((not 0 len (mod a len-mod))) len-mod)))
++ sub ~/ %sub |=([a=@ b=@] ^-(@ (add a (neg b))))
++ mul ~/ %mul |=([a=@ b=@] ^-(@ (mod (^mul (mod a len-mod) (mod b len-mod)) len-mod)))
:: +abs: |a| as a two's-complement value (abs(min) wraps back to min).
++ abs |=(a=@ ^-(@ ?:(=(1 (msb a)) (neg a) (mod a len-mod))))
++ abs ~/ %abs |=(a=@ ^-(@ ?:(=(1 (msb a)) (neg a) (mod a len-mod))))
:: +overflow: would a + b overflow a CHECKED signed add? Retained for
:: callers that want to detect rather than wrap; +add no longer uses it.
++ overflow
Expand All @@ -60,6 +62,7 @@
:: +div: signed division, truncating toward zero (C / Hoon `div` style).
:: Division by zero crashes (as `div` does). div(min, -1) wraps to min.
++ div
~/ %div
|= [a=@ b=@]
^- @
=/ sa (msb a)
Expand All @@ -70,13 +73,15 @@
:: +rem: signed remainder, sign follows the DIVIDEND (C `%` / truncated).
:: a == (add (mul (div a b) b) (rem a b)).
++ rem
~/ %rem
|= [a=@ b=@]
^- @
=/ r (mod (abs a) (abs b))
?: =(0 (msb a)) (mod r len-mod)
(neg r)
:: +pow: a to a NON-NEGATIVE integer power n (raw @), by modular squaring.
++ pow
~/ %pow
|= [a=@ n=@]
^- @
=/ base (mod a len-mod)
Expand All @@ -94,13 +99,14 @@
(dec (bex len))
:: Comparisons, by two's-complement order (negatives below non-negatives).
++ gth
~/ %gth
|= [a=@ b=@]
?: =(1 (mix (msb a) (msb b)))
=(0 (msb a))
(^gth a b)
++ lth |=([a=@ b=@] (gth b a))
++ lte |=([a=@ b=@] !(gth a b))
++ gte |=([a=@ b=@] !(gth b a))
++ lth ~/ %lth |=([a=@ b=@] (gth b a))
++ lte ~/ %lte |=([a=@ b=@] !(gth a b))
++ gte ~/ %gte |=([a=@ b=@] !(gth b a))
--
::
:: +twoc: bloq-keyed facade over +twid (width = 2^bloq). Logic lives in
Expand Down
23 changes: 23 additions & 0 deletions libmath/tools/twoc/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
# `/lib/twoc` jet validation harnesses

Standalone C programs that check the arithmetic in the `/lib/twoc` jet
(`libmath/vere/noun/jets/i/twoc.c`) against an independent two's-complement
reference, BEFORE any ship build. Each copies the jet's arithmetic VERBATIM and
compares to a from-scratch reference built on hardware `__int128` / `mpz`.

- `twoc_native_test.c` — the native paths (`c3_d` for wid<=64, `__int128` for
wid<=128). Widths 8/16/17/32/64/100, all sign edges (min, -1, max, 0), 12
ops. **8280 checks.**
- `twoc_gmp_test.c` — the GNU MP path, cross-checked vs the same hardware
reference at widths 8..100, plus hand-verified 200-bit cases GMP-only.
**5692 checks.**

## Run

```
cc -O2 -Wall -Wno-unused-function twoc_native_test.c -o /tmp/t && /tmp/t
cc -O2 -Wall -Wno-unused-function -I/opt/homebrew/include -L/opt/homebrew/lib \
twoc_gmp_test.c -lgmp -o /tmp/tg && /tmp/tg
```

Both print `checks=<N> fails=0` on success.
106 changes: 106 additions & 0 deletions libmath/tools/twoc/twoc_gmp_test.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
// Validate the GMP (arbitrary-precision) path of twoc.c.
// The _tg_* helpers + op blocks are copied VERBATIM from the jet; we run them
// at widths 8..100 and compare to the independent hardware u128/s128 reference
// (the same one that validated the native path), then hand-check a few >128
// cases that no hardware type can hold.
#include <stdio.h>
#include <gmp.h>

typedef unsigned long long c3_d;
typedef unsigned char c3_t;
typedef unsigned __int128 u128;
typedef __int128 s128;
#define c3y 0
#define c3n 1

// ---- VERBATIM helpers from twoc.c ----
static inline void _tg_mask(mpz_t r, const mpz_t a, c3_d wid){ mpz_fdiv_r_2exp(r,a,wid); }
static inline c3_t _tg_sign(const mpz_t m, c3_d wid){ return mpz_tstbit(m,wid-1)?c3y:c3n; }
static inline void _tg_signed(mpz_t s, const mpz_t m, c3_d wid){
if ( _tg_sign(m,wid)==c3y ){ mpz_t p; mpz_init(p); mpz_ui_pow_ui(p,2,wid);
mpz_sub(s,m,p); mpz_clear(p);
} else mpz_set(s,m);
}
// ---- op blocks mirroring the jet GMPBLOCKs (operate on in-range inputs) ----
static void G_add(mpz_t r,mpz_t a,mpz_t b,c3_d w){ mpz_add(r,a,b); _tg_mask(r,r,w); }
static void G_sub(mpz_t r,mpz_t a,mpz_t b,c3_d w){ mpz_sub(r,a,b); _tg_mask(r,r,w); }
static void G_mul(mpz_t r,mpz_t a,mpz_t b,c3_d w){ mpz_mul(r,a,b); _tg_mask(r,r,w); }
static void G_neg(mpz_t r,mpz_t a,c3_d w){ mpz_neg(r,a); _tg_mask(r,r,w); }
static void G_abs(mpz_t r,mpz_t a,c3_d w){ mpz_t m; mpz_init(m); _tg_mask(m,a,w);
if(_tg_sign(m,w)==c3y){ mpz_neg(m,m); _tg_mask(m,m,w);} mpz_set(r,m); mpz_clear(m); }
static void G_div(mpz_t r,mpz_t a,mpz_t b,c3_d w){ mpz_t am,bm,as,bs;
mpz_inits(am,bm,as,bs,NULL); _tg_mask(am,a,w); _tg_mask(bm,b,w);
_tg_signed(as,am,w); _tg_signed(bs,bm,w); mpz_tdiv_q(as,as,bs); _tg_mask(as,as,w);
mpz_set(r,as); mpz_clears(am,bm,as,bs,NULL); }
static void G_rem(mpz_t r,mpz_t a,mpz_t b,c3_d w){ mpz_t am,bm,as,bs;
mpz_inits(am,bm,as,bs,NULL); _tg_mask(am,a,w); _tg_mask(bm,b,w);
_tg_signed(as,am,w); _tg_signed(bs,bm,w); mpz_tdiv_r(as,as,bs); _tg_mask(as,as,w);
mpz_set(r,as); mpz_clears(am,bm,as,bs,NULL); }
static void G_pow(mpz_t r,mpz_t a,c3_d n,c3_d w){ mpz_t am,nm,mod;
mpz_inits(am,nm,mod,NULL); _tg_mask(am,a,w); mpz_set_ui(nm,n);
mpz_ui_pow_ui(mod,2,w); mpz_powm(am,am,nm,mod); mpz_set(r,am);
mpz_clears(am,nm,mod,NULL); }
static c3_t G_gth(mpz_t a,mpz_t b,c3_d w){ mpz_t as,bs; mpz_inits(as,bs,NULL);
_tg_signed(as,a,w); _tg_signed(bs,b,w); int c=mpz_cmp(as,bs);
mpz_clears(as,bs,NULL); return c>0?c3y:c3n; }

// ---- independent hardware reference (wid <= 100) ----
static u128 R_mask(int w){ return (((u128)1)<<w)-1; }
static s128 R_sgn(u128 p,int w){ p&=R_mask(w);
if((p>>(w-1))&1) return (s128)p-(s128)(((u128)1)<<w); return (s128)p; }
static u128 R_pat(s128 v,int w){ return ((u128)v)&R_mask(w); }

static void u128_to_mpz(mpz_t r,u128 v){ mpz_set_ui(r,(c3_d)(v>>64));
mpz_mul_2exp(r,r,64); mpz_t lo; mpz_init_set_ui(lo,(c3_d)v); mpz_add(r,r,lo); mpz_clear(lo); }
static u128 mpz_to_u128(const mpz_t v){ mpz_t t,lo; mpz_inits(t,lo,NULL);
mpz_fdiv_r_2exp(lo,v,64); mpz_fdiv_q_2exp(t,v,64);
u128 r=((u128)mpz_get_ui(t)<<64)|(u128)mpz_get_ui(lo); mpz_clears(t,lo,NULL); return r; }

static long long FAILS=0, TOTS=0;
static void chku(const char*op,int w,u128 got,u128 ref){ TOTS++;
if(got!=ref){ FAILS++; if(FAILS<=40)
fprintf(stderr,"FAIL %s w=%d got=%llx%016llx ref=%llx%016llx\n",op,w,
(c3_d)(got>>64),(c3_d)got,(c3_d)(ref>>64),(c3_d)ref); } }

int main(void){
int widths[]={8,16,17,32,64,100};
mpz_t a,b,r; mpz_inits(a,b,r,NULL);
for(int wi=0;wi<6;wi++){ int w=widths[wi]; u128 m=R_mask(w), half=((u128)1)<<(w-1);
u128 vals[]={0,1,2,3,half-1,half,half+1,m-1,m,(u128)5<<(w/2),m/3,(m/7)|half};
int NV=12;
for(int i=0;i<NV;i++)for(int j=0;j<NV;j++){
u128 ua=vals[i]&m, ub=vals[j]&m; u128_to_mpz(a,ua); u128_to_mpz(b,ub);
G_add(r,a,b,w); chku("add",w,mpz_to_u128(r),(ua+ub)&m);
G_sub(r,a,b,w); chku("sub",w,mpz_to_u128(r),R_pat(R_sgn(ua,w)-R_sgn(ub,w),w));
G_mul(r,a,b,w); chku("mul",w,mpz_to_u128(r),(ua*ub)&m);
chku("gth",w,G_gth(a,b,w),(R_sgn(ua,w)>R_sgn(ub,w))?c3y:c3n);
if((ub&m)!=0){ s128 sa=R_sgn(ua,w),sb=R_sgn(ub,w); s128 q=sa/sb, rr=sa-q*sb;
G_div(r,a,b,w); chku("div",w,mpz_to_u128(r),R_pat(q,w));
G_rem(r,a,b,w); chku("rem",w,mpz_to_u128(r),R_pat(rr,w)); }
}
for(int i=0;i<NV;i++){ u128 ua=vals[i]&m; u128_to_mpz(a,ua);
G_neg(r,a,w); chku("neg",w,mpz_to_u128(r),R_pat(-R_sgn(ua,w),w));
G_abs(r,a,w); chku("abs",w,mpz_to_u128(r),(R_sgn(ua,w)<0)?R_pat(-R_sgn(ua,w),w):(ua&m));
c3_d ns[]={0,1,2,3,5,8,13};
for(int k=0;k<7;k++){ c3_d n=ns[k]; u128 acc=1&m,base=ua&m;
for(c3_d t=0;t<n;t++) acc=(acc*base)&m;
G_pow(r,a,n,w); chku("pow",w,mpz_to_u128(r),acc); }
}
}
// hand-checked wide cases (w=200): -1 + -1 = -2; -1 * -1 = 1; min/-1 = min; max+1 = min
{ c3_d w=200; mpz_t one,negone,res,exp; mpz_inits(one,negone,res,exp,NULL);
mpz_set_ui(one,1); _tg_mask(negone,/*=-1 as 2^w-1*/ negone,w);
mpz_ui_pow_ui(negone,2,w); mpz_sub_ui(negone,negone,1); // 2^200 - 1 == -1
G_add(res,negone,negone,w); mpz_ui_pow_ui(exp,2,w); mpz_sub_ui(exp,exp,2);
TOTS++; if(mpz_cmp(res,exp)!=0){FAILS++; fprintf(stderr,"FAIL wide add -1+-1\n");}
G_mul(res,negone,negone,w); TOTS++; if(mpz_cmp_ui(res,1)!=0){FAILS++; fprintf(stderr,"FAIL wide mul -1*-1\n");}
mpz_t mn; mpz_init(mn); mpz_ui_pow_ui(mn,2,w-1); // min = 2^199
G_div(res,mn,negone,w); TOTS++; if(mpz_cmp(res,mn)!=0){FAILS++; fprintf(stderr,"FAIL wide div min/-1\n");}
mpz_t mx; mpz_init(mx); mpz_ui_pow_ui(mx,2,w-1); mpz_sub_ui(mx,mx,1); // max=2^199-1
G_add(res,mx,one,w); TOTS++; if(mpz_cmp(res,mn)!=0){FAILS++; fprintf(stderr,"FAIL wide max+1\n");}
mpz_clears(one,negone,res,exp,mn,mx,NULL);
}
mpz_clears(a,b,r,NULL);
printf("checks=%lld fails=%lld\n",TOTS,FAILS);
return FAILS?1:0;
}
Loading