開平法

// 開平法

int number = 2; // ルートを求めたい整数。現状1から99までの整数
int digit = 100; // 求めたい有効桁数

void setup(){
  
  size(1000,1000,P2D);

  int[] root_n = new int[digit]; // root_n[0]は1から9までの整数、その他は0から9までの整数。0番目が最大桁の値。
  int[] left_num = new int[digit+1]; // 開平法の左側に出てくる数。各要素は0から9までの整数。0番目が最小桁の値。
  int[] right_num = new int[2*digit+1]; // 開平法の右側に出てくる数（差の値に次の桁を加えたもの）。各要素は0から9までの整数。0番目が最小桁の値。
  int[] right_num2 = new int[2*digit+1]; // 開平法の右側に出てくる数（left_num × left_num[0]）。各要素は0から9までの整数。0番目が最小桁の値。
  int[] right_num3 = new int[2*digit+1]; // 開平法の右側に出てくる数（差の値）。各要素は0から9までの整数。0番目が最小桁の値。
  int[] forComp = new int[2*digit+1]; // right_num3と比較する値。各要素は0から9までの整数。0番目が最小桁の値。
  int left_digit = 1;
  int right_digit = 1;
  int right_digit2 = 0;
  int right_digit3 = 0;
  int forComp_digit = 0;
  
  // left_num, right_numを初期化
  for(int i=0; i<digit+1; i++){
    left_num[i] = 0;
  }
  for(int i=0; i<2*digit+1; i++){
    right_num[i] = 0;
    right_num2[i] = 0;
    right_num3[i] = 0;
    forComp[i] = 0;
  }
  
  // rootの整数部分を算出
  for(int i=1; i<10; i++){
    if( (i+1) * (i+1) > number ){
      root_n[0] = i;
      break;
    }
  }
  
  // left_numを更新
  left_num[0] = 2 * root_n[0];
  if( left_num[0] >= 10 ){ 
    left_num[1] = left_num[0] / 10;
    left_num[0] = left_num[0] % 10;
    left_digit = 2;
  }  
  // right_numを更新
  right_num[0] = number - root_n[0] * root_n[0];
  if( right_num[0] >= 10 ){ 
    right_num[1] = right_num[0] / 10;
    right_num[0] = right_num[0] % 10;
    right_digit = 2;
  }
  
  for(int k=1; k<digit; k++){
      
    // left_numの配列を1桁ずつずらす
    for(int i=0; i<left_digit; i++){
      left_num[left_digit-i] = left_num[left_digit-1-i];
    }
    left_digit++;
  
    // right_numの配列を2桁ずつずらす
    for(int i=0; i<right_digit; i++){
      right_num[right_digit+1-i] = right_num[right_digit-1-i];
    }
    right_num[0] = 0;
    right_num[1] = 0;
    right_digit += 2;
  
    // root_nのk桁目の値を計算する
    left_num[0] = 1; // 0ではなく1を代入　比較は

    if( left_digit > right_digit ){
      root_n[k] = 0;
      continue;
    }
    boolean flag = false;
    if( left_digit == right_digit ){
      for(int i=0; i<right_digit; i++){
        if( left_num[left_digit-1-i] > right_num[right_digit-1-i] ){
          flag = true;
          break;
        } else if( left_num[left_digit-1-i] < right_num[right_digit-1-i] ){
          break;
        }
        
      }
      if( flag ){
        root_n[k] = 0;
        left_num[0] = 0;
        continue;
      }
    }
    
    for(int i=1; i<10; i++){
      left_num[0] = i;
    
      // left_num × left_num[0] = right_num2
      int carry = 0;
      for(int j=0; j<left_digit; j++){
        right_num2[j] = left_num[j] * left_num[0] + carry;
        carry = right_num2[j] / 10;
        right_num2[j] = right_num2[j] % 10;
      }
      right_digit2 = left_digit;
      if(carry > 0){
        right_digit2++;
        right_num2[right_digit2-1] = carry;
      }
    
      // right_num - right_num2 = right_num3
      int borrow = 0;
      for(int j=0; j<right_digit; j++){
        if( right_num[j] - borrow < right_num2[j] ){
          right_num3[j] = 10 + right_num[j] - borrow - right_num2[j];
          borrow = 1;
        } else {
          right_num3[j] = right_num[j] - borrow - right_num2[j];
          borrow = 0;
        }
      }
      right_digit3 = right_digit;
      for(int j=0; j<right_digit; j++){
        if( right_num3[right_digit-1-j] > 0 ){
          break;
        } else {
          right_digit3--;
        }
      }
      
      // 比較用の値を算出
      forComp_digit = left_digit;
      for(int j=0; j<left_digit; j++){
        forComp[j] = left_num[j];
      }
      forComp[0] = 2 * forComp[0] + 1;
      carry  = 0;
      for(int j=0; j<forComp_digit; j++){
        int temp = forComp[j] + carry;
        forComp[j] = temp % 10;
        carry = temp / 10;
      }
      if( carry > 0 ){
        forComp_digit++;
        forComp[forComp_digit-1] = carry;
      }     
      
      flag = false;
      if( forComp_digit > right_digit3 ){
        flag = true;
      }
      if( forComp_digit == right_digit3 ){
        for(int j=0; j<right_digit3; j++){
          if( forComp[forComp_digit-1-j] > right_num3[right_digit3-1-j] ){
            flag = true;
            break;
          } else if( forComp[forComp_digit-1-j] < right_num3[right_digit3-1-j] ){
            break;
          }
        }
      }
      if( flag ) {
        // rootのk桁目の値を確定
        root_n[k] = i;
        // right_numにright_num3を代入
        right_digit = right_digit3;
        for(int j=0; j<right_digit3; j++){
          right_num[j] = right_num3[j];
        }
        // left_numを更新
        left_num[0] = 2 * left_num[0];
        carry  = 0;
        for(int j=0; j<left_digit; j++){
          int temp = left_num[j] + carry;
          left_num[j] = temp % 10;
          carry = temp / 10;
        }
        if( carry > 0 ){
          left_digit++;
          left_num[left_digit-1] = carry;
        }
        break;
      }
    }     
  }

  print("root_n ");
  for(int i=0; i<digit; i++){
    print(root_n[i]);
  }
  
  // 開平法の結果を使って、作図
  background(255,255,255);
  int diameter = width / digit;
  int shift = diameter / 2;
  strokeWeight(diameter);
  for(int i=0; i<digit; i++){
    for(int j=0; j<digit; j++){
      if(j<i){
        stroke(25 * root_n[i], 25 * root_n[i], 25 * root_n[i]);
      } else {
        stroke(25 * root_n[j], 25 * root_n[j], 25 * root_n[j]);
      }
      point(i * diameter + shift, j * diameter + shift);
    }
  }

  save("extractionOfSquareRoot.jpg");  
  
}