text:
\subsection{Integration}
This section decribes all the functions in the file ``ptwXY\_integration.c''.

\subsubsection{ptwXY\_f\_integrate}
This function returns the integral between two points using \highlight{interpolation}.
\setargumentNameLengths{interpolation}
\CallingC{nfu\_status ptwXY\_f\_integrate(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{ptwXY\_interpolation interpolation,}
    \addArgument{double x1,}
    \addArgument{double y1,}
    \addArgument{double x2,}
    \addArgument{double y2,}
    \addArgument{double *value );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to a \highlight{ptwXYPoints} object.}
    \argumentBox{interpolation}{The interpolation between the two points.}
    \argumentBox{x2}{The x-value of the lower point.}
    \argumentBox{y2}{The y-value of the lower point}
    \argumentBox{x2}{The x-value of the upper point.}
    \argumentBox{y2}{The y-value of the upper point}
    \argumentBox{value}{On return, the value of the integral.}
    \vskip 0.05 in \noindent

\subsubsection{ptwXY\_integrate}
This function returns the integral of \highlight{ptwXY} from \highlight{domainMin} to \highlight{domainMax}.
\setargumentNameLengths{domainMax}
\CallingC{nfu\_status ptwXY\_integrate(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double domainMin,}
    \addArgument{double domainMax,}
    \addArgument{double *integral);}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{domainMin}{The lower limit of integration.}
    \argumentBox{domainMax}{The upperlimit of integration.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int_{\rm domainMin}^{\rm domainMax} f(x) dx$.

\subsubsection{ptwXY\_integrateDomain}
This function returns the integral of \highlight{ptwXY} over its domain.
\setargumentNameLengths{integral}
\CallingC{nfu\_status ptwXY\_integrateDomain(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double *value );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int f(x) dx$ over the domain of \highlight{ptwXY}.

\subsubsection{ptwXY\_normalize}
This function multiplies each y-value of \highlight{ptwXY} by a constant so that its integral is then normalized to 1.
\setargumentNameLengths{ptwXY}
\CallingC{nfu\_status ptwXY\_normalize(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \vskip 0.05 in \noindent

\subsubsection{ptwXY\_integrateDomainWithWeight\_x}
This function returns the integral of \highlight{ptwXY} weighted by x over its domain.
\setargumentNameLengths{integral}
\CallingC{nfu\_status ptwXY\_integrateDomainWithWeight\_x(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double *integral );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int x f(x) dx$ over the domain of \highlight{ptwXY}.

\subsubsection{ptwXY\_integrateWithWeight\_x}
This function returns the integral of \highlight{ptwXY} weighted by x from domainMin to domainMax.
\setargumentNameLengths{domainMax}
\CallingC{nfu\_status ptwXY\_integrateWithWeight\_x(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double domainMin,}
    \addArgument{double domainMax,}
    \addArgument{double *integral );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{domainMin}{The lower limit of the integration.}
    \argumentBox{domainMax}{The upper limit of the integration.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int_{\rm domainMin}^{\rm domainMax} x f(x) dx$ over the domain of \highlight{ptwXY}.

\subsubsection{ptwXY\_integrateDomainWithWeight\_sqrt\_x}
This function returns the integral of \highlight{ptwXY} weighted by $\sqrt{ x }$ over its domain.
\setargumentNameLengths{integral}
\CallingC{nfu\_status ptwXY\_integrateDomainWithWeight\_sqrt\_x(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double *integral );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int \sqrt{ x } f(x) dx$ over the domain of \highlight{ptwXY}.

\subsubsection{ptwXY\_integrateWithWeight\_sqrt\_x}
This function returns the integral of \highlight{ptwXY} weighted by x from domainMin to domainMax.
\setargumentNameLengths{domainMin}
\CallingC{nfu\_status ptwXY\_integrateWithWeight\_sqrt\_x(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{double domainMin,}
    \addArgument{double domainMax,}
    \addArgument{double *integral );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{domainMin}{The lower limit of the integration.}
    \argumentBox{domainMax}{The upper limit of the integration.}
    \argumentBox{integral}{On return, the value of the integral.}
    \vskip 0.05 in \noindent
The return value is $\int_{\rm domainMin}^{\rm domainMax} \sqrt{ x } f(x) dx$ over the domain of \highlight{ptwXY}.

\subsubsection{ptwXY\_groupOneFunction}
This function integrates \highlight{ptwXY} between each pair of consecutive points in \highlight{groupBoundaries} and
returns each integral's value as an element of the returned \highlight{ptwXPoints}.
\setargumentNameLengths{groupBoundaries}
\CallingC{ptwXPoints *ptwXY\_groupOneFunction(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY,}
    \addArgument{ptwXPoints *groupBoundaries,}
    \addArgument{ptwXY\_group\_normType normType,}
    \addArgument{ptwXPoints *norm );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{groupBoundaries}{A list of x-values.}
    \argumentBox{normType}{The type of normalization to apply to integration.}
    \argumentBox{norm}{A list of normalizations to be applied when normType is \highlight{ptwXY\-\_group\-\_normType\-\_norm}.}
    \vskip 0.05 in \noindent
Let \highlight{groupBoundaries} contain $n$ x-values with $x_i < x_{i+1}$. The returned \highlight{ptwXPoints} will contain $n-1$
values $I_i$ such that
\begin{equation}
    I_i = { 1 \over n_i } \int_{x_i}^{x_{i+1}} f(x) dx
\end{equation}
where $n_i$ is determined by \highlight{normType} as,
\begin{description}
    \item[ptwXY\_group\_normType\_none:] $n_i = 1$.
    \item[ptwXY\_group\_normType\_dx:] $n_i = x_{i+1} - x_i$.
    \item[ptwXY\_group\_normType\_norm:] $n_i = $ the $(i-1)^{th}$ element of norm.
\end{description}

\subsubsection{ptwXY\_groupTwoFunctions}
This function integrates the product of \highlight{ptwXY1} and \highlight{ptwXY2} between each pair of consecutive points in \highlight{groupBoundaries} and
returns each integral's value as an element of the returned \highlight{ptwXPoints}.
\setargumentNameLengths{groupBoundaries}
\CallingC{ptwXPoints *ptwXY\_groupTwoFunctions(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY1,}
    \addArgument{ptwXYPoints *ptwXY2,}
    \addArgument{ptwXPoints *groupBoundaries,}
    \addArgument{ptwXY\_group\_normType normType,}
    \addArgument{ptwXPoints *norm );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY1}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{ptwXY2}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{groupBoundaries}{A list of x-values.}
    \argumentBox{normType}{The type of normalization to apply to integration.}
    \argumentBox{norm}{A list of normalizations to be applied when normType is \highlight{ptwXY\-\_group\-\_normType\-\_norm}.}
    \vskip 0.05 in \noindent
Let \highlight{groupBoundaries} contain $n$ x-values with $x_i < x_{i+1}$. The returned \highlight{ptwXPoints} will contain $n-1$
values $I_i$ such that
\begin{equation}
    I_i = { 1 \over n_i } \int_{x_i}^{x_{i+1}} f(x) \; g(x) dx
\end{equation}
where $n_i$ is determined by \highlight{normType} as,
\begin{description}
    \item[ptwXY\_group\_normType\_none:] $n_i = 1$.
    \item[ptwXY\_group\_normType\_dx:] $n_i = x_{i+1} - x_i$.
    \item[ptwXY\_group\_normType\_norm:] $n_i = $ the $(i-1)^{th}$ element of norm.
\end{description}

\subsubsection{ptwXY\_groupThreeFunctions}
This function integrates the product \highlight{ptwXY1}, \highlight{ptwXY2} and \highlight{ptwXY3} between each pair of
consecutive points in \highlight{groupBoundaries} and
returns each integral's value as an element of the returned \highlight{ptwXPoints}.
\setargumentNameLengths{groupBoundaries}
\CallingC{ptwXPoints *ptwXY\_groupThreeFunctions(}{statusMessageReporting *smr,
    \addArgument{ptwXYPoints *ptwXY1,}
    \addArgument{ptwXYPoints *ptwXY2,}
    \addArgument{ptwXYPoints *ptwXY3,}
    \addArgument{ptwXPoints *groupBoundaries,}
    \addArgument{ptwXY\_group\_normType normType,}
    \addArgument{ptwXPoints *norm );}}
    \argumentBox{smr}{The \highlight{statusMessageReporting} instance to record errors.}
    \argumentBox{ptwXY1}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{ptwXY2}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{ptwXY3}{A pointer to the \highlight{ptwXYPoints} object.}
    \argumentBox{groupBoundaries}{A list of x-values.}
    \argumentBox{normType}{The type of normalization to apply to integration.}
    \argumentBox{norm}{A list of normalizations to be applied when normType is \highlight{ptwXY\-\_group\-\_normType\-\_norm}.}
    \vskip 0.05 in \noindent
Let \highlight{groupBoundaries} contain $n$ x-values with $x_i < x_{i+1}$. The returned \highlight{ptwXPoints} will contain $n-1$
values $I_i$ such that
\begin{equation}
    I_i = { 1 \over n_i } \int_{x_i}^{x_{i+1}} f(x) g(x) h(x) dx
\end{equation}
where $n_i$ is determined by \highlight{normType} as,
\begin{description}
    \item[ptwXY\_group\_normType\_none:] $n_i = 1$.
    \item[ptwXY\_group\_normType\_dx:] $n_i = x_{i+1} - x_i$.
    \item[ptwXY\_group\_normType\_norm:] $n_i = $ the $(i-1)^{th}$ element of norm.
\end{description}


text:
% Part: methods
% Chapter: proofs
% Section: resources

\documentclass[../../../include/open-logic-section]{subfiles}

\begin{document}

\olfileid{mth}{prf}{res}

\section{Other Resources}

There are many books on how to do proofs in mathematics which may be
useful. Check out \emph{How to Read and do Proofs: An Introduction to
Mathematical Thought Processes} \citep{Solow2013} and \emph{How to
Prove It: A Structured Approach} \citep{Velleman2019} in particular.
The
\href{http://www.people.vcu.edu/~rhammack/BookOfProof/BookOfProof.pdf}{\emph{Book
of Proof}} \citep{Hammack2013} and
\href{https://scholarworks.gvsu.edu/books/7/}{\emph{Mathematical
Reasoning}} \citep{Sandstrum2019} are books on proof that are freely
available online.  Philosophers might find
\emph{More
Precisely: The Math you need to do Philosophy} \citep{Steinhart2018} to
be a good primer on mathematical reasoning.

There are also various shorter guides to proofs available on the
internet; e.g.,
\href{https://math.berkeley.edu/~hutching/teach/proofs.pdf}{``Introduction
  to Mathematical Arguments''} \citep{Hutchings2003} and
\href{http://http://eugeniacheng.com/wp-content/uploads/2017/02/cheng-proofguide.pdf}{``How to
  write proofs''} \citep{Cheng2004}.

\subsection{Motivational Videos}

Feel like you have no motivation to do your homework? Feeling down? These
videos might help!

\begin{itemize}
\item \url{https://www.youtube.com/watch?v=ZXsQAXx_ao0}
\item \url{https://www.youtube.com/watch?v=BQ4yd2W50No}
\item \url{https://www.youtube.com/watch?v=StTqXEQ2l-Y}
\end{itemize}

\end{document}

ext: tex

max_stars_repo_name: jzc/OpenLogic

max_issues_repo_name: jzc/OpenLogic

max_forks_repo_name: jzc/OpenLogic


text:
\subsection{Programming Style and Conventions}\label{sec:ProgrammingStyleConventions}

Everyone has their own preferences and style, and the nature of a project such as COMPAS will reflect that. However, there is a need to suggest some guidelines for programming style, naming conventions etc. Following is a description of some of the elements of programming style and naming conventions used to develop COMPAS v2. These may evolve over time.

\subsubsection{Object-Oriented Programming}\label{sec:OOP}

COMPAS is written in C++, an object-oriented programming (OOP) language, and OOP concepts and conventions should apply throughout the code. There are many texts and web pages devoted to understanding C++ and OOP -- following is a brief description of the key OOP concepts:

\paragraph{Abstraction}\label{sec:Abstraction}\mbox{}

\medskip
For any entity, product, or service, the goal of abstraction is to handle the complexity of the implementation by hiding details that don't need to be known in order to use, or consume, the entity, product, or service. In the OOP paradigm, hiding details in this way enables the consumer to implement more complex logic on top of the provided abstraction without needing to understand the hidden implementation details and complexity. (There is no suggestion that consumers shouldn't understand the implementation details, but they shouldn't need to in order to consume the entity, product, or service).

Abstraction in C++ is achieved via the use of \textit{objects} -- an object is an instance of a \textit{class}, and typically corresponds to a real-world object or entity (in COMPAS, usually a star or binary star). An object maintains the state of an object (via class member variables), and provides all necessary means of changing the state of the object (by exposing public class member functions (methods)). A class may expose public functions to allow consumers to determine the value of class member variables (``getters''), and to set the value of class member variables (``setters'').

\paragraph{Encapsulation}\label{sec:Encapsulation}\mbox{}

\medskip
Encapsulation binds together the data and functions that manipulate the data in an attempt to keep both safe from outside interference and accidental misuse. An encapsulation paradigm that does not allow calling code to access internal object data and permits access through functions only is a strong form of abstraction. C++ allows developers to enforce access restrictions explicitly by defining class member variables and functions as \textit{private}, \textit{protected}, or \textit{public}. These keywords are used throughout COMPAS to enforce encapsulation.

There are very few circumstances in which a consumer should change the value of a class member variable directly (via the use of a setter function) -- almost always consumers should present new situational information to an object (via a public member function), and allow the object to respond to the new information. For example, in COMPAS, there should be almost no reason for a consumer of a star object to directly change (say) the radius of the star -- the consumer should inform the star object of new circumstances or events, and allow the star object to respond to those events (perhaps changing the value of the radius of the star). Changing a single class member variable directly introduces the possibility that related class member variables (e.g. other attributes of stars) will not be changed accordingly. Moreover, developers changing the code in the future should, in almost all cases, expect that the state of an object is maintained consistently by the object, and that there should be no unexpected side-effects caused by calling non class-member functions. In short, changing the state of an object outside the object is potentially unsafe and should be avoided where possible.

\paragraph{Inheritance}\label{sec:Inheritance}\mbox{}

\medskip
Inheritance allows classes to be arranged in a hierarchy that represents \textit{is-a-type-of} relationships. All \textit{non-private} class member variables and functions of the parent (base) class are available to the child (derived) class (and, therefore, child classes of the child class). This allows easy re-use of the same procedures and data definitions, in addition to describing real-world relationships in an intuitive way. C++ allows multiple inheritance -- a class may inherit from multiple parent classes.

Derived classes can define additional class member variables (using the \textit{private}, \textit{protected}, and \textit{public} access restrictions), which will be available to any descendent classes (subject to inheritance rules), but will only be available to ancestor classes via the normal access methods (getters and setters).

\paragraph{Polymorphism}\label{sec:Polymorphism}\mbox{}

\medskip
Polymorphism means \textit{having many forms}. In OOP, polymorphism occurs when there is a hierarchy of classes and they are related by inheritance.

Following the discussion above regarding inheritance, in the OOP paradigm, and C++ specifically, derived classes can override methods defined by ancestor classes, allowing a derived class to implement functions specific to its circumstances. This means that a call to a class member function will cause a different function to be executed depending on the type of object that invokes the function. Descendent classes of a class that has overridden a base class member function inherit the overridden function (but can override it themselves).

COMPAS makes heavy use of inheritance and polymorphism, especially for the implementation of the different stellar types.

\subsubsection{Programming Style}\label{sec:ProgrammingStyle}\mbox{}

\medskip
The goal of coding to a suggested style is readability and maintainability -- if many developers implement code in COMPAS with their own coding style, readability and maintainability will be more difficult than if a consistent style is used throughout the code. Strict adherence isn't really necessary, but it will make it easier on all COMPAS developers if the coding style is consistent throughout.

\subparagraph{Comments}\label{sec:Comments}\mbox{}

\medskip
An old, but good, rule-of-thumb is that any file that contains computer code should be about one-third code, one-third comments, and one-third white space. Adhering to this rule-of-thumb just makes the code a bit easier on the eye, and provides some description (at least of the intention) of the implementation.

\subparagraph{Braces}\label{sec:Braces}\mbox{}

\medskip
The placement of braces in C++ code (actually, any code that uses braces to enclose scope) is a contentious issue, with many developers having long-held, often dogmatic preferences.  COMPAS (so far) uses the K\&R style ("the one true brace style") - the style used in the original Unix kernel and Kernighan and Ritchie's book \textit{The C Programming Language}.

The K\&R style puts the opening brace on the same line as the control statement:

\tabto{2em}while (x == y) \lcb \\
\tabto{4em}something(); \\
\tabto{4em}somethingelse(); \\
\rcb

Note also the space between the keyword \textit{while} and the opening parenthesis, surrounding the == operator, and between the closing parenthesis and the opening brace. Spaces here helps with code readability. Surrounding all arithmetic operators with spaces is preferred.

\subparagraph{Indentation}\label{sec:Indentation}\mbox{}

\medskip
There is ongoing debate in the programming community as to whether indentation should be achieved using spaces or tabs (strange, but true{\dots}).  The use of spaces is more common.  COMPAS (so far) has a mix of both -- whatever is convenient (pragmatism is your friend...). Unfortunately a mix of spaces and tabs doesn't work well with some editors - we should settle on one method and try to stick to it.

COMPAS (mostly) uses an indentation size of 4 spaces - again we should settle on a soze and stick to it.

\subparagraph{Function Parameters}\label{sec:FunctionParameters}\mbox{}

\medskip
In most cases, function parameters should be input only -- meaning that the values of function parameters should not be changed by the function. Anything that needs to be changed and returned to the caller should be returned as a functional return. There are a few exceptions to this in COMPAS -- all were done for performance reasons, and are documented in the code.

To avoid unexpected side-effects, developers should expect (in most cases) that any variables they pass to a function will remain unchanged -- all changes should be returned as a functional return.

\subparagraph{Performance \& Optimisation}\label{sec:PerformanceOptimisation}\mbox{}

\medskip
In general, COMPAS developers should code for performance -- within reason. Bear in mind that many functions will be called many, many thousands of times (in some cases, millions) in one execution of the program.

\begin{itemize}
    \item{Avoid calculating values inside loops that could be calculated once outside the loop.}

    \item{Try to use constants where possible.}

    \item{Use multiplication in preference to functions such as \textit{pow()} and \textit{sqrt()} (note that \textit{pow()} is very expensive computationally; \textit{sqrt()} is expensive, but much less expensive than \textit{pow()}).}

    \item{Don't optimise to the point that readability and maintainability is compromised. Bear in mind that most compilers are good at optimising, and are very forgiving of less-than-optimally-written code (though they are not miracle workers...).}
\end{itemize}

\subsubsection{Naming Conventions}\label{sec:NamingConventions}\mbox{} \\

COMPAS (so far) uses the following naming conventions:

\begin{itemize}
    \item{All variable names should be in camelCase -- don't use underscore\_to\_separate\_words.}

    \item{Function names should be in camelCase, beginning with an uppercase letter.  Function names should be descriptive.}

    \item{Class member variable names are prefixed with "m\_", and the character immediately following the prefix should be uppercase (in most cases -- sometimes, for well-known names or words that are always written in lowercase, lowercase might be used).}

    \item{Local variable names are just camelCase, beginning with a lowercase letter (again, with the caveat that sometimes, for well-known names or words that are always written in uppercase, uppercase might be used).}

    \item{Function parameter names are prefixed with ``p\_'', and the character immediately following the prefix should be uppercase (again, with the caveat that sometimes, for well-known names or words that are always written in lowercase, lowercase might be used).}
\end{itemize}


text:
\chapterimage{chapter_head_2.pdf}

\chapter{Classes}

\section{Classes Should Be Small!}\index{Classes!Classes Should Be Small!}

The fact that classes should be small is not an issue to be discusses, the question, during review, process is \textit{how small should it be?} With functions we measured size by counting physical lines. With classes we use a different measure. We count \textit{responsibilities}.

The name of a class should describe what responsibilities it fulfills. If we cannot derive a concise name for a class, then it's likely too large and has too many responsibilities.

We should also be able to write a brief description of the class in about 25 words,
without using the words "if", "and", "or", or "but".

\begin{marker}
A class's single responsibility should be shown by a self-descriptive name and a short class description.
\end{marker}

\subsection{The Single Responsibility Principle}

\begin{definition}[Single Responsibility Principle (SRP)]
A class or module should have one and only one, reason to change.
\end{definition}

Trying to identify responsibilities (reasons to change) often helps us recognize and
create better abstractions in our code. During code review, you could ask author to break a class into separate ones by methods

\begin{marker}
Review to make sure classes satisfies SRP and ask the author to break classes into multiples if not.
\end{marker}

\subsubsection{Why Do We Often Violate SRP?}

Getting software to work and making software clean are two very different activities. Most of us have limited room in our heads, so we focus on getting our code to work more than organization and cleanliness. This is wholly appropriate. Maintaining a separation of concerns is just as important in our programming activities as it is in our programs

The problem is that too many of us think that we are done once the program works. We fail to switch to the other concern of organization and cleanliness.

\textbf{This is an important reason we do code review - Once we review to make sure code works, we would want to do a second review to make sure code is also clean}

\begin{marker}
There should be at least 2 review process:

\begin{enumerate}
    \item make sure code works
    \item make sure code is clean
\end{enumerate}
\end{marker}

\subsection{Cohesion}

\textbf{Classes should have a small number of instance variables}. Each of the methods of a class should manipulate one or more of those variables. In general the more variables a method manipulates the more cohesive that method is to its class. A class in which each variable is used by each method is maximally cohesive.

In general it is neither advisable nor possible to create such maximally cohesive classes; on the other hand, we would like cohesion to be high. When cohesion is high, it means that the methods and variables of the class are co-dependent and hang together as a logical whole.

Since a class should have small number of instance variables, the strategy of keeping functions small and keeping parameter lists short can sometimes lead to a proliferation of instance variables that are used by a subset of methods. When this happens, it almost always means that there is at least one other class trying to get out of the larger class. You should try to separate the variables and methods into two or more classes such that the new classes are more cohesive.

\begin{marker}
Keep the number of instance variables small and breaks a class apart when there are many
\end{marker}

\subsubsection{Application of Cohesion - Keeping Classes Small}

Consider a large function with many variables declared within it. Let's say you
want to extract one small part of that function into a separate function. However, the code you want to extract uses four of the variables declared in the function. Must you pass all four of those variables into the new function as arguments?

Not at all! If we promoted those four variables to instance variables of the class, then we could extract the code without passing any variables at all. It would be easy to break the function up into small pieces.

Unfortunately, this also means that our classes lose cohesion because they accumulate more and more instance variables that exist solely to allow a few functions to share them. But wait! If there are a few functions that want to share certain variables, doesn't that make them a class in their own right? Of course it does. When classes lose cohesion, split them!

So breaking a large function into many smaller functions often gives us the opportunity to split several smaller classes out as well. This gives our program a much better organization and a more transparent structure.

\begin{marker}
During review process, see if splitting methods could lead to smaller classes.
\end{marker}

\section{Organizing for Change}\index{Classes!Organizing for Change}

max_stars_repo_name: QubitPi/jersey-fundamentals

max_issues_repo_name: QubitPi/jersey-fundamentals

max_forks_repo_name: QubitPi/jersey-fundamentals


text:
\documentclass{llncs}

\usepackage[T1]{fontenc}
%\usepackage[utf8x]{inputenc}
\usepackage[english]{babel}
\usepackage{lmodern}
\usepackage{mathtools}
%\usepackage{fullpage}
\usepackage{graphicx}
\usepackage{xspace}
\usepackage{tabularx}
\usepackage[lf]{ebgaramond}
\usepackage{biolinum}
\usepackage[cmintegrals,cmbraces]{newtxmath}
\usepackage{ebgaramond-maths}
\usepackage{sectsty}
\usepackage[noend]{algpseudocode}
\usepackage{algorithm}
\usepackage{algorithmicx}
\usepackage{multicol}
\usepackage[bottom]{footmisc}
\usepackage{caption}
%\captionsetup{font=small}
\captionsetup{labelfont={sf,bf}}
\usepackage{CJKutf8}
\def\DevnagVersion{2.17}
\usepackage{devanagari}

\PassOptionsToPackage{hyphens}{url}\usepackage{hyperref}

\makeatletter
  % Recover some math symbols that were masked by eb-garamond, but do not
  % have replacement definitions.
  \DeclareSymbolFont{ntxletters}{OML}{ntxmi}{m}{it}
  \SetSymbolFont{ntxletters}{bold}{OML}{ntxmi}{b}{it}
  \re@DeclareMathSymbol{\leftharpoonup}{\mathrel}{ntxletters}{"28}
  \re@DeclareMathSymbol{\leftharpoondown}{\mathrel}{ntxletters}{"29}
  \re@DeclareMathSymbol{\rightharpoonup}{\mathrel}{ntxletters}{"2A}
  \re@DeclareMathSymbol{\rightharpoondown}{\mathrel}{ntxletters}{"2B}
  \re@DeclareMathSymbol{\triangleleft}{\mathbin}{ntxletters}{"2F}
  \re@DeclareMathSymbol{\triangleright}{\mathbin}{ntxletters}{"2E}
  \re@DeclareMathSymbol{\partial}{\mathord}{ntxletters}{"40}
  \re@DeclareMathSymbol{\flat}{\mathord}{ntxletters}{"5B}
  \re@DeclareMathSymbol{\natural}{\mathord}{ntxletters}{"5C}
  \re@DeclareMathSymbol{\star}{\mathbin}{ntxletters}{"3F}
  \re@DeclareMathSymbol{\smile}{\mathrel}{ntxletters}{"5E}
  \re@DeclareMathSymbol{\frown}{\mathrel}{ntxletters}{"5F}
  \re@DeclareMathSymbol{\sharp}{\mathord}{ntxletters}{"5D}
  \re@DeclareMathAccent{\vec}{\mathord}{ntxletters}{"7E}

  % Change font for algorithm label.
  \renewcommand\ALG@name{\sffamily\bfseries Algorithm}
\makeatother

\allsectionsfont{\sffamily}
\pagestyle{plain}

\makeatletter
\renewenvironment{abstract}{%
      \list{}{\advance\topsep by0.35cm\relax\small
      \leftmargin=1cm
      \labelwidth=\z@
      \listparindent=\z@
      \itemindent\listparindent
      \rightmargin\leftmargin}\item[\hskip\labelsep
                                    \textsf{\textbf{\abstractname}}]}
    {\endlist}
\newenvironment{extranote}{%
      \list{}{\advance\topsep by0.35cm\relax\small
      \leftmargin=1cm
      \labelwidth=\z@
      \listparindent=\z@
      \itemindent\listparindent
      \rightmargin\leftmargin}\item[\hskip\labelsep
                                    \textsf{\textbf{Note.}}]}
    {\endlist}
\makeatother

\spnewtheorem{mtheorem}{Theorem}{\sffamily\bfseries}{\itshape}
\spnewtheorem*{mproof}{Proof}{\sffamily\bfseries\itshape}{\rmfamily}

%\newcommand{\GF}{\mathrm{\textit{GF}}}
\newcommand{\GF}{GF}
\newcommand{\QR}{QR}
\newcommand{\bN}{\mathbb{N}}
\newcommand{\bZ}{\mathbb{Z}}
\newcommand{\bR}{\mathbb{R}}
\newcommand{\bF}{\mathbb{F}}
\newcommand{\bG}{\mathbb{G}}
\newcommand{\neutral}{\mathbb{O}}
\newcommand{\vol}{\text{vol}}
\newcommand{\bitlength}{\text{len}}
\newcommand{\smod}[1]{\,\,\,(\text{mod}^{\pm} #1)}
\newcommand{\smodnospace}[1]{(\text{mod}^{\pm} #1)}

\newcommand{\indentA}{\hbox to 1.5em{\hfill}}
\newcommand{\indentB}{\hbox to 3em{\hfill}}
\newcommand{\indentC}{\hbox to 4.5em{\hfill}}
\newcommand{\indentD}{\hbox to 6em{\hfill}}

\raggedbottom

\begin{document}

\title{\textsf{Optimized Binary GCD for Modular Inversion}}

\author{Thomas Pornin}
\institute{NCC Group, \email{thomas.pornin@nccgroup.com}}

\maketitle
\noindent\makebox[\textwidth]{23 August 2020} % FIXME

\begin{abstract}
In this short note, we describe a practical optimization of the
well-known extended binary GCD algorithm, for the purpose of computing
modular inverses. The method is conceptually simple and is applicable
to all odd moduli (including non-prime moduli). When implemented
for inversion in the field of integers modulo the prime $2^{255}-19$, on
a recent x86 CPU (Coffee Lake core), we compute the inverse in 6253
cycles, with a fully constant-time implementation.
\end{abstract}

\begin{extranote}
An initial version of this note used a slightly more aggressive
algorithm ($k$ iterations of the inner loop instead of $k-1$) which was
wrong. The gap in the proof was discovered by Tim Taubert, and further
discussions showed a counterexample for which the algorithm failed.
The current version of the algorithm, and the revised proof, are believed
correct. A side effect of the fix is to enable extra optimizations in
the implementation of the inner loop, leading to substantially better
performance, now at 6253 cycles (from an original 7490 cycles).
\end{extranote}

\section{Introduction}\label{sec:intro}

We study here the following problem: given an odd integer $m \geq 3$ and
a value $y \in \bZ_m$, compute its modular inverse $x$ which is such
that $xy = 1 \bmod m$. This is a notoriously expensive operation,
compared with addition and multiplication in $\bZ_m$. Its cost is the
main reason why, for instance, operations on elliptic curves are often
done with projective coordinates or other similar systems where
coordinates are handled as fractions, so that all costly inversions can
be mutualized into a final inversion, once the whole computation is
finished.

There are some applications of modular inversion where the value to
invert and/or the modulus is secret. In elliptic curve cryptography, the
modulus $m$ is normally prime and publicly known, but the value to
invert is private. In RSA key pair generation, modular inversion is
needed to compute the private exponent, and the CRT reconstruction
factor; in these cases, the value to invert and the modulus are both
secret, and when computing the private exponent from the public
exponent, the modulus is not even prime. In all situations where secret
values are handled, a \emph{constant-time} implementation is needed
(i.e. one with an execution time and a memory access pattern that
do not vary depending on secret data).

When $m$ is prime, a simple method is to use Fermat's little theorem:
\begin{equation*}
    1/y = y^{m-2} \mod m
\end{equation*}
This method requires $O(\log m)$ multiplications modulo $m$, hence it
usually has a cost cubic in the size of $m$ (assuming a quadratic cost
for multiplication, which is the usual case for moduli used in
cryptographic applications). When $m = 2^{255}-19$, and implementing on
a recent 64-bit x86 CPU, modular multiplications are quite fast, and
this is favourable to Fermat's little theorem; our implementation of
this inversion method, on an Intel Core i5-8259U at 2.30~GHz (Coffee
Lake core), can perform a constant-time inversion in 9175 cycles. The
goal is to find something faster.

Apart from Fermat's little theorem, the main method to compute modular
inverses is the Extended Euclidean Algorithm for GCD. While Euclid's
algorithm was presented by Euclid himself\cite{Euc300BC}, the extension
to compute values $u$ and $v$ such that $au + bv = \text{GCD}(a,b)$ was
first described by \={A}ryabha\d{t}a\cite{Ary499}. The binary GCD is a
variant of Euclid's algorithm that performs only comparisons,
subtractions and divisions by 2 (i.e. right shifts), and is therefore
more amenable to fast and constant-time implementations than the general
divisions in Euclid's algorithm; that algorithm was first presented,
albeit a bit cryptically, in the Nine Chapters on the Mathematical
Art\cite{Chi1}. A more accessible description was published by Stein in
1967\cite{Ste1967}. A variant called the plus-minus algorithm was
presented by Brent and Kung in 1983\cite{BreKun1983}. More recently, a
new (and quite complex) algorithm was described by Bernstein and
Yang\cite{BerYan2019}. See also \cite{NatSar2018} for extensive
benchmarks for inversions modulo primes with the $2^n - r$ format (for a
small $r$). Asymptotically fast but not constant-time GCD algorithms are
presented in~\cite{Mol2008}.

Bernstein and Yang claim to hold the current speed record for computing
inverses modulo $2^{255}-19$ on recent 64-bit x86 CPUs (Skylake cores
and similar variants). The field of integers modulo the prime $p =
2^{255}-19$ is used by some well-known elliptic curves, Curve25519 and
Edwards25519, thus highlighting the relevance of this particular
modulus. Bernstein and Yang obtain a performance of 8520 cycles on a
Coffee Lake core\footnote{I am deeply grateful to Bo-Yin Yang, who sent
me a copy of their code so that I may benchmark it on the same test
system as the one used for the algorithm described in this note. This
speed is almost identical to the 8543 cycles on a Kaby Lake core, as
reported in the published article; this makes sense, since Coffee Lake
and Kaby Lake use the same internal scheduler and execution units, and
differ only in external pieces that should not matter much for a
compact, L1-cache only computational piece of code.}. With the method
presented here, we do better: 6253 cycles.

Our algorithm was added to the BearSSL library\cite{BearSSL} in 2018; it
was also deployed as part of the implementation of key pair generation
in the reference code for the Falcon post-quantum signature
scheme\cite{Falcon}. This note formally describes the algorithm and
proves its correctness.

It shall be said that gaining even 2000~cycles on inversion modulo
$2^{255}-19$ on recent x86 CPUs will not yield substantial performance
gains in elliptic curve implementations that use that field: a curve
point multiplication will typically use close to 100000~cycles or more,
and a 2000-cycle gain will be hardly noticeable, even in the improbable
situation of a system where raw elliptic curve performance is the
bottleneck. In practice, on large server and desktop systems,
cryptographic operations, especially asymmetric cryptographic
operations, use only a relatively small proportion of the available
power, and a 2\% gain on these operations will not really matter.
However, on smaller architectures with less efficient multipliers, or
when using larger moduli, the performance ratio with Fermat's little
theorem is expected to increase, making the optimization worthwhile. In
particular, on small microcontrollers (e.g. ARM Cortex M0+), the
inversion algorithm presented in this note may be fast enough that a
curve point multiplication routine using window optimizations would
benefit from using it in order to normalize the window points to affine
coordinates, for faster subsequent point additions.

\section{Optimizing the Extended Binary GCD}

Algorithm~\ref{alg:bingcd:classic} describes the classic extended
binary GCD.

\begin{algorithm}[H]
    \caption{\ \ Extended Binary GCD (classic algorithm)\label{alg:bingcd:classic}}
    \begin{algorithmic}[1]
        \Require{Odd modulus $m$ ($m \geq 3$, $m \bmod 2 = 1$) and
        value to invert $y$ ($0 \leq y < m$)}
        \Ensure{$1/y \bmod m$ (if $\text{GCD}(y,m) = 1$), or zero}
        \State{$a \leftarrow y$, $u\leftarrow 1$,
            $b\leftarrow m$, $v\leftarrow 0$}
        \While{$a \neq 0$}
            \If{$a = 0 \bmod 2$}
                \State{$a \leftarrow a/2$}
                    \Comment{$a$ is even, so this division is exact.}
                \State{$u \leftarrow u/2 \bmod m$}
            \Else
                \If{$a < b$}
                    \State{$(a, u, b, v) \leftarrow (b, v, a, u)$}
                        \Comment{Conditional swap to ensure $a\geq b$.}
                \EndIf
                \State{$a \leftarrow (a - b)/2$}
                    \Comment{$a$ and $b$ are odd, so this division is exact.}
                \State{$u \leftarrow (u - v)/2 \bmod m$}
            \EndIf
        \EndWhile
        \If{$b \neq 1$}
            \State{\Return 0 (value $y$ is not invertible)}
                \Comment{$b$ contains $\text{GCD}(y, m)$ at this point.}
        \EndIf
        \State{\Return $v$}
    \end{algorithmic}
\end{algorithm}

In this algorithm, values $a$ and $b$ are nonnegative integers, while
$u$ and $v$ are integers modulo $m$. It is not hard to see that:
\begin{itemize}
    \item $b$ is always an odd integer;
    \item $\text{GCD}(a, b)$ is constant throughout the algorithm;
    \item when $a = 0$ (exit condition), $b$ contains the GCD of $y$
    and $m$;
    \item divisions by 2 to compute the new value of $a$ are exact;
    \item each iteration maintains the invariants $a = uy \bmod m$ and
    $b = vy \bmod m$;
    \item the sum of the lengths, in bits, of $a$ and $b$ is decremented
    by at least one at each iteration.
\end{itemize}
Therefore, if $m$ is a $n$-bit integer, then at most $2n-1$ iterations
are needed to reach the exit condition. In order to have a constant-time
implementation, we need to \emph{always} perform $2n-1$ iterations; if
$a = 0$, the extra iterations will not modify $b$ or $v$, so they won't
modify the result\footnote{As will be explained later on, the last
iteration needs not be done if the input $y$ is known to be
invertible.}. A constant-time implementation will also need to perform
each iteration in a constant-time way, replacing all conditions with
bitwise combinations of values.

The extended binary GCD is quadratic in the length of the modulus: it
has $O(n)$ iterations, and each iteration involves a fixed number of
operations which can be computed in time $O(n)$ (comparisons,
subtractions, halvings). If implemented as is, its performance is
comparable to Fermat's little theorem for $m = 2^{255}-19$, but somewhat
worse for platforms with fast large multiplications; for larger moduli,
it becomes better. It also properly handles non-prime moduli.

We will now present two complementary optimization methods that apply to
this algorithm and substantially improve performance in practice.

\paragraph{Optimization 1: mutualize updates to \emph{u} and \emph{v}:}

In the description above, $u$ and $v$ record the modifications applied
to $a$ and $b$, but no decision is taken based on their values. We can
thus delay and group the updates to $u$ and $v$:
\begin{itemize}

    \item Initialize the \emph{update factors}:
    $f_0 \leftarrow 1$, $g_0 \leftarrow 0$,
    $f_1 \leftarrow 0$, $g_1 \leftarrow 1$ \\
    The update factors are such that the new value of $(u, v)$ should be:
    \begin{equation*}
        (u, v) \leftarrow (f_0 u + g_0 v, f_1 u + g_1 v)
    \end{equation*}
    We update $f_0$, $f_1$, $g_0$ and $g_1$ instead of $u$ and $v$ in
    the algorithm. Once every $r$ iterations, while the update factors
    are still small (e.g. they fit in a single register each), we apply
    them to $u$ and $v$ using multiplication opcodes, which is much
    faster than making $r$ linear updates. After each such update, the
    update factors are reinitialized.

    \item When swapping $a$ with $b$ (in case $a$ and $b$ are both odd,
    and $a < b$), exchange $f_0$ with $f_1$ and $g_0$ with $g_1$.

    \item When subtracting $b$ from $a$, subtract $f_1$ from $f_0$ and
    $g_1$ from $g_0$.

    \item At each iteration, we should divide $f_0$ and $g_0$ by 2
    (corresponding to the division of $a$ by 2), but this would entail
    making them non-integral. Instead, we multiply $f_1$ and $g_1$ by 2.

\end{itemize}

The effect of the last item is that the update factors, after $r$
iterations, are in fact themselves wrong by a factor of exactly $2^r$.
This is not much of a problem; we can perform a division of the final
result $v$ by $2^{2n-1}$, after the $2n-1$ iterations (if $m$ is known
in advance, this division is merely a multiplication by a precomputed
constant; even in the general case, this division is the same thing as a
pair of Montgomery reductions, and has a cost no higher than that of a
couple of modular multiplications).

The update factors are signed integers, and we can bound them: after
$t \geq 1$ iterations, we have:
\begin{equation*}
    \begin{gathered}
    -2^t < f_0, g_0 < +2^t \\
    -2^t < f_1, g_1 \leq +2^t
    \end{gathered}
\end{equation*}
This is easily seen to be true after one iteration: the only possible
values for $(f_0, g_0)$ are then $(1, 0)$ ($a$ was even), $(1, -1)$ ($a$
was odd, and $b$ was lower than $a$) or $(-1, 1)$ ($a$ was odd and lower
than $b$, and a swap occurred before the subtraction); $(f_1, g_1)$ will
be either $(0, 2)$ or $(2, 0)$. Similarly, if the relations are verified
after $t$ iterations, then it can be seen that all three possible
outcomes of the next iteration will lead to the relations still being
verified after $t+1$ iterations.

A consequence is that it is possible to encode update factors as
\emph{unsigned} integers of $r+1$ bits, by storing $f_0 + 2^r - 1$,
$f_1 + 2^r - 1$... In particular, if the current architecture has
registers of $2r+2$ bits, then two update factors may be stored in the
same register. This will be used in our most optimized implementation.

\paragraph{Optimization 2: use approximations of \emph{a} and \emph{b}:}

While the update factors allow grouping the updates of $u$ and $v$, and
performing them efficiently with the help of integer multiplication
opcodes, the values of $a$ and $b$ are still large. However, the update
factors \emph{record} the operations which have been performed on $a$
and $b$, to be mimicked by $u$ and $v$; thus, if we could compute the
update factors \emph{first}, we may also apply them to $a$ and $b$,
thereby grouping and optimizing the updates to $a$ and $b$ too.

The issue here is that all update decisions depend on the values of $a$
and $b$; specifically, whether $a$ is even or odd (thus using the low
bit of $a$), and whether $a$ is greater than $b$ (which requires mostly
looking at the high bits of $a$ and $b$). It can easily be seen that for
$k$ iterations, the value of the low bit of $a$ will only depend on the
initial $k$ low bits of $a$ and $b$. The optimization idea is to make
\emph{approximations} of $a$ and $b$ that will fit in a single register
each, and can tell us in a mostly correct way whether $a$ is greater
than $b$ or not.

Suppose that the local architecture offers $2k$-bit registers. Let
$n = \text{max}(\bitlength(a), \bitlength(b))$, where $\bitlength(x)$ is
the length of $x$ in bits (the smallest integer $i$ such that
$x < 2^i$). Thus, $n$ is the current length of the largest of $a$ and $b$.
If $n \leq 2k$, then the values $a$ and $b$ already fit into a single
register each, and we can use the classic extended GCD algorithm to
compute the exact update factors. Otherwise, extract the $k-1$ low bits
and the $k+1$ top bits of $a$ and $b$, and assemble them into $\bar a$ and
$\bar b$:
\begin{align*}
    \bar a &= (a \bmod 2^{k-1}) + 2^{k-1} \lfloor a / 2^{n-k-1} \rfloor \\
    \bar b &= (b \bmod 2^{k-1}) + 2^{k-1} \lfloor b / 2^{n-k-1} \rfloor
\end{align*}
Then, apply the extended binary GCD on $\bar a$ and $\bar b$ for $k-1$
iterations, yielding update factors $f_0$, $g_0$, $f_1$ and $g_1$.
Finally, apply these factors to update $a$, $b$, $u$ and $v$:
\begin{align*}
    (a, b) &\leftarrow ((f_0 a + g_0 b) / 2^{k-1}, (f_1 a + g_1 b) / 2^{k-1}) \\
    (u, v) &\leftarrow (f_0 u + g_0 v, f_1 u + g_1 v)
\end{align*}

The divisions by $2^{k-1}$ are exact: while $\bar a$ and $\bar b$ are
approximations of $a$ and $b$, their low $k-1$ bits match exactly, and
the application of the update factors clears the $k-1$ low bits of $a$
and $b$.

The computed update factors are still approximations; they may be
slightly wrong, i.e. fail to shrink $a$ and $b$ as much as expected.
They may also lead to a negative value for $a$ or $b$; in case a
negative $a$ is obtained, it suffices to negate it, and also to negate
$u$ (this time modulo $m$), and similarly for $b$. In practice, it is
convenient to:
\begin{enumerate}
    \item update $a$ with the update factors $f_0$ and $g_0$;
    \item if the new value is negative, negate it, and also negate
    $f_0$ and $g_0$;
    \item apply the (possibly negated) update factors $f_0$ and $g_0$
    to compute the new value of $u$.
\end{enumerate}

The tricky point is computing how many iterations are needed in the
worst case; for a constant-time implementation, we need an absolute
bound. Appendix~\ref{sec:iterbound} is dedicated to proving that each
sequence of $k-1$ iterations ensures that $\bitlength(a) + \bitlength(b)$
is reduced by at least $k-1$; therefore, $2\bitlength(m)-1$ iterations are
sufficient (just like the classic binary GCD algorithm). This bound
is optimal: if $y = 2^{\bitlength(m)-1}$, then all $2\bitlength(m)-1$
iterations are needed.

It may be convenient, for implementation purposes, to slightly raise the
total number of iterations so that it is a multiple of $k-1$.
Alternatively, the last iterations may be specialized away, which may
allow for some extra optimizations, especially when the modulus $m$ is
prime; this will be detailed later on.

The optimized algorithm is described below
(algorithm~\ref{alg:bingcd:optimized}).

\begin{algorithm}[H]
    \caption{\ \ Extended Binary GCD (optimized algorithm)\label{alg:bingcd:optimized}}
    \begin{algorithmic}[1]
        \Require{Odd modulus $m$ ($m \geq 3$, $m \bmod 2 = 1$),
        value $y$ ($0 \leq y < m$), and $k > 1$}
        \Ensure{$1/y \bmod m$ (if $\text{GCD}(y,m) = 1$)}
        \State{$a \leftarrow y$, $u\leftarrow 1$,
            $b\leftarrow m$, $v\leftarrow 0$}
        \For{$1 \leq i \leq \lceil (2\bitlength(m)-1)/k \rceil$}
            \State{$n \leftarrow \text{max}(\bitlength(a), \bitlength(b), 2k)$}
            \State{$\bar a \leftarrow (a \bmod 2^{k-1}) + 2^{k-1} \lfloor a / 2^{n-k-1} \rfloor$}
                \Comment{$\bar a < 2^{2k}$}
            \State{$\bar b \leftarrow (b \bmod 2^{k-1}) + 2^{k-1} \lfloor b / 2^{n-k-1} \rfloor$}
                \Comment{$\bar b < 2^{2k}$}
            \State{\label{alg:bingcd:optimized:inner}$f_0 \leftarrow 1$, $g_0 \leftarrow 0$, $f_1 \leftarrow 0$, $g_1 \leftarrow 1$}
            \For{$1 \leq j \leq k-1$}
                \If{$\bar a = 0 \bmod 2$}
                    \State{$\bar a \leftarrow \bar a/2$}
                \Else
                    \If{$\bar a < \bar b$}
                        \State{$(\bar a, \bar b) \leftarrow (\bar b, \bar a)$}
                        \State{$(f_0, g_0, f_1, g_1) \leftarrow (f_1, g_1, f_0, g_0)$}
                    \EndIf
                    \State{$\bar a \leftarrow (\bar a - \bar b)/2$}
                    \State{$(f_0, g_0) \leftarrow (f_0 - f_1, g_0 - g_1)$}
                \EndIf
                \State{$(f_1, g_1) \leftarrow (2f_1, 2g_1)$}
            \EndFor
            \State{$(a, b) \leftarrow ((af_0 + bg_0) / 2^{k-1}, (af_1 + bg_1) / 2^{k-1})$}
            \If{$a < 0$}
                \State{$(a, f_0, g_0) \leftarrow (-a, -f_0, -g_0)$}
            \EndIf
            \If{$b < 0$}
                \State{$(b, f_1, g_1) \leftarrow (-b, -f_1, -g_1)$}
            \EndIf
            \State{$(u, v) \leftarrow (uf_0 + vg_0 \bmod m, uf_1 + vg_1 \bmod m)$}
        \EndFor
        \State{$v \leftarrow v/2^{(k-1)\lceil (2\bitlength(m)-1)/(k-1) \rceil} \bmod m$}
        \If{$b \neq 1$}
            \State{\Return 0 (value $y$ is not invertible)}
        \EndIf
        \State{\Return $v$}
    \end{algorithmic}
\end{algorithm}

When $y$ is known in advance to be either 0 or an invertible number
modulo $m$ (this is always the case when $m$ is prime), then the last
two outer iterations can be specialized and simplified, because, at
that point:
\begin{itemize}

    \item If $y = 0$, then the update factors after $k$ inner
    iterations will be $f_0 = 1$, $g_0 = 0$, $f_1 = 0$ and $g_1 = 2^{k-1}$,
    regardless of the value of $\bar b$.

    \item Otherwise, it is known that $\bitlength(a) + \bitlength(b) <
    2k$, which means that $\bar a = a$ and $\bar b = b$; the computed
    update factors will then be exact, and will remain exact even if
    more than $k-1$ iterations are performed. Depending on the
    implementation of the routines that apply the update factors, it may
    be possible to merge the last two outer iterations.

    \item At the end of the final iteration, it is not necessary to
    compute the new values of $a$, $b$ and $u$, only $v$ (since the
    update factors are exact at that point, they cannot lead to a
    negative $b$).

    \item The last inner loop iteration can be omitted because when the
    GCD is $1$, the last iteration that starts with a non-zero $a$ must
    have $a = 1$ and $b = 1$, and $v$ already contains the correct
    result. This reduces the total number of required iterations to
    $2\bitlength(m) - 2$, i.e. 508 for a 255-bit modulus.

\end{itemize}

In the general case, when the modulus $m$ is not necessarily prime, and
a non-invertible but distinct from zero input $y$ is possible, these
optimizations on the final steps can still be applied, but an extra
verification step is required. It suffices to verify that the computed
value, when multiplied with the original operand, yields the value 1.
The cost of this extra check is small (typically less than 2\% of the
overall cost, for a modulus of about 256 bits).

\section{Implementation and Benchmarks}

We implemented our method in the case of the prime modulus
$p = 2^{255}-19$. Our code is available at:
\begin{center}
    \url{https://github.com/pornin/bingcd}
\end{center}
Implementation is in C, with some inline assembly, as well as calls to
intrinsic functions. It should compile and run on any Linux or similar
system, using GCC or Clang, and an x86 CPU with BMI2 and ADX extensions
(in the Intel line, this means using a Haswell core or newer).

Extraction of the approximated words $\bar a$ and $\bar b$ is done with
constant-time selection primitives to obtain the highest non-zero limbs
in $a \text{\textsc{\ or\ }} b$ (bitwise Boolean \textsc{or}), then
the \verb+lzcnt+ opcode to count the number of leading zeros. This
opcode is available on all x86 CPUs that also feature the BMI2 opcodes.
This extraction is constant-time under the assumption that bitwise shifts
execute in time independent of the shift count, an assumption that is
valid on all relevant x86 CPUs to date\footnote{The last Intel CPU on
which shifts were not constant-time with regard to the shift count was
the Pentium IV, in the early 2000s, but that CPU does not offer
\texttt{lzcnt}.}.

Updates to $a$, $b$, $u$ and $v$ use $64\times 64\rightarrow 128$
multiplication opcodes, accessed through the non-standard
\verb+unsigned __int128+ type. Large integers are represented over 256
bits as four 64-bit limbs; carry propagation uses the \verb+adc+ and
\verb+sbb+ opcodes, accessed through the \verb+_addcarry_u64()+ and
\verb+_subborrow_u64()+ intrinsic functions.

The whole point of the optimization described in this note is to make
the core of the inner loop as simple and fast as possible. In our
implementation, we use two versions of that core, implemented using
inline assembly. The first version uses six 64-bit registers, to
$\bar a$, $\bar b$, and the four update factors:
\begin{verbatim}
        # Copy f0, g0, f1, g1, xa and xb into r10..r15
        movq   %rax, %r10
        movq   %rbx, %r11
        movq   %rcx, %r12
        movq   %rdx, %r13
        movq   %rsi, %r14
        movq   %rdi, %r15

        # Conditional swap if xa < xb
        cmpq   %rdi, %rsi
        cmovb  %r15, %rsi
        cmovb  %r14, %rdi
        cmovb  %r12, %rax
        cmovb  %r10, %rcx
        cmovb  %r13, %rbx
        cmovb  %r11, %rdx

        # Subtract xb from xa
        subq   %rdi, %rsi
        subq   %rcx, %rax
        subq   %rdx, %rbx

        # If xa was even, override the operations above
        testl  $1, %r14d
        cmovz  %r10, %rax
        cmovz  %r11, %rbx
        cmovz  %r12, %rcx
        cmovz  %r13, %rdx
        cmovz  %r14, %rsi
        cmovz  %r15, %rdi

        # Now xa is even; apply shift
        shrq   $1, %rsi
        addq   %rcx, %rcx
        addq   %rdx, %rdx
\end{verbatim}
Only simple instructions that can be executed in several distinct CPU
ports, thus amenable to parallelism, are used. Conditional swaps and
selection are performed with \verb+cmov+; this instruction is
constant-time, has a 1-cycle latency, and two of them can be executed in
a single cycle on a Coffee Lake core, thus each sequence of six
\verb+cmov+ conceptually needs only 3 cycles to execute. The
non-conditional \verb+mov+ are even more amenable to parallelism (four
ports are eligible, leading to up to four \verb+mov+ executing in the
same cycle) and some may be in fact eliminated by the register renaming
unit. The \verb+sub+ and \verb+add+ opcodes can similarly execute on
four ports. Note that multiplications by 2 of $f_1$ and $g_1$
(\verb+rcx+ and \verb+rdx+ in the code above) are performed with
\verb+add+ instead of \verb+shl+, because the latter may execute on only
two ports and was measured to slightly decrease performance.

Counting each \verb+cmov+ at 0.5~cycle, \verb+sub+, \verb+add+,
\verb+cmp+ and \verb+test+ at 0.25~cycle each, and \verb+shr+ at
0.5~cycle (using the reported reciprocal throughputs from Agner Fog's
comprehensive tables\cite{FogTables}), the sequence above requires a
minimum of 8.25~cycles, not counting the \verb+mov+ opcodes. Simple
measurements show that this sequence, when used in a loop, requires 9
cycles per iteration (on our test system), which means that some of the
\verb+mov+ opcodes are eliminated and replaced with internal free data
routing.

The second version of the core loop is even faster:
\begin{verbatim}
        # Copy f0, g0, f1, g1, xa and xb into r10, r12, r14, r15
        movq   %rax, %r10
        movq   %rcx, %r12
        movq   %rsi, %r14
        movq   %rdi, %r15

        # Conditional swap if xa < xb
        cmpq   %rdi, %rsi
        cmovb  %r15, %rsi
        cmovb  %r14, %rdi
        cmovb  %r12, %rax
        cmovb  %r10, %rcx

        # Subtract xb from xa
        subq   %rdi, %rsi
        subq   %rcx, %rax
        addq   %rdx, %rax

        # If xa was even, override the operations above
        testl  $1, %r14d
        cmovz  %r10, %rax
        cmovz  %r12, %rcx
        cmovz  %r14, %rsi
        cmovz  %r15, %rdi

        # Now xa is even; apply shift
        shrq   $1, %rsi
        addq   %rcx, %rcx
        subq   %rdx, %rcx
\end{verbatim}
This alternate implementation keeps two update factors in the same
register. Specifically, register \verb+rax+ contains $f_0 + 2^{31} - 1 +
2^{32} (g_0 + 2^{31} - 1)$; \verb+rcx+ similarly contains $f_1$ and
$g_1$. The constant $(2^{31} - 1)(2^{32} + 1)$ is kept in \verb+rdx+.
Grouping values by pairs reduces the number of instructions, especially
for conditional swaps and moves. The addition of the constant $2^{31}-1$
to each update factor ensures that the stored values remain positive;
thus, there will be no unwanted carry propagating from the low to high
halves of the registers. The alternate loop appears to use, on average,
about 6.17~cycles per iteration.

In total, the 508 iterations represent about 3260~cycles by themselves;
all the other operations (extraction of the properly shifted words,
updates to $a$, $b$, $u$ and $v$, final multiplication by $2^{-508}$)
use slightly less than 3000~cycles. This highlights the importance of
optimizing the inner loop as much as possible.

Measures are made on an Intel Core i5-8295U CPU, clocked at 2.3~GHz.
Host system is Linux (Ubuntu 20.04, kernel 5.4.0-42). C compiler is
Clang-10.0.0 and compilation uses optimization flags
\verb+-O3 -march=native+. TurboBoost is disabled (in the BIOS), so there
is no frequency scaling above 2.3~GHz. SMT (HyperThreading) is disabled.
Clock cycles are obtained with the \verb+rdtsc+ opcode before and after
a sequence of 100 dependent inversions (i.e. each inversion works over
the output of the previous one); a serializing \verb+lfence+ opcode is
also used to ensure that the CPU scheduler will run \verb+rdtsc+ at the
right place in the sequence of instructions. The bench program first
runs some extensive unit tests, to verify correction of the
computations, and also act as a ``warm-up'' to make sure the CPU has
gone out of energy-saving mode and has reached its nominal frequency.
2000 measures are then performed; the first 1000 are discarded (they are
used to make sure caches are populated and branch prediction is properly
trained), then the median of the remaining 1000 is reported\footnote{The
deactivation of TurboBoost and the warm-up steps ensure that the CPU
runs at exactly the frequency used by the TSC counter. We also tried
using the performance counters, namely the CPU\_CLK.UNHALTED counter
which is supposed to count ``true'' cycles, and can be read with
\texttt{rdpmc} at index \texttt{0x40000001}; this does not yield
different results. Our test program uses \texttt{rdtsc} because it is
readable by non-root userland processes by default. Regardless, some
measurement variance remains, of unclear provenance. It does not seem to
be correlated with the value being inverted (since we observe that
variance even when using the same values again and again), and thus
should not contradict our claim of constant-time behaviour. It is known
that modern complex CPUs are full of accumulated heuristics, many of
them undocumented.}.

On this system, our inversion routine runs in 6253~cycles. This is
better than the previously published record\cite{BerYan2019}, which runs
in 8520~cycles on the same system. For comparison purposes, we also
implemented the usual inversion method with Fermat's little theorem
(FLT); for this, we use handwritten assembly routines for
multiplications, squarings, and sequences of squarings modulo
$2^{255}-19$. These routines are inspired from, and very similar to, the
ones described in~\cite{OliLopHisFazRod2017} and \cite{NatSar2018},
though our FLT inversion routine ends up running in 9175~cycles, which
is marginally faster than the previous record (9301~cycles,
in~\cite{NatSar2018}): apparently, our multiple squaring routine is
faster by about half a cycle per iteration. We don't really know why. In
any case, the extended binary GCD is still faster.

A noteworthy feature of our code is that it can work with other moduli,
with the same performance; supported moduli have the format $2^{255}-r$
for any odd $r$ between 1 and $2^{31}-1$ (this limit is imposed by the
way we perform modular reduction after a multiplication). If the modulus
is not prime, the FLT inversion method no longer works, but the extended
binary GCD still returns the correct result (in that case, our code
includes the extra verification step, to also cover non-invertible inputs;
that extra verification step costs about 100~cycles).

\section{Conclusion}

We described a simple and practical optimization to the classic extended
binary GCD algorithm. In the specific case of inversion modulo
$2^{255}-r$ (for a small $r$), this optimization happens to lead to a
routine which is faster than the usual method based on Fermat's little
theorem, and also slightly faster than the previously best reported
inversion routine for that class of moduli, the recent (and more
complicated) algorithm from Bernstein and Yang. This does not mean that
our optimization makes binary GCD the best algorithm in all cases; it
depends on many parameters, such as modulus size and format and target
implementation architecture. We also did not try to apply it to
polynomial GCDs, though a similar optimization might also provide
practical benefits there. One can view our method as an improvement in
\emph{data locality}, thereby reducing data routing cost (in our primary
example, by keeping values in single registers, thereby avoiding RAM
traffic in the innermost loop).

The same optimization technique should, conceptually, apply to the
computation of the Jacobi symbol\cite{Coh1993}. The Jacobi symbol
is useful to test the quadratic residue status of a finite field element,
which in turn is part of the operations used for constant-time hashing
into elliptic curves\cite{DraftHashToCurve09}.

\vspace{2ex}
\noindent\textsf{\textbf{Acknowledgements:}} We thank Eric Schorn
and Paul Bottinelli, who reviewed this article, and Tim Taubert,
who found the gap in the proof in a previous version.

\begin{thebibliography}{20}

\bibitem{Ary499}
%\={A}ryabha\d{t}a,
{\dn aAy\0BV},
\emph{\={A}ryabha\d{t}\={\i}ya}, 499.

\bibitem{BearSSL}
\emph{BearSSL, a smaller SSL/TLS library},\\
\url{https://www.bearssl.org/}

\bibitem{BerYan2019}
D.~Bernstein and B.-Y. Yang,
\emph{Fast constant-time gcd computation and modular inversion},\\
\url{https://gcd.cr.yp.to/papers.html#safegcd}

\bibitem{BreKun1983}
R.~Brent and H.~Kung,
\emph{Systolic VLSI arrays for linear-time GCD computation},
VLSI 1983, pp.~145-154, 1983.

\bibitem{Chi1}
\begin{CJK*}{UTF8}{gbsn}九章算术\end{CJK*}, chapter 1, section 6,
1st century AD.

\bibitem{Coh1993}
H.~Cohen,
\emph{A Course in Computational Algebraic Number Theory},
Springer, pp.~29-31, 1993.

\bibitem{Euc300BC}
Euclid,
\emph{Elements}, book 7, propositions 1 and 2, c.~300~BC.

\bibitem{Falcon}
\emph{Falcon: Fast-Fourier Lattice-based Compact Signatures over NTRU},\\
\url{https://falcon-sign.info/}

\bibitem{DraftHashToCurve09}
A.~Faz-Hernandez, S.~Scott, N.~Sullivan, R.~Wahby and C.~Wood,
\emph{Hashing to Elliptic Curves},\\
\url{https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-09}

\bibitem{FogTables}
A.~Fog,
\emph{Instruction tables: Lists of instruction latencies, throughputs
and micro-operation breakdowns for Intel, AMD and VIA CPUs},\\
\url{https://www.agner.org/optimize/instruction_tables.pdf}

\bibitem{Mol2008}
N.~Möller,
\emph{On Schönhage's algorithm and subquadratic integer GCD computation},
Mathematics of Computation, vol.~77, pp.~589-607, 2008.

\bibitem{NatSar2018}
K.~Nath and P.~Sarkar,
\emph{Efficient Arithmetic In (Pseudo-)Mersenne Prime Order Fields},\\
\url{https://eprint.iacr.org/2018/985}

\bibitem{OliLopHisFazRod2017}
T.~Oliveira, J.~López, H.~Hisil, A.~Faz-Hernández and F.~Rodríguez-Henríquez,\\
\emph{How to \hbox{(pre-)compute} a ladder -- improving the performance of
X25519 and X448}, SAC 2017, Lecture Notes in Computer Science,
vol.~10719, pp.~172-191.

\bibitem{Ste1967}
J.~Stein,
\emph{Computational problems associated with Racah algebra},
Journal of Computational Physics, vol.~1, issue~3, pp.~397–405, 1967.

\end{thebibliography}

\appendix

\newpage
\section{Bounding the Number of Iterations}\label{sec:iterbound}

In this section, we prove that the inner loop of
algorithm~\ref{alg:bingcd:optimized} ensures that the $k-1$ iterations
always reduce $\bitlength(a)+\bitlength(b)$ by at least $k-1$ bits.

Each excution of the inner loop consists of $k-1$ iterations. We number
iterations from $0$ to $k-2$; this means that iteration $t$ happens just
after $t$ iterations have completed. We will note $x_t$ the value of the
quantity $x$ (for any of the quantities we manipulate, namely $a$, $b$,
$\bar a$ and $\bar b$) when \emph{exiting} iteration $t-1$ and
\emph{entering} iteration $t$.

\subsection{Maximum Error Bound}\label{sec:iterbound:error}

We first need to show that $\bar a$ and $\bar b$ remain ``good
approximations'' of $a$ and $b$ over the $k-1$ iterations, i.e.,
that the upper bits of $\bar a$ and $\bar b$ still match or are
close to the upper bits of $a$ and $b$.

First, notice that after $t$ iterations of the inner loop
($0\leq t\leq k-1$), the following inequalities hold:
\begin{align*}
    |f_0| + |g_0| &\leq 2^t \\
    |f_1| + |g_1| &\leq 2^t
\end{align*}
They are obviously correct at the start ($t = 0$) since, at that point,
$f_0 = g_1 = 1$ and $f_1 = g_0 = 0$. Then, at each successive iteration:
\begin{itemize}

    \item If $\bar a$ is even, then:
    \begin{equation*}
        \begin{gathered}
            |f_0| + |g_0| \leq 2^t < 2^{t+1} \\
            |2f_1| + |2g_1| = 2(|f_1| + |g_1|) \leq 2^{t+1}
        \end{gathered}
    \end{equation*}

    \item If $\bar a$ is odd, then the conditional swap preserves the
    inequalities, and the subtraction implies that:
    \begin{align*}
        |f_0 - f_1| + |g_0 - g_1| &\leq |f_0| + |f_1| + |g_0| + |g_1| \\
                                  &\leq |f_0| + |g_0| + |f_1| + |g_1| \\
                                  &\leq 2^{t+1}
    \end{align*}

\end{itemize}

Suppose that the inner loop started with $a$ and $b$, approximated into
$\bar a$ and $\bar b$. We have $n = \text{max}(\bitlength(a),
\bitlength(b))$. If $n \leq 2k$, then $\bar a = a$ and $\bar b = b$, and
there is no approximation; the algorithm then proceeds exactly as the
classic algorithm (algorithm~\ref{alg:bingcd:classic}) and ensures a
gain of at least 1 bit per iteration. We thus consider thereafter that
$n > 2k$.

Considering the situation after $t$ iterations of the inner loop
($0\leq t\leq k-1$), we define:
\begin{align*}
    \bar a_t &= (\bar a f_0 + \bar b g_0)/2^t \\
    \bar b_t &= (\bar a f_1 + \bar b g_1)/2^t \\
    a_t &= (a f_0 + b g_0)/2^t \\
    b_t &= (a f_1 + b g_1)/2^t
\end{align*}
These are, respectively, the values of $\bar a$ and $\bar b$ after $t$
iterations, and the corresponding values of $a$ and $b$, should the
update factors be applied at that point. Note that, by construction, the
divisions by $2^t$ are exact.

Note that $\bar a$ and $\bar b$ are approximations of $a$ and $b$ in the
sense that bits $k-1$ to $2k-1$ of $\bar a$ (respectively $\bar b$) are
exactly the same as bits $n-k-1$ to $n-1$ of $a$ (respectively $b$). This
can be expressed with the following inequalities:
\begin{align*}
    |a - 2^{n-2k} \bar a| &< 2^{n-k-1} \\
    |b - 2^{n-2k} \bar b| &< 2^{n-k-1}
\end{align*}

Consider now the following:
\begin{align*}
    a_t - 2^{n-2k} \bar a_t
        &= \frac{(af_0 + bg_0) - 2^{n-2k}(\bar a f_0 + \bar b g_0)}{2^t} \\
        &= \frac{(a - 2^{n-2k} \bar a)f_0 + (b - 2^{n-2k} \bar b)g_0}{2^t}
\end{align*}
Therefore, applying the triangular inequality:
\begin{align*}
    |a_t - 2^{n-2k} \bar a_t|
        &\leq 2^{-t} (|2^{n-k-1}| |f_0| + |2^{n-k-1}| |g_0|) \\
        &\leq 2^{n-k-1-t} (|f_0| + |g_0|)
\end{align*}
Since $|f_0| + |g_0| \leq 2^t$, we then get that:
\begin{equation*}
    |a_t - 2^{n-2k} \bar a_t| \leq 2^{n-k-1}
\end{equation*}
and similary for $b$:
\begin{equation*}
    |b_t - 2^{n-2k} \bar b_t| \leq 2^{n-k-1}
\end{equation*}

This means that, throughout the $k$ iterations, the upper bits ($k-1$ to
$2k-1$) of $\bar a$ and $\bar b$ remain a close approximation of the
upper bits ($n-k-1$ to $n-1$) of the values $a$ and $b$ that they
represent, i.e. the values we would get by applying the update factors
right away.

\subsection{Length Reduction at Divergence Point}\label{sec:iterbound:diverge}

As in the previous section, consider the values $a$ and $b$ at the start
of an inner loop. The $k-1$ iterations of the inner loop will perform
halvings, swaps and subtractions based on the bits in $\bar a$ and $\bar
b$, used as representations of $a$ and $b$. The halvings and
subtractions rely on the $k-1$ low bits of $\bar a$ and $\bar b$, which
are always exact; only the conditional swaps may differ between the
classic and optimized algorithms, since, in the latter case, the
relative ordering of $\bar a$ and $\bar b$ might not match that of $a$
and $b$.

If, throughout the $k-1$ iterations, the inner loop makes the exact same
decisions that it would take with the actual updated values of $a$ and
$b$, then it follows the same steps as the classic algorithm
(algorithm~\ref{alg:bingcd:classic}) on the same inputs, and thus
obtains the same reduction in length of $a$ and $b$, i.e. at least $k-1$
bits in total. Otherwise, we call the \emph{divergence point} the
iteration $d$ ($0\leq d < k-1$) at which the inner loop of
algorithm~\ref{alg:bingcd:optimized} takes a different decision from
that would have been taken by the classic algorithm. Let us observe the
situation at that point, using the notations from
section~\ref{sec:iterbound:error}:

\begin{itemize}

    \item The algorithm has $\bar a_d$ and $\bar b_d$, which are both odd.
    Without loss of generality, we can assume that $\bar a_d \geq \bar b_d$.
    The inner loop of algorithm~\ref{alg:bingcd:optimized} decides to
    subtract $\bar b_d$ from $\bar a_d$.

    \item However, the actual values $a_d$ and $b_d$ (updated values of
    $a$ and $b$) would be such that $a_d < b_d$. The subtraction will
    then yield a negative value.

\end{itemize}

As we saw in section~\ref{sec:iterbound:error}, $\bar a_d$ and $\bar b_d$
are good approximations of $a_d$ and $b_d$, respectively. Write that:
\begin{equation*}
    2^{n-2k} (\bar a_d - \bar b_d) - (a_d - b_d)
        = (2^{n-2k} \bar a_d - a_d) + (b_d - 2^{n-2k} \bar b_d)
\end{equation*}
and thus:
\begin{equation*}
    |2^{n-2k} (\bar a_d - \bar b_d) - (a_d - b_d)|
        \leq |2^{n-2k} \bar a_d - a_d| + |b_d - 2^{n-2k} \bar b_d|
        \leq 2^{n-k}
\end{equation*}
Since $2^{n-2k}(\bar a_d - \bar b_d) \geq 0$ and $a_d - b_d < 0$, and their
difference is no greater in absolute value than $2^{n-k}$, it follows that:
\begin{equation*}
    0 \leq \bar a_d - \bar b_d < 2^k
\end{equation*}
Therefore, at the divergence point $d$, the inner loop will compute a new
value $\bar a_{d+1} = (\bar a_d - \bar b_d)/2$ which will be strictly lower
than $2^{k-1}$, hence such that all its $k+1$ upper bits will be zero.
Similarly:
\begin{equation*}
    0 > a_d - b_d \geq -2^{n-k}
\end{equation*}
which implies that $|a_{d+1}| \leq 2^{n-k-1}$. Note that $|a_{d+1}|$ may
be equal to $2^{n-k-1}$ only if $\bar a_d = \bar b_d$, in which case $\bar
a_{d+1} = 0$ and all subsequent iterations of the inner loop will decide
to divide $\bar a_t$ (and thus $a_t$) by 2, without any further
subtraction.

\subsection{Minimal Length Reduction}\label{sec:iterbound:minred}

We again consider a situation at the start of an inner loop, such that:
\begin{itemize}

    \item One of $a$ and $b$ (or both) has length at least $2k+1$ bits.

    \item The inner loop execution includes a divergence point at
    iteration $d$.

\end{itemize}
In all other inner loop cases, the execution follows the steps of the
classic binary GCD algorithm, thereby ensuring a total length reduction
of $(a,b)$ by at least $k-1$ bits.

After the divergence point, the following properties apply:
\begin{itemize}

    \item Dividing $a_t$ by 2 does not change the sign of $a_t$.

    \item If $a_t > 0$ and $b_t < 0$, then $a_t - b_t > 0$.

    \item If $a_t < 0$ and $b_t > 0$, then $a_t - b_t < 0$.

\end{itemize}
Therefore, in all subsequent steps, one of the values will be negative,
and the other one will be positive. It cannot happen that two strictly
negative values, or two nonnegative values are obtained, unless $a$
reaches zero (in which case the algorithm has converged).

Moreover, if $|a_t| \geq 2^{n-k-1}$ and $|b_t| < 2^{n-k-1}$, then:
\begin{equation*}
    \left| \frac{a_t - b_t}{2} \right| \leq \frac{|a_t| + |b_t|}{2} < |a_t|
\end{equation*}
and if $|a_t| < 2^{n-k-1}$ and $|b_t| < 2^{n-k-1}$, then:
\begin{equation*}
    \left| \frac{a_t - b_t}{2} \right| \leq \frac{|a_t| + |b_t|}{2} < 2^{n-k-1}
\end{equation*}
As remarked in section~\ref{sec:iterbound:error}, if the divergence
point yields $a_{d+1} = -2^{n-k-1}$, then all subsequent iterations will
use $\bar a_t = 0$ and none of them will try to perform a subtraction
(or a swap); otherwise, $|a_{d+1}| < 2^{n-k-1}$. It follows that, after
the point of divergence, even if an iteration fails to reduce the sum
of the lengths of $a_t$ and $b_t$, it will not allow any value $a$ or
$b$ above $2^{n-k-1}$ to grow, or any value below $2^{n-k-1}$ to regrow
beyond $n-k-1$ bits.

Using these properties, we analyze the possible cases, depending
on the initial lengths of the values.

\vspace{2ex}
\noindent\textsf{\textbf{Case 1:}} $\bitlength(a) = n$ and
$\bitlength(b) \leq n-k$

This situation means that, when starting the inner loop, $b$ was shorter
than $a$ by at least $k$ bits.

In that case, $\bar b < 2^k$ and $\bar a \geq 2^{2k - 1}$. The first
iteration will see that $\bar a > \bar b$ and there will be no
swap. In fact, after $t$ iterations with no swap, the minimal value
for $\bar a_t$ will be $(\bar a - (2^t-1)\bar b)/2^t$. For
$t \leq k-2$, we have:
\begin{align*}
    \bar a - (2^t - 1)\bar b &\geq 2^{2k-1} - (2^t - 1)(2^k - 1) \\
        &\geq 2^{2k-1} - (2^{k-2} - 1)(2^k - 1) \\
        &>    2^{2k-1} - 2^{2k-2} \\
        &>    2^{2k-2}
\end{align*}
This implies that $\bar a_t > 2^{2k-2-t} \geq 2^{k} > \bar b_t$.
Therefore, the next iteration will still perform no swap. In total,
there will be no swap over the $k-1$ iterations. It follows that, after
$k-1$ iterations, the reduction factors are such that $f_0 = 1$, $0\geq
g_0\geq -(2^{k-1}-1)$, $f_1 = 0$, and $g_1 = 2^{k-1}$. The update on $a$
will compute $(af_0 + bg_0)/2^{k-1}$; since $a > 2^{k-1} b$ and $|g_0| <
2^{k-1}$, the new value of $a$ cannot be negative. Similarly, the new
value of $b$ cannot be negative since $f_1$ and $g_1$ are nonnnegative.
As was pointed out previously, after a divergence point, values $a_t$
and $b_t$ have opposite signs. We concluded that, in this case, there is
no divergence point, and a reduction of at least $k-1$ bits is achieved,
since the optimized algorithm follows the exact same steps as the
classic algorithm.

\vspace{2ex}
\noindent\textsf{\textbf{Case 2:}} $\bitlength(a) \leq n-k$ and
$\bitlength(b) = n$

This case is similar to the previous case. It may start with some
divisions of $a$ by 2; if $s$ such steps happen before reaching an odd
$\bar a$, then we can apply the analysis of the previous case, replacing
$k-1$ with $k-1-s$: a swap occurs at iteration $s$, and after the swap,
the new $a$ (previously $b$) has length $n$ bits while the new $b$
(previously $a$) has length at most $n-k-s$ bits. The subsequent
$k-1-s$ iterations will each reduce $a$ by at least one bit each, but
there will be no further swap and no divergence point will be reached.

\vspace{2ex}
\noindent\textsf{\textbf{Case 3:}} $\bitlength(a) = n$ and
$\bitlength(b) \geq n - k + 1$, or $\bitlength(a) \geq n - k + 1$
and $\bitlength(a) = n$

Suppose that when the divergence point ($d$) is reached,
$b_d < 2^{n-k-1}$. Since the divergence point is such that $a_d < b_d$,
this implies that $a_d < 2^{n-k-1}$ as well. Therefore, the two values
have already been reduced by at least $k+3$ bits in total at this
point. Since operations beyond the divergence point won't allow values
to regrow beyond $2^{n-k-1}$ in absolute value, it follows that the
total reduction will be at least $k+3$ bits.

We now assume that when the divergence point is reached,
$\bitlength(b_d) = n - k + j$ for some $j\geq 0$. Since $a_d < b_d$
(this is a divergence point), it implies that $\bitlength(a_d) \leq n -
k + j$ as well. Since one of the starting values had length $n$ bits,
this means that at least $k - j$ bits of reduction have been achieved at
that point. Then, the new value of $a$ is $a_{d+1}$ which is such that
$0 > a_{d+1} \geq -2^{n-k-1}$. If $a_{d+1} \neq -2^{n-k-1}$, then this
represents an extra reduction of at least $j$ bits in length, hence $k -
j + j = k$ bits at least in total. As remarked above, subsequent steps
won't allow any value of more than $k-1$ bits to grow, so this reduction
level cannot be compromised; after the $k-1$ iterations, total reduction
will still be at least $k$ bits.

The remaining case is when $a_{d+1} = -2^{n-k-1}$, exactly. At that
point, we have a reduction of at least $k - 1$ bits. Moreover, since
$\bar a_{d+1} = 0$, subsequent steps, if any, will divide $a$ by 2
but leave $b$ unchanged; there will be no loss on that reduction
level. We thus obtain at least $k-1$ bits of reduction.

\vspace{2ex}
\noindent\textsf{\textbf{Conclusion:}} In all cases, we saw that
a reduction of at least $k-1$ bits is achieved, which completes the proof.

\vspace{2ex}
\noindent\textsf{\textbf{Commentary:}} It may be surprising that even
using approximations, we do not need more iterations than the classic
binary GCD, for which the worst case bound $2\bitlength(m)-1$ is already
optimal. Intuitively, the reason is that we have some extra ``hidden''
iterations: the conditional negation of $a$ and $b$ (after their
respective updates every $k-1$ iterations) is the point where the
approximation is resolved. In that sense, the optimized algorithm really
needs $k$ iterations in the worst case to achieve a $k-1$-bit reduction;
the last of these $k$ iterations masquerades as a pair of conditional
subtractions.

\end{document}

max_stars_repo_name: sipa/bingcd

max_issues_repo_name: sipa/bingcd

max_forks_repo_name: sipa/bingcd


text:
\documentclass[]{csaund_resume-openfont}
\usepackage{fancyhdr}

\pagestyle{fancy}
\fancyhf{}

\rfoot{Page \thepage \hspace{1pt}}
\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     LAST UPDATED DATE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \lastupdated September 17, 2018

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     TITLE NAME
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\namesection{Carolyn}{Saund}{ \urlstyle{same}\href{http://carolynsaund.me}{carolynsaund.me} \\
\href{mailto:carolyn.saund@gmail.com}{carolyn.saund@gmail.com} | +44.759.557.9616
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     COLUMN ONE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{minipage}[t]{0.30\textwidth}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     EDUCATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Education}

\subsection{Tufts University}
\descript{BSc in Computer Science, Brain and Cognitive Science}
\location{May 2015 | Somerville, MA}
\sectionsep

\subsection{University of Canterbury}
\descript{BS in Computer Science}
\location{June 2013 - Jan 2014 | Christchurch, NZ}
College of Engineering \\
\location{ GPA: 4.0 / 4.0 }
\sectionsep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     LINKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     COURSEWORK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Coursework}
\subsection{Undergraduate}
Advanced Topics in Machine Learning \\
Computational Models of Cognitive Science \\
Stochastic Systems \\
Linear Algebra \\
Machine Architecture and Assembly \\
Advanced Topics in Artificial Intelligence \\
Programming Languages \\
Computational Theory \\
\sectionsep

\subsection{Continuing}
Advanced Seminar in Cognitive Architectures \\
Neural Networks for Machine Learning

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     SKILLS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Skills}
\subsection{Languages}
\location{Over 5000 lines:}
Javascript \textbullet{} Python \textbullet{} C \\
C++ \textbullet{} Java  \\
\location{Over 500 lines:}
LISP \textbullet{} CSS \textbullet{} HTTP \textbullet{} R \\
\vspace{\topsep} % Hacky fix for awkward extra vertical space
\subsection{Workflow/Frameworks}
git/github \textbullet{} NodeJS  \textbullet{} React \textbullet{} SVN \textbullet{} \\
 Unix/Linux \textbullet{} Jenkins/TDD \textbullet{} CI releases \\
 Stash \textbullet{} Phabricator/arcanist \textbullet{} Gradle \\
 MongoDB \textbullet{} Electron
\sectionsep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     COLUMN TWO
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{minipage}
\hfill
\begin{minipage}[t]{0.66\textwidth}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     EXPERIENCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Experience}
\runsubsection{Cogito} \\

\descript{Platform Team Lead}
\location{Jan 2018 - August 2018 | Boston MA}
\vspace{\topsep} %
\begin{tightemize}
\item Managed Agile team of 10 engineers
\item Oversaw top-to-bottom implementation of three major user-facing features
\end{tightemize}

\descript{Senior Software Developer }
\location{August 2017 - Jan 2018 | Boston MA}
\begin{tightemize}
\item Originally Implemented and manage release pipeline.
\item Full-stack development of Role-Based User Management (RBUM) APIs.
\end{tightemize}
\sectionsep

\runsubsection{Jibo} \\
\descript{Robot SDK Engineer }
\location{May 2015 - Dec 2016 | Boston, MA}
\begin{tightemize}
\item Employee 20(/110) Designed a series of tutorials from "Getting Started" to advanced usage still used to onboard internal and external developers.
\item Architected data APIs for dev tools to manage robot sensory input from cameras and microphones.
\item Implemented frontend UI and developer interfaces to build high-level skills.
\item Won first company culture award for "outstanding cultural contributions."
\end{tightemize}
\descript{Character AI Engineer }
\location{Jan 2017 - August 2017 | Boston, MA}
%\vspace{\topsep} % Hacky fix for awkward extra vertical space
\begin{tightemize}
\item Designed and implemented robotic bi-axial Emotion System and corresponding impacts on robot behavior.
\item Single-handedly architected and maintained robot listening pipeline (main interface for interaction).
\end{tightemize}
\descript{ Software Consultant }
\location{August 2017 - Present | Boston, MA}
%\vspace{\topsep} % Hacky fix for awkward extra vertical space
\begin{tightemize}
\item Continuing work to rearchitect APIs in cloud-based architecture
\end{tightemize}
\sectionsep

\runsubsection{Spotify} \\
\descript{Software Engineering Intern }
\location{May 2014 - Sep 2014 | NY, NY}
%\vspace{\topsep} % Hacky fix for awkward extra vertical space
\begin{tightemize}
\item Full stack development of internal desktop application.
\item Developed web APIs and GUI for music annotation tools
\end{tightemize}
\sectionsep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Public Speaking and Awards
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Projects}
\runsubsection{Drosophila}
\begin{tightemize}
\item Designed and programmed arduino "badge" with  genetic algorithms for LED display to recombine when near another badge.
\end{tightemize}
\sectionsep
\runsubsection{The Neuron}
\begin{tightemize}
\item Programmed LED strips around metal sculpture to respond based on microphone input, compute audio intensity threshold, and "fire" another LED strip with enough sound intensity.
\end{tightemize}
\sectionsep

\end{minipage}

\pagebreak

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Community Service
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Leadership and Community Service}
\runsubsection{Cambridge Inventors Club}
\descript{| Volunteer}
\location{Jan 2016 - Present | Cambridge, MA}
Work with children age 8-13 on engineering fundamentals, such as iterative design, retrospectives, and cooperation
\sectionsep

\runsubsection{Boston Ultimate Disk Alliance}
\descript{| Coach}
\location{March 2016 - Sept 2016 | Jamaica Plain, MA}
First coach on new high school girls team, focused on empowering women in the Boston ultimate community
\sectionsep

\runsubsection{Girls Who Code}
\descript{| Instructor}
\location{Sept 2015 - June 2016 | Belmont, MA}
Instructed three levels of middle- and high-school aged students from introduction to basic data strcutures. Supplemented given course materials with additional lesson plans to emphasize CS fundamentals.
\sectionsep

\runsubsection{Women in Computer Science}
\descript{| President}
\location{Jan 2013 - Jan 2015 | Somerville, MA}
\begin{tightemize}
\item Coordinated logistics of club hackathons, speakers, and outreach.
\item Led executive board in creating WiCS mission and club constitution.
\end{tightemize}
\sectionsep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Other Extracurriculars and Awards
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Extracurriculars}
\runsubsection{Awards}
\descript{ }
\begin{tightemize}
\item Culture Award (First Recipient) | Jibo Inc. Award for outstanding cultural contributions to the company.
\item Girl Scout Gold Award | Highest Honor of Girl Scouts of America
\end{tightemize}
\sectionsep

\runsubsection{Public Speaking Appearances}
\descript {}
\begin{tightemize}
\item SXSW, Robot Petting Zoo (presentation \& booth)
\item Women Engineers who Code Conference, (technical workshop \& imposter syndrome seminar)
\item Lesbians Who Tech, presentation on emotional future of robotics
\item Boston ML Meetup, presenter on Innovations in Human-Computer Interfaces
\item MLHack Championships, (technical presentation \& hackathon judge)
\item Tufts WICS Conference, Alumni Pannel
\end{tightemize}
\sectionsep

\runsubsection{Memberships \& Activities}
\descript{ }
\location{Current}
\begin{tightemize}
\item Cognitive Science Society
\item Project Management and Leadership Book Group
\item AI and Technological Futurism Book Group (Organizer)
\item Div II Ultimate Frisbee (Captain)
\end{tightemize}
\location{Former (Undergraduate)}
\begin{tightemize}
\item Vice President of Academic Development, Alpha Omicron Pi
\item Tufts Mountain Club, Trip Coordinator
\item Midnight Judges Jazz Combo
\end{tightemize}
\location{Misc.}
\begin{tightemize}
\item 20+ countries visited
\item WWOOF Worker (Germany)
\item Avid century bike-rider and marathoner
\end{tightemize}
\sectionsep

\end{document}  \documentclass[]{article}