ICS Lab 攻略

Lab 1: 大整数运算

可移植性

PRIu64

CSAPP Page 47

main.c 中 test 函数所使用的宏 PRIu64：

void test(uint64_t a, uint64_t b, uint64_t m) {
  #define U64 "%" PRIu64
  printf(U64 " * " U64 " mod " U64 " = " U64 "\n", a, b, m, multimod(a, b, m));
}

64 位会被展开为 lu，使用 gcc -E main.c | indent - | vim - 观察：

#5 "main.c"
     uint64_t multimod (uint64_t, uint64_t, uint64_t);

     void test (uint64_t a, uint64_t b, uint64_t m)
{

  printf ("%"
#9 "main.c" 3 4
      "l" "u"
#9 "main.c"
      " * " "%"
#9 "main.c" 3 4
      "l" "u"
#9 "main.c"
      " mod " "%"
#9 "main.c" 3 4
      "l" "u"
#9 "main.c"
      " = " "%"
#9 "main.c" 3 4
      "l" "u"
#9 "main.c"
      "\n", a, b, m, multimod (a, b, m));
}

int
main ()
{
  test (123, 456, 789);
  test (123, 456, -1ULL);
  test (-2ULL, -2ULL, -1ULL);
}

32 位会被展开为 llu，使用 gcc -E -m32 main.c | indent - | vim - 观察：

#5 "main.c"
     uint64_t multimod (uint64_t, uint64_t, uint64_t);

     void test (uint64_t a, uint64_t b, uint64_t m)
{

  printf ("%"
#9 "main.c" 3 4
      "ll" "u"
#9 "main.c"
      " * " "%"
#9 "main.c" 3 4
      "ll" "u"
#9 "main.c"
      " mod " "%"
#9 "main.c" 3 4
      "ll" "u"
#9 "main.c"
      " = " "%"
#9 "main.c" 3 4
      "ll" "u"
#9 "main.c"
      "\n", a, b, m, multimod (a, b, m));
}

int
main ()
{
  test (123, 456, 789);
  test (123, 456, -1ULL);
  test (-2ULL, -2ULL, -1ULL);
}

总结来说：

在 32 位平台

typedef unsigned long long int  uint64_t;

在 64 位平台

typedef unsigned long int  uint64_t;

ULL

常量后缀，代表 unsigned long long int，一定为 64-bit

参考实现与测试

使用 Python 进行，文件名为 cheat.py：

#!/usr/bin/python3

import sys

assert len(sys.argv) == 4

a = int(sys.argv[1])
b = int(sys.argv[2])
m = int(sys.argv[3])

assert 1 <= a < pow(2,64)
assert 1 <= b < pow(2,64)
assert 1 <= m < pow(2,64)

print(a*b%m)

注意修改 cheat.py 文件的可执行性。

之后在 main 函数中添加 cheat 函数：

char* cheat(uint64_t a, uint64_t b, uint64_t m, char *buf) {
    char cmd[200];
    sprintf(cmd, "./cheat.py " U64 " " U64 " " U64, a, b, m);
    FILE *fp = popen(cmd, "r");
    assert(fp);
    fscanf(fp, "%s", buf);
    //printf("popen() returns: %s\n", buf);
    pclose(fp);
    return buf;
}

test 函数调用了 cheat 函数，通过字符串比较的方式判断正确性：

void test(uint64_t a, uint64_t b, uint64_t m) {
    //printf(U64 " * " U64 " mod " U64 " = " U64 "\n", a, b, m, multimod(a, b, m));
    char buf[100];
    char real[100];
    strcpy(buf, cheat(a, b, m, buf));
    uint64_t res = multimod(a, b, m);
    sscanf(real, U64, &res);
    assert(strcmp(buf, real));
}

这里的 U64 定义如下：

#define U64 "%" PRIu64

最后在 main 函数中使用随机生成的 64 位无符号整数进行测试：

