Which algorithm does Mathematica use for double-precision N[*,16] computation of the AiryAi and AiryAiPrime functions? Are there any references? Any TOMS number or libraries?