SJTU Xflops 一轮考核小记

#include <time.h>
#include <stdlib.h>
#include <stdio.h>

static int check_with(int N, int a, int b, int o);
static int cmp_int(const void* l, const void* r);

int _main() {
  const int T=1000;
  const int MAX=4096;
  srand(time(0));
  for(int i=0; i<T; i++) {
    #define R (rand()%MAX)
    int nums[4]={R,R,R,R};
    #undef R
    qsort(nums,4,4,cmp_int);

    const int N=nums[3];
    const int a=nums[0], b=nums[2]-a, o=nums[1];
    if (check_with(N,a,b,o)!=0)
      return -1;
    printf("%d/%d check passed!\n", i, T);
  }
  return 0;
}

static void print_u64_le(unsigned long x) {
  for(int i=0; i<16; i++) {
    char c = x&15;
    if(c<10) c+='0';
    else c+='a'-10;
    putchar(c);
    x>>=4;
  }
}

static void print_buf(unsigned long *b) {
  for(int i=-1; i<3; i++) {
    print_u64_le(b[i]);
    putchar(' ');
  }
}

static int check_with(int N, int a, int b, int o) {
  bit_vector_t *bv1 = bit_vector_new(N), *bv2=bit_vector_new(N), *bv3=bit_vector_new(N);

  for (int i=0; i<(bv1->bit_sz+7)/8; i++){
    bv1->buf[i] = rand();
  }
  int or=o%b;
  printf("Rotate %d bits on bit %d-%d in %d bit bit_vector (%d-%d-%d)\n", or,a,a+b,N, a, a+or, a+b);
  memcpy(bv2->buf, bv1->buf, (N+7)/8);
  memcpy(bv3->buf, bv1->buf, (N+7)/8);
  const size_t bit_right_real = o%b;
  reverse_the_bit_vector(bv1, a, a+b-1);
  reverse_the_bit_vector_slow(bv2, a, a+b-1);
  for(int i=0; i<N; i++) {
    bool b1=bit_vector_get(bv1,i), b2=bit_vector_get(bv2,i);
    // printf("check on bit %d(%x/%x)\n",i, bv2->buf[i/8], bv1->buf[i/8]);
    if(b1!=b2) {
      printf("check failed on bit %d(%d/%d)\n",i, bit_vector_get(bv1,i), bv1->buf[i/8]);
      unsigned long *b1 = (unsigned long*)((size_t)(bv1->buf+i/8)&(~7ul));
      unsigned long *b2 = (unsigned long*)((size_t)(bv2->buf+i/8)&(~7ul));
      unsigned long *b3 = (unsigned long*)((size_t)(bv3->buf+i/8)&(~7ul));
      int ii = (i&(~63))-64;
      printf("        %-8d         %-8d         %-8d         %-8d         \n", ii, ii+64, ii+128, ii+192);
      printf("got:    ");
      print_buf(b1);
      putchar('\n');
      printf("expect: ");
      print_buf(b2);
      putchar('\n');
      printf("origin: ");
      print_buf(b3);
      putchar('\n');
      return -1;
    }
  }
  bit_vector_free(bv1); bit_vector_free(bv2); bit_vector_free(bv3);
  return 0;
}

static int cmp_int(const void* l, const void* r) {
  int f = *((int*)l);
  int s = *((int*)r);
  if (f > s) return  1;
  if (f < s) return -1;
  return 0;
}

use core::slice;
type Float4 = wide::f32x4;

#[no_mangle]
pub extern "C" fn gradient_descent(points: *mut f32, n: u32, m: u32, eta: f32, params: &[f32;4]) {
    if m==0 {return;}
    let p = [4.*params[0]*eta, 3.*params[1]*eta, 2.*params[2]*eta, params[3]*eta];
    let n = n as usize;
    assert!(points as usize % 8 == 0); // alignment check
    // just assume 4 threads
    let slice_16 = unsafe{ slice::from_raw_parts_mut(points as *mut [Float4;4], n/16) };
    parallel_for(slice_16, 4, &|v| run_gd_16(v, m, &p));
    let slice_4 = unsafe{ slice::from_raw_parts_mut((points.add(n/16*16)) as *mut Float4, (n%16)/4) };
    parallel_for(slice_4, 4, &|v| run_gd_4(v, m, &p));
    let slice_1 = unsafe{ slice::from_raw_parts_mut((points.add(n/4*4)) as *mut f32, n%4) };
    parallel_for(slice_1, 4, &|v| run_gd(v, m, &p));
}

pub fn parallel_for<T: Send>(slice: &mut [T], n_threads: usize, f: &(dyn Fn(&mut T)+Sync)) {
    let chunk_size = (slice.len()+n_threads-1)/n_threads;
    if chunk_size==0 {return;}
    std::thread::scope(|s| {
        for ch in slice.chunks_mut(chunk_size) {
            s.spawn(move || { for i in ch { f(i); } });
        }
    });
}

pub fn run_gd(px: &mut f32, rounds: u32, params: &[f32;4]) {
    let [a,b,c,d] = *params;
    let mut x = *px;
    for _ in 0..rounds {
        let x2 = x*x;
        let ax = a*x;
        x -= (ax*x2+b*x2)+(c*x+d);
    }
    *px = x;
}

pub fn run_gd_4(px: &mut Float4, rounds: u32, params: &[f32;4]) {
    let [a, b, c, d] = params.map(|v| Float4::splat(v));
    let mut x = *px;
    for _ in 0..rounds {
        let x2 = x*x;
        let ax = a*x;
        x -= (ax*x2+b*x2)+(c*x+d);
    }
    *px = x;
}

pub fn run_gd_16(px: &mut [Float4;4], rounds: u32, params: &[f32;4]) {
    let [a, b, c, d] = params.map(|v| Float4::splat(v));
    let mut x = *px;
    for _ in 0..rounds {
        for i in 0..4 {
            let (x2, ax) = (x[i]*x[i], a*x[i]);
            x[i] -= (ax*x2+b*x2)+(c*x[i]+d);
        }
    }
    *px = x;
}