uint64_t rand_uint64_slow() {
    uint64_t r = 0;
    for (int i=0;i<64;++i)
        r=r*2+rand()%2;
    if(r!=0)
        return r;
    else
        return 1;
}


int main(int argc, char *argv[]) {
    assert(argc>=2);
    char *end;
    long count = strtol(argv[1],&end,10);
    srand(time(0));
    for(long i=0;i<count;++i){
        uint64_t a = rand_uint64_slow();
        uint64_t b = rand_uint64_slow();
        uint64_t c = rand_uint64_slow();
        test(a, b, c);
    }
    printf("AC!\n");
}

实际的实现

参考机组，即无符号整数的乘法和除法：

#include <stdint.h>
#include <assert.h>

int mulRes[256];
int modRes[64];

int innerAdd64(int *a, int *b) { // a = a + b (64-bit)
    int carry = 0;
    for(int i=63;i>=0;--i){
        int res = a[i] ^ b[i] ^ carry;
        carry = (a[i] & carry) | (b[i] & carry) | (a[i] & b[i]);
        a[i] = res;
    }
    return carry;
}

int innerAdd128(int *a, int *b) { // a = a + b (128-bit)
    int carry = 0;
    for(int i=127;i>=0;--i){
        int res = a[i] ^ b[i] ^ carry;
        carry = (a[i] & carry) | (b[i] & carry) | (a[i] & b[i]);
        a[i] = res;
    }
    return carry;
}


void multi(uint64_t a, uint64_t b) {
    for(int i=64;i<128;++i)
        mulRes[i]=b>>(127-i)&1;
    int temp[64];
    for(int i=0;i<64;++i){
        // prepare temp
        if(mulRes[127]==1){
            for(int j=0;j<64;++j)
                temp[j]=a>>(63-j)&1;
        }else{
            for(int j=0;j<64;++j)
                temp[j]=0;
        }

        // modify mulRes
        int carry = innerAdd64(mulRes, temp); // carry!

        // shift right
        for(int i=127;i>=1;--i)
            mulRes[i] = mulRes[i-1];
        mulRes[0]=carry;
    }
}

int compare(int *a, int *b) { // 128-bit
    for(int i=0;i<128;++i){
        if(a[i]<b[i])
            return 0; // a < b
        else if(a[i]>b[i])
            return 1; // a > b
    }
    return 1; // a == b
}


void div(uint64_t m) {
    // prepare
    for(int i=0;i<128;++i)
        mulRes[i+128]=mulRes[i];
    for(int i=0;i<128;++i)
        mulRes[i]=0;

    // ori m
    int ori[128];
    for(int i=0;i<64;++i)
        ori[i+64]=m>>(63-i)&1;
    for(int i=0;i<64;++i){
        ori[i]=0;
    }

    // negate m
    int index=63;
    for(;index>=0;--index)
        if((m>>(63-index)&1)==1)
            break;
    assert(index>=0);
    int negate[128];
    for(int i=64+index-1;i>=64+0;--i)
        negate[i]=!(m>>(64+63-i)&1);
    for(int i=64+63;i>=64+index;--i)
        negate[i]=(m>>(64+63-i)&1);
    for(int i=0;i<64;++i)
        negate[i]=1;

    // iterate
    for(int i=0;i<128;++i){
        if(compare(mulRes,ori)){
            innerAdd128(mulRes,negate);
            // shift left
            for(int j=0;j<255;++j)
                mulRes[j]=mulRes[j+1];
            mulRes[255]=1;
        }
        else{
            for(int j=0;j<255;++j)
                mulRes[j]=mulRes[j+1];
            mulRes[255]=0;
        }
    }

    // last
    if(compare(mulRes,ori))
        innerAdd128(mulRes,negate);
    for(int i=0;i<64;++i)
        modRes[i]=mulRes[i+64];

}


uint64_t multimod(uint64_t a, uint64_t b, uint64_t m) {
    for(int i=0;i<64;++i)
        modRes[i]=0;
    for(int i=0;i<256;++i)
        mulRes[i]=0;
    assert(m != 0);
    multi(a, b);
    div(m);
    uint64_t res = 0;
    for(int i=0;i<64;++i){
        if(modRes[i]==1)
            res+=1;
        if(i!=63)
            res*=2;
    }
    return res;
}

无符号整数的乘法注意在右移时需要考虑进位。

无符号整数的除法则是 128 位除法，最后取商的低 64 位作为结果。

submit

在完成所有工作后，在 ~/ics-workbench 进行 git add . 和 git commit。

然后修改 ~/ics-workbench/multimod 中的 Makefile 文件，添加学号和姓名：

NAME := $(shell basename $(PWD))
all: $(NAME)-64 $(NAME)-32
export MODULE := Lab1
export STUID = XXX
export STUNAME = XXX
include ../Makefile

最后在 ~/ics-workbench/multimod 中 make submit 即可。

卷积

TODO

一个部分正确的 Ο(1) 实现

ICS EX

Lab 2: x86-64 内联汇编

终于等到了……

Ref.

How to Use Inline Assembly Language in C Code

https://www.ibiblio.org/gferg/ldp/GCC-Inline-Assembly-HOWTO.html

GDB

工欲善其事，必先利其器

https://zhuanlan.zhihu.com/p/32843449

经过一番体验，最终选择了 gdbgui。

使用 pip 安装完成后，就像使用 gdb 一样键入：

vgalaxy@vgalaxy-VirtualBox:~/Templates$ gdbgui ./a.out

便会弹出浏览器窗口，有时候需要刷新。

下面是一些快捷键：

Run: r
Continue: c
Next: n or right arrow
Step: s or down arrow
Up: u or up arrow
Next Instruction: m
Step Instruction: ,

对寄存器的支持似乎不太行，可以使用 expressions 窗口观察

报错 OSError: [Errno 98] Address already in use 后的处理

sudo lsof -i:5000
sudo kill -9 [process_id]

实验要求

在本实验中，所有函数的所有部分都必须使用 inline assembly 实现。除定义临时变量 (可以赋常量初值) 和 return 返回临时变量/参数外，不得使用任何表达式/条件/循环等 C 语言成分。

以 64 位为例，编译优化级别应使得参数传递为 x86-64 风格（寄存器）

一些调试友好的选项：

vgalaxy@vgalaxy-VirtualBox:~/ics-workbench/asm$ gcc -g -Og -Wall -Werror main.c asm-impl.c

asm_add

int64_t asm_add(int64_t a, int64_t b) {
    int64_t res = 0;
    asm (
            "addq %1, %2; "
            "movq %2, %0; "
            : "=r"(res)
            : "r"(a), "r"(b)
            :
        );
    return res;
}

先简单熟悉一下。

asm_popcnt

int asm_popcnt(uint64_t x) {
  int res = 0;
  asm (
          "mov $0x0, %%ecx; "
          "mov $0x0, %%edx; "
          ".byte 0xeb, 0x03; " // jmp
          "add $0x1, %%ecx; "
          "cmp $0x3f, %%ecx; "
          ".byte 0x7f, 0x0f; " // jg
          "movq %1, %%rax; "
          "shr %%cl, %%rax; "
          "test $0x1, %%al; "
          ".byte 0x74, 0xee; " // je
          "add $0x1, %%edx; "
          ".byte 0xeb, 0xe9; " // jmp
          "movl %%edx, %0; "
          : "=r"(res)
          : "r"(x)
          : "ecx", "edx", "rax"
      );
  return res;
}

其中寄存器 ecx 作为循环变量，寄存器 edx 保存中间结果。

由于跳转指令为相对跳转，故根据指令长度显式编码。

注意需要有显式的 return 语句，否则 OJ 会报错 control reaches end of non-void function。

asm_memcpy

先看源代码反汇编的结果：

0000000000001070 <memcpy@plt>:
    1070:  f3 0f 1e fa            endbr64
    1074:  f2 ff 25 55 2f 00 00   bnd jmp *0x2f55(%rip)        # 3fd0 <memcpy@GLIBC_2.14>
    107b:  0f 1f 44 00 00         nopl   0x0(%rax,%rax,1)

0000000000001169 <asm_memcpy>:
    1169:  f3 0f 1e fa            endbr64
    116d:  48 83 ec 08            sub    $0x8,%rsp
    1171:  e8 fa fe ff ff         call   1070 <memcpy@plt>
    1176:  48 83 c4 08            add    $0x8,%rsp
    117a:  c3                     ret

似乎直接调用 memcpy@plt 有亿点困难。

所以参考之前的实现：

void *memcpy(void *s1, const void *s2, size_t n) {
  unsigned char *us1=(unsigned char *)s1;
  const unsigned char *us2=(const unsigned char *)s2;
  for(;n>0;++us1,++us2,--n)
      *us1=*us2;
  return s1;
}

给出对应的内联汇编：

void *asm_memcpy(void *dest, const void *src, size_t n) {
    void *res = NULL;
    asm (
            "movq %1, %%rax; "
            ".byte 0xeb, 0x11; " // jmp
            "movzbl (%2), %%ecx; "
            "movb %%cl, (%%rax); "
            "addq $0x1, %%rax; "
            "addq $0x1, %2; "
            "subq $0x1, %3; "
            "test %3, %3; "
            ".byte 0x75, 0xea; " // jne
            "movq %1, %0; "
            : "=r"(res)
            : "r"(dest), "r"(src), "r"(n)
            : "rax", "ecx"
        );
    return res;
}

一个小疑问，对参数的写为什么不需要在 output operands 中指出？

asm_setjmp / asm_longjmp

Ref.

Intro & i386: https://zhuanlan.zhihu.com/p/82492121

x86-64: https://blog.csdn.net/dog250/article/details/89742140

struct: https://blog.codingnow.com/2005/12/typedef_struct_array.html

实现及原理

先看函数签名：

int     asm_setjmp(asm_jmp_buf env);
void    asm_longjmp(asm_jmp_buf env, int val);

这里的 asm_jmp_buf 是一个自定义数据结构：

typedef struct {
    uint64_t rsp;
    uint64_t rip;
    uint64_t rbp;
    uint64_t rbx;
    uint64_t r12;
    uint64_t r13;
    uint64_t r14;
    uint64_t r15;
} context[1];

#define asm_jmp_buf context

有两处需要注意的地方：

需要保存的环境

STFM

The  setjmp()  function  saves various information about the calling environment (typically, the stack pointer, the instruction pointer, possibly the  values  of other  registers  and  the  signal  mask)  in  the  buffer  env for later use by longjmp().  In this case, setjmp() returns 0.

可知保存寄存器 rsp 和 rip 就可以了，剩下的 6 个是被调用者保存寄存器（虽然讲不出什么道理）。

结构体的定义

在使用函数时，类似这样：

asm_jmp_buf buf;
int r = asm_setjmp(buf);

参数是直接传入到函数中的，若将 context[1] 改为 context，便成了值传递，不满足需求。使用 context[1]，可以在函数调用时让寄存器 rdi 保存 context 数组的地址，这样就可以通过寄存器 rdi 来保存环境了。

下面看实现：

int asm_setjmp(asm_jmp_buf env) {
    int res = 0;
    asm volatile (
            "movq %%rsp, (%1); "
            "movq (%%rsp), %%rdx; "
            "movq %%rdx, (8)(%1); " // rip: return addr
            "movq %%rbp, (16)(%1); "
            "movq %%rbx, (24)(%1); "
            "movq %%r12, (32)(%1); "
            "movq %%r13, (40)(%1); "
            "movq %%r14, (48)(%1); "
            "movq %%r15, (56)(%1); "
            "movl $0x0, %0; "
            : "=r"(res)
            : "r"(env)
            : "rdx", "rax"
        );
    return res;
}

void asm_longjmp(asm_jmp_buf env, int val) {
    asm volatile (
            "movl %1, %%eax; "
            "test %1, %1; "
            ".byte 0x75, 0x4; " // jne
            "add $0x1, %%rax; "
            "movq (%0), %%rsp; "
            "movq (8)(%0), %%rdx; "
            "movq %%rdx, (%%rsp); "
            "movq (16)(%0), %%rbp; "
            "movq (24)(%0), %%rbx; "
            "movq (32)(%0), %%r12; "
            "movq (40)(%0), %%r13; "
            "movq (48)(%0), %%r14; "
            "movq (56)(%0), %%r15; "
            :
            : "r"(env), "r"(val)
            : "rdx", "rax"
        );
}

注意 asm_longjmp 中要特判 val 是否为 0：

If the programmer mistakenly passes the value 0 in val, the "fake" return will instead return 1.

还写了一个很有感觉的测试：

#include "asm.h"
#include <assert.h>
#include <stdio.h>

const char src[] = "HELL...O";
char dest[10];

asm_jmp_buf env;

void f(int n) {
  if (n >= 18) asm_longjmp(env, n); // 某个条件达成时，恢复快照
  printf("Call f(%d)\n", n);
  f(n + 1);
}

void recursion() {
  int r = asm_setjmp(env); // 快照
  if (r == 0) {
    f(1);
  } else { // longjmp goes here
    printf("Recursion reaches %d\n", r);
  }
}

int main() {
  asm_jmp_buf buf;
  int r = asm_setjmp(buf);
  if (r == 0) {
    assert(asm_add(1234, 5678) == 6912);
    assert(asm_popcnt(0x0123456789abcdefULL) == 32);

    asm_memcpy(dest, src, 7);
    printf("%s\n", dest);

    asm_longjmp(buf, 123);
  } else {
    assert(r == 123);
    recursion();
    printf("PASSED.\n");
  }
}

递归的十八层地狱……

通过 gdbgui 来分析一下执行过程：

从 main 进入 asm_setjmp 后，返回地址为 0x555555555248，即 (%rsp)，此时 rsp 的值为 0x7fffffffddc8。我们保存寄存器 rsp 到环境，并将解引用的值存入环境中的 rip。返回值为 0。
从 main 进入 asm_longjmp 后，rsp 的值为 0x7fffffffddc8（与环境中的恰好相同，由于是同层），我们将环境中 rip 的值写 (%rsp)，从而改变了 main 为 asm_longjmp 准备的返回地址，于是 asm_longjmp 便返回到了 asm_setjmp 返回的地方。由于寄存器 rax 的值也被改写了，所以此时 r 的值为 123。
recursion 进行递归调用下的测试，环境中 rsp 为 0x7fffffffddb8，rip 为 0x5555555551fb。多次递归后，实际 rsp 的值在减少。最终到达 asm_longjmp 后，将 rsp 又改为 0x7fffffffddb8，相当于将低地址的栈全部 unwind 了，并改写返回地址，于是便成功进行了 nonlocal goto。

协程

使用 setjmp / longjmp 可以实现异常处理外，还可以实现协程：

https://www.cnblogs.com/sewain/p/14360853.html

#include <stdio.h>
#include <unistd.h>
#include <setjmp.h>

typedef int     BOOL;
#define TRUE    1
#define FALSE   0

// 用来存储主程和协程的上下文的数据结构
typedef struct _Context_ {
    jmp_buf mainBuf;
    jmp_buf coBuf;
} Context;

// 上下文全局变量
Context gCtx;

// 恢复
#define resume() \
    if (0 == setjmp(gCtx.mainBuf)) \
    { \
        longjmp(gCtx.coBuf, 1); \
    }

// 挂起
#define yield() \
    if (0 == setjmp(gCtx.coBuf)) \
    { \
        longjmp(gCtx.mainBuf, 1); \
    }

// 在协程中执行的函数
void coroutine_function(void *arg)
{
    while (TRUE)  // 死循环
    {
        printf("\n*** coroutine: working \n");
        // 模拟耗时操作
        for (int i = 0; i < 10; ++i)
        {
            fprintf(stderr, ".");
            usleep(1000 * 200);
        }
        printf("\n*** coroutine: suspend \n");

        // 让出 CPU
        yield();
    }
}

// 启动一个协程
// 参数1：func 在协程中执行的函数
// 参数2：func 需要的参数
typedef void (*pf)(void *);
BOOL start_coroutine(pf func, void *arg)
{
    // 保存主程的跳转点
    if (0 == setjmp(gCtx.mainBuf))
    {
        func(arg); // 调用函数
        return TRUE;
    }

    return FALSE;
}

int main()
{
    // 启动一个协程
    start_coroutine(coroutine_function, NULL);

    while (TRUE) // 死循环
    {
        printf("\n=== main: working \n");

        // 模拟耗时操作
        for (int i = 0; i < 10; ++i)
        {
            fprintf(stderr, ".");
            usleep(1000 * 200);
        }

        printf("\n=== main: suspend \n");

        // 放弃 CPU，让协程执行
        resume();
    }

    return 0;
}

编译选项如下：

gcc -g -Og -D_FORTIFY_SOURCE=0 a.c

对于其中的 -D_FORTIFY_SOURCE，参阅 https://stackoverflow.com/questions/13517526/difference-between-gcc-d-fortify-source-1-and-d-fortify-source-2

部分执行结果如下：

*** coroutine: working
..........
*** coroutine: suspend

=== main: working
..........
=== main: suspend

*** coroutine: working
..........
*** coroutine: suspend

其原理并不难看懂。

Duff’s Device

起源于循环展开，却可以来实现协程……

https://zhuanlan.zhihu.com/p/284223705

https://www.chiark.greenend.org.uk/~sgtatham/coroutines.html

https://mthli.xyz/coroutines-in-c/

Lab 3: 性能调优

Premature optimization is the root of all evil.

Inner Timer

在 main.c 中：

#include <stdint.h>
#include <stdio.h>

#include <sys/time.h>

int *sieve(int n);

int main() {
  struct timeval t1,t2;
  gettimeofday(&t1, NULL);
  int *primes = sieve(9000000);
  gettimeofday(&t2, NULL);
  for (int *cur = primes, i = 0; *cur; cur++, i++) {
    printf("%8d", *cur);
    if (i % 8 == 7 || !*(cur + 1)) printf("\n");
  }
  uint64_t use = (t2.tv_sec - t1.tv_sec) * 1000000 + (t2.tv_usec - t1.tv_usec);
  printf("inner timer: %lu\n", use);
}

以微秒为单位。

Outer Timer

基本的计时功能：

$ time -p ./a.out

Trace

系统调用计时和库函数调用记录：

$ strace -T ./a.out
$ ltrace python3 --help

Profiler

每隔一段时间对程序的执行采样。

例如如果我们希望了解哪个函数使用了最多的时间，可以想到如下的实现：

借助系统的机制，每秒给进程发送若干 (例如 100) 次“中断”，中断到来后，会跳转到一个“中断”处理程序。
在“中断”处理程序中，我们遍历函数的调用栈，获得当前的调用链并记录。
程序执行完毕后，根据采样的得到的信息，就能推断出哪些函数 (近似) 消耗了多少时间。

Linux 为我们提供了 perf tools 系列的性能诊断工具，安装如下：

$ sudo apt install linux-tools-common
$ sudo apt install linux-tools-5.11.0-16-generic
$ sudo apt install linux-cloud-tools-5.11.0-16-generic

然后以管理员权限运行：

$ sudo perf record ./a.out

此时在当前目录下会生成 perf.data 文件，我们键入：

$ sudo perf report

就可以进行一番探索了……

此外，可以用 perf stat 来查看一次运行的统计信息……

For more:

RTFM

https://www.brendangregg.com/perf.html

Makefile

.PHONY: time strace perf-report perf-stat perf-stat-detail remove

time: $(NAME)-64
    time -p ./$(NAME)-64

strace: $(NAME)-64
    strace -T ./$(NAME)-64

perf-report: $(NAME)-64
    sudo perf record ./$(NAME)-64
    sudo perf report

perf-stat: $(NAME)-64
    sudo perf stat ./$(NAME)-64

perf-stat-detail: $(NAME)-64
    sudo perf stat -e cycles,instructions,L1-dcache-loads,L1-dcache-load-misses,LLC-loads,LLC-load-misses ./$(NAME)-64

remove:
    rm -f ./$(NAME)-64
    sudo rm -f perf.data*

具体的优化

Set N = 9000000

原始

#include <stdbool.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>

#define N 10000000

static bool is_prime[N];
static int  primes[N];

int *sieve(int n) {
  assert(n + 1 < N);
  for (int i = 0; i <= n; i++)
    is_prime[i] = true;

  for (int i = 2; i <= n; i++) {
    for (int j = i + i; j <= n; j += i) {
      is_prime[j] = false;
    }
  }

  int *p = primes;
  for (int i = 2; i <= n; i++)
    if (is_prime[i]) {
      *p++ = i;
    }
  *p = 0;
  return primes;
}

perf stat 查看运行时间：

       1.192843152 seconds time elapsed

       0.332289000 seconds user
       0.132915000 seconds sys

筛至平方根

小优化：

       1.090751154 seconds time elapsed

       0.183597000 seconds user
       0.154388000 seconds sys

欧拉筛

#include <stdbool.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>

#define N 10000000

static bool is_prime[N] = {1, 1};
static int  primes[N];

int *sieve(int n) {
  assert(n + 1 < N);

  int primes_index = 0;

  for (int i = 2; i <= n; i++) {
    if (!is_prime[i]) primes[primes_index++] = i;
    for (int j = 0; j < primes_index; ++j) {
      if (i * primes[j] > N) break;
      is_prime[i * primes[j]] = true;
      if (i % primes[j] == 0) break;
    }
  }

  primes[primes_index] = 0;
  return primes;
}

运行时间似乎并不乐观：

       1.035647067 seconds time elapsed

       0.084943000 seconds user
       0.169887000 seconds sys

只筛奇数

#include <stdbool.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>

#define N 10000000

static bool is_prime[N] = {1, 1, 0};
static int  primes[N] = {2};

int *sieve(int n) {
  assert(n + 1 < N);

  int primes_index = 1;

  for (int i = 3; i <= n; i += 2) {
    if (!is_prime[i]) primes[primes_index++] = i;
    for (int j = 0; j < primes_index; ++j) {
      if (i * primes[j] > N) break;
      is_prime[i * primes[j]] = true;
      if (!(i % primes[j])) break;
    }
  }

  primes[primes_index] = 0;
  return primes;
}

还是差一点：

       0.978246240 seconds time elapsed

       0.065922000 seconds user
       0.174501000 seconds sys

放弃了……

ICS Lab 攻略

TOC

ICS Lab 攻略

Lab 1: 大整数运算

可移植性

PRIu64

ULL

参考实现与测试

实际的实现

submit

卷积

一个部分正确的 Ο(1) 实现

Lab 2: x86-64 内联汇编

Ref.

GDB

实验要求

asm_add

asm_popcnt

asm_memcpy

asm_setjmp / asm_longjmp

实现及原理

协程

Duff’s Device

Lab 3: 性能调优

Inner Timer

Outer Timer

Trace

Profiler

Makefile

具体的优化

原始

筛至平方根

欧拉筛

只筛奇